Merge pull request #1617 from vlcvboyer/windpower

Take wind into account during power estimation
This commit is contained in:
Mark Liversedge
2015-12-04 16:06:44 +00:00

View File

@@ -24,6 +24,10 @@
#include <algorithm>
#include <QVector>
#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; i<ride->dataPoints().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="<<p->kphd<<" , V="<<V<<" , m="<<M<<" , Pa="<<(p->kphd > 1 ? 1 : p->kphd*V*M);
double watts = (afCm * V * (Ka * (vw * vw) + Frg + V * CrDyn))+(p->kphd > 1 ? 1 : p->kphd*V*M);
ride->command->setPointValue(i, RideFile::watts, watts > 0 ? (watts > 1000 ? 1000 : watts) : 0);
//qDebug()<<"w="<<p->watts<<", Ka="<<Ka<<", CwaRi="<<CwaRider<<", slope="<<p->slope<<", v="<<p->kph<<" Cwa="<<(CwaRider + CwaBike);
// qDebug() << "watts = "<<p->watts;
// qDebug() << " " << afCm * V * Ka * (vw * vw) << " = afCm(=" << afCm << ") * V(=" << V << ") * Ka(="<<Ka<<") * (vw^2(=" << V+W << "^2))";
// qDebug() << " " << afCm * V * Frg << " = afCm * V * Frg(=" << Frg << ")";
// qDebug() << " " << afCm * V * V * CrDyn << " = afCm * V^2 * CrDyn(=" << CrDyn << ")";
// qDebug() << " " << p->kphd*V*M << " = kphd(=" << p->kphd << ") * V * M(=" << M << ")";
// qDebug() << " Ka="<<Ka<<", CwaRi="<<CwaRider<<", slope="<<p->slope<<", v="<<p->kph<<" Cwa="<<(CwaRider + CwaBike);
} else {
ride->command->setPointValue(i, RideFile::watts, 0);
}