From fc223b3b77cb8e0b1c367c9d31382d40bc26a5bb Mon Sep 17 00:00:00 2001 From: Vianney Boyer Date: Fri, 4 Dec 2015 08:46:03 +0100 Subject: [PATCH] take headwind into account for power estimation --- src/FixDerivePower.cpp | 94 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 87 insertions(+), 7 deletions(-) diff --git a/src/FixDerivePower.cpp b/src/FixDerivePower.cpp index fa77c2763..6fe37ed65 100644 --- a/src/FixDerivePower.cpp +++ b/src/FixDerivePower.cpp @@ -24,6 +24,10 @@ #include #include +#ifndef MATHCONST_PI +#define MATHCONST_PI 3.141592653589793238462643383279502884L /* pi */ +#endif + // Config widget used by the Preferences/Options config panes class FixDerivePower; class FixDerivePowerConfig : public DataProcessorConfig @@ -32,11 +36,17 @@ class FixDerivePowerConfig : public DataProcessorConfig friend class ::FixDerivePower; protected: + QVBoxLayout *layoutV; QHBoxLayout *layout; + QHBoxLayout *windLayout; QLabel *bikeWeightLabel; QDoubleSpinBox *bikeWeight; QLabel *crrLabel; QDoubleSpinBox *crr; + QLabel *windSpeedLabel; + QDoubleSpinBox *windSpeed; + QLabel *windHeadingLabel; + QDoubleSpinBox *windHeading; public: FixDerivePowerConfig(QWidget *parent) : DataProcessorConfig(parent) { @@ -44,13 +54,19 @@ class FixDerivePowerConfig : public DataProcessorConfig HelpWhatsThis *help = new HelpWhatsThis(parent); parent->setWhatsThis(help->getWhatsThisText(HelpWhatsThis::MenuBar_Edit_EstimatePowerValues)); - layout = new QHBoxLayout(this); + layoutV = new QVBoxLayout(this); + layout = new QHBoxLayout(); + windLayout = new QHBoxLayout(); layout->setContentsMargins(0,0,0,0); + layoutV->setContentsMargins(0,0,0,0); + windLayout->setContentsMargins(0,0,0,0); setContentsMargins(0,0,0,0); bikeWeightLabel = new QLabel(tr("Bike Weight (kg)")); crrLabel = new QLabel(tr("Crr")); + windSpeedLabel = new QLabel(tr("Wind (kph)")); + windHeadingLabel = new QLabel(tr(", heading")); bikeWeight = new QDoubleSpinBox(); bikeWeight->setMaximum(99.9); @@ -64,11 +80,33 @@ class FixDerivePowerConfig : public DataProcessorConfig crr->setSingleStep(0.0001); crr->setDecimals(4); + windSpeed = new QDoubleSpinBox(); + windSpeed->setMaximum(99.9); + windSpeed->setMinimum(0); + windSpeed->setSingleStep(0.1); + windSpeed->setDecimals(1); + + windHeading = new QDoubleSpinBox(); + windHeading->setMaximum(180.0); + windHeading->setMinimum(-179.0); + windHeading->setSingleStep(1); + windHeading->setDecimals(0); + layout->addWidget(bikeWeightLabel); layout->addWidget(bikeWeight); layout->addWidget(crrLabel); layout->addWidget(crr); + windLayout->addWidget(windSpeedLabel); + windLayout->addWidget(windSpeed); + windLayout->addWidget(windHeadingLabel); + windLayout->addWidget(windHeading); + + layoutV->addLayout(layout); + layoutV->addLayout(windLayout); + layout->addStretch(); + layoutV->addStretch(); + windLayout->addStretch(); } //~FixDerivePowerConfig() {} // deliberately not declared since Qt will delete @@ -81,7 +119,12 @@ class FixDerivePowerConfig : public DataProcessorConfig "weight to compound total mass, it should " "include apparel, shoes, etc\n\n" "CRR parameter is the coefficient of rolling " - "resistance, it depends on tires and surface"))); + "resistance, it depends on tires and surface\n" + "wind speed shall be indicated in kph\n" + "wind heading (origin) unit is degrees " + "from -179 to +180 (-90=W, 0=N, 90=E, 180=S)\n" + "Note: if the ride file already contain wind data\n" + " it will be overridden if wind is entered manually"))); } void readConfig() { @@ -89,6 +132,8 @@ class FixDerivePowerConfig : public DataProcessorConfig bikeWeight->setValue(MBik); double Crr = appsettings->value(NULL, GC_DPDP_CRR, "0.0031").toDouble(); crr->setValue(Crr); + windSpeed->setValue(0.0); + windHeading->setValue(0.0); } void saveConfig() { @@ -130,12 +175,18 @@ FixDerivePower::postProcess(RideFile *ride, DataProcessorConfig *config=0) // get settings double MBik; // Bike weight kg double CrV; // Coefficient of Rolling Resistance + double windSpeed; // wind speed + double windHeading; //wind direction if (config == NULL) { // being called automatically MBik = appsettings->value(NULL, GC_DPDP_BIKEWEIGHT, "9.5").toDouble(); CrV = appsettings->value(NULL, GC_DPDP_CRR, "0.0031").toDouble(); + windSpeed = 0.0; + windHeading = 0.0; } else { // being called manually MBik = ((FixDerivePowerConfig*)(config))->bikeWeight->value(); CrV = ((FixDerivePowerConfig*)(config))->crr->value(); + windSpeed = ((FixDerivePowerConfig*)(config))->windSpeed->value(); // kph + windHeading = ((FixDerivePowerConfig*)(config))->windHeading->value() / 180 * MATHCONST_PI; // rad } // if its already there do nothing ! @@ -144,11 +195,12 @@ FixDerivePower::postProcess(RideFile *ride, DataProcessorConfig *config=0) // no dice if we don't have alt and speed if (!ride->areDataPresent()->alt || !ride->areDataPresent()->kph) return false; - // Power Estimation Constants + // Power Estimation Constants (mostly constants...) double hRider = ride->getHeight(); //Height in m double M = ride->getWeight(); //Weight kg double T = 15; //Temp degC in not in ride data - double W = 0; //Assume no wind + double W = 0; // headwind (from records or based on wind parameters entered manually) + double bearing = 0.0; //cyclist direction used to compute headwind double cCad=.002; double afCd = 0.62; double afSin = 0.89; @@ -170,25 +222,53 @@ FixDerivePower::postProcess(RideFile *ride, DataProcessorConfig *config=0) && ride->areDataPresent()->km) { for (int i=0; idataPoints().count(); i++) { RideFilePoint *p = ride->dataPoints()[i]; + + // compute bearing in order to calculate headwind + if (i>=1) + { + RideFilePoint *prevPoint = ride->dataPoints()[i-1]; + + // ensure a movement occurred and valid lat/lon in order to compute cyclist direction + if ( (prevPoint->lat != p->lat || prevPoint->lon != p->lon ) + && (prevPoint->lat != 0 || prevPoint->lon != 0 ) + && (p->lat != 0 || p->lon != 0 ) ) + bearing = atan2(cos(p->lat)*sin(p->lon - prevPoint->lon), + cos(prevPoint->lat)*sin(p->lat)-sin(prevPoint->lat)*cos(p->lat)*cos(p->lon - prevPoint->lon)); + } + // else keep previous bearing (or 0 at the beginning) + + // wind parameters to be considered + if (windSpeed || windHeading) // if wind parameters were entered manually then use it and override rideFile headwind + W = cos(bearing - windHeading) * windSpeed * 0.27777777777778; // Absolute wind speed relative to cyclist orientation (m/s) + else if (ride->areDataPresent()->headwind) // otherwise use headwind from rideFile records (typ. weather forecast included in FIT files) + W = (p->headwind - p->kph) * 0.27777777777778; // Compute absolute wind speed relative to cyclist orientation (m/s) + else + W = 0.0; //otherwise assume no wind + // Estimate Power if not in data double cad = ride->areDataPresent()->cad ? p->cad : 85.00; if (cad > 0) { if (ride->areDataPresent()->temp) T = p->temp; double Slope = atan(p->slope * .01); - double V = p->kph * 0.27777777777778; //Speed m/s + double V = p->kph * 0.27777777777778; // Cyclist speed m/s double CrDyn = 0.1 * cos(Slope); double CwaRider, Ka; double Frg = 9.81 * (MBik + M) * (CrEff * cos(Slope) + sin(Slope)); - double vw=V+W; + double vw=V+W; // Wind speed against cyclist = cyclist speed + wind speed CwaRider = (1 + cad * cCad) * afCd * adipos * (((hRider - adipos) * afSin) + adipos); Ka = 176.5 * exp(-p->alt * .0001253) * (CwaRider + CwaBike) / (273 + T); //qDebug()<<"acc="<kphd<<" , V="<kphd*V*M << " = kphd(=" << p->kphd << ") * V * M(=" << M << ")"; + // qDebug() << " Ka="<