/* * Copyright (c) 2009 Andy M. Froncioni (me@andyfroncioni.com) * * This program is free software; you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by the Free * Software Foundation; either version 2 of the License, or (at your option) * any later version. * * This program is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for * more details. * * You should have received a copy of the GNU General Public License along * with this program; if not, write to the Free Software Foundation, Inc., 51 * Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ #include "Aerolab.h" #include "AerolabWindow.h" #include "Context.h" #include "Athlete.h" #include "IntervalItem.h" #include "RideFile.h" #include "RideItem.h" #include "Settings.h" #include "Units.h" #include "Colors.h" #include #include #include #include #include #include #include #include #include #define PI M_PI static inline double max(double a, double b) { if (a > b) return a; else return b; } static inline double min(double a, double b) { if (a < b) return a; else return b; } /*---------------------------------------------------------------------- * Interval plotting *--------------------------------------------------------------------*/ class IntervalAerolabData : public QwtSeriesData { public: Aerolab *aerolab; Context *context; IntervalAerolabData ( Aerolab *aerolab, Context *context ) : aerolab( aerolab ), context(context) { } double x( size_t ) const; double y( size_t ) const; size_t size() const; //virtual QwtData *copy() const; void init(); IntervalItem *intervalNum( int ) const; int intervalCount() const; virtual QPointF sample(size_t i) const; virtual QRectF boundingRect() const; }; /* * HELPER FUNCTIONS: * intervalNum - returns a pointer to the nth selected interval * intervalCount - returns the number of highlighted intervals */ // ------------------------------------------------------------------------------------------------------------ // note this is operating on the children of allIntervals and not the // intervalWidget (QTreeWidget) -- this is why we do not use the // selectedItems() member. N starts a one not zero. // ------------------------------------------------------------------------------------------------------------ IntervalItem *IntervalAerolabData::intervalNum ( int number ) const { int highlighted = 0; const QTreeWidgetItem *allIntervals = context->athlete->allIntervalItems(); for ( int ii = 0; ii < allIntervals->childCount(); ii++) { IntervalItem *current = (IntervalItem *) allIntervals->child( ii ); if ( current == NULL) { return NULL; } if ( current->isSelected() == true ) { ++highlighted; } if ( highlighted == number ) { return current; } } return NULL; } // ------------------------------------------------------------------------------------------------------------ // how many intervals selected? // ------------------------------------------------------------------------------------------------------------ int IntervalAerolabData::intervalCount() const { int highlighted = 0; if ( context->athlete->allIntervalItems() != NULL ) { const QTreeWidgetItem *allIntervals = context->athlete->allIntervalItems(); for ( int ii = 0; ii < allIntervals->childCount(); ii++) { IntervalItem *current = (IntervalItem *) allIntervals->child( ii ); if ( current != NULL ) { if ( current->isSelected() == true ) { ++highlighted; } } } } return highlighted; } /* * INTERVAL HIGHLIGHTING CURVE * IntervalAerolabData - implements the qwtdata interface where * x,y return point co-ordinates and * size returns the number of points */ // The interval curve data is derived from the intervals that have // been selected in the Context leftlayout for each selected // interval we return 4 data points; bottomleft, topleft, topright // and bottom right. // // the points correspond to: // bottom left = interval start, 0 watts // top left = interval start, maxwatts // top right = interval stop, maxwatts // bottom right = interval stop, 0 watts // double IntervalAerolabData::x ( size_t number ) const { // for each interval there are four points, which interval is this for? // interval numbers start at 1 not ZERO in the utility functions double result = 0; int interval_no = number ? 1 + number / 4 : 1; // get the interval IntervalItem *current = intervalNum( interval_no ); if ( current != NULL ) { double multiplier = context->athlete->useMetricUnits ? 1 : MILES_PER_KM; // which point are we returning? //qDebug() << "number = " << number << endl; switch ( number % 4 ) { case 0 : result = aerolab->bydist ? multiplier * current->startKM : current->start/60; // bottom left break; case 1 : result = aerolab->bydist ? multiplier * current->startKM : current->start/60; // top left break; case 2 : result = aerolab->bydist ? multiplier * current->stopKM : current->stop/60; // bottom right break; case 3 : result = aerolab->bydist ? multiplier * current->stopKM : current->stop/60; // top right break; } } return result; } double IntervalAerolabData::y ( size_t number ) const { // which point are we returning? double result = 0; switch ( number % 4 ) { case 0 : result = -5000; // bottom left break; case 1 : result = 5000; // top left - set to out of bound value break; case 2 : result = 5000; // top right - set to out of bound value break; case 3 : result = -5000; // bottom right break; } return result; } size_t IntervalAerolabData::size() const { return intervalCount() * 4; } QPointF IntervalAerolabData::sample(size_t i) const { return QPointF(x(i), y(i)); } QRectF IntervalAerolabData::boundingRect() const { return QRectF(-5000, 5000, 10000, 10000); } //********************************************** //** END IntervalAerolabData ** //********************************************** Aerolab::Aerolab( AerolabWindow *parent, Context *context ): QwtPlot(parent), context(context), parent(parent), rideItem(NULL), smooth(1), bydist(true), autoEoffset(true) { crr = 0.005; cda = 0.500; totalMass = 85.0; rho = 1.236; eta = 1.0; eoffset = 0.0; constantAlt = false; insertLegend(new QwtLegend(), QwtPlot::BottomLegend); setCanvasBackground(Qt::white); static_cast(canvas())->setFrameStyle(QFrame::NoFrame); setXTitle(); setAxisTitle(yLeft, tr("Elevation (m)")); setAxisScale(yLeft, -300, 300); setAxisTitle(xBottom, tr("Distance (km)")); setAxisScale(xBottom, 0, 60); veCurve = new QwtPlotCurve(tr("V-Elevation")); altCurve = new QwtPlotCurve(tr("Elevation")); // get rid of nasty blank space on right of the plot veCurve->setYAxis( yLeft ); altCurve->setYAxis( yLeft ); intervalHighlighterCurve = new QwtPlotCurve(); intervalHighlighterCurve->setBaseline(-5000); intervalHighlighterCurve->setYAxis( yLeft ); intervalHighlighterCurve->setSamples( new IntervalAerolabData( this, context ) ); intervalHighlighterCurve->attach( this ); //XXX broken this->legend()->remove( intervalHighlighterCurve ); // don't show in legend grid = new QwtPlotGrid(); grid->enableX(false); grid->attach(this); configChanged(CONFIG_APPEARANCE); } void Aerolab::configChanged(qint32) { // set colors setCanvasBackground(GColor(CPLOTBACKGROUND)); QPen vePen = QPen(GColor(CAEROVE)); vePen.setWidth(appsettings->value(this, GC_LINEWIDTH, 0.5).toDouble()); veCurve->setPen(vePen); QPen altPen = QPen(GColor(CAEROEL)); altPen.setWidth(appsettings->value(this, GC_LINEWIDTH, 0.5).toDouble()); altCurve->setPen(altPen); QPen gridPen(GColor(CPLOTGRID)); gridPen.setStyle(Qt::DotLine); grid->setPen(gridPen); QPen ihlPen = QPen( GColor( CINTERVALHIGHLIGHTER ) ); ihlPen.setWidth(1); intervalHighlighterCurve->setPen( ihlPen ); QColor ihlbrush = QColor(GColor(CINTERVALHIGHLIGHTER)); ihlbrush.setAlpha(40); intervalHighlighterCurve->setBrush(ihlbrush); // fill below the line //XXX broken this->legend()->remove( intervalHighlighterCurve ); // don't show in legend } void Aerolab::setData(RideItem *_rideItem, bool new_zoom) { // HARD-CODED DATA: p1->kph double vfactor = 3.600; double m = totalMass; double small_number = 0.00001; rideItem = _rideItem; RideFile *ride = rideItem->ride(); veArray.clear(); altArray.clear(); distanceArray.clear(); timeArray.clear(); if( ride ) { const RideFileDataPresent *dataPresent = ride->areDataPresent(); //setTitle(ride->startTime().toString(GC_DATETIME_FORMAT)); if( dataPresent->watts ) { // If watts are present, then we can fill the veArray data: const RideFileDataPresent *dataPresent = ride->areDataPresent(); int npoints = ride->dataPoints().size(); double dt = ride->recIntSecs(); veArray.resize(dataPresent->watts ? npoints : 0); altArray.resize(dataPresent->alt || constantAlt ? npoints : 0); timeArray.resize(dataPresent->watts ? npoints : 0); distanceArray.resize(dataPresent->watts ? npoints : 0); // quickly erase old data veCurve->setVisible(false); altCurve->setVisible(false); // detach and re-attach the ve curve: veCurve->detach(); if (!veArray.empty()) { veCurve->attach(this); veCurve->setVisible(dataPresent->watts); } // detach and re-attach the ve curve: bool have_recorded_alt_curve = false; altCurve->detach(); if (!altArray.empty()) { have_recorded_alt_curve = true; altCurve->attach(this); altCurve->setVisible(dataPresent->alt || constantAlt ); } // Fill the virtual elevation profile with data from the ride data: double t = 0.0; double vlast = 0.0; double e = 0.0; arrayLength = 0; foreach(const RideFilePoint *p1, ride->dataPoints()) { if ( arrayLength == 0 ) e = eoffset; timeArray[arrayLength] = p1->secs / 60.0; if ( have_recorded_alt_curve ) { if ( constantAlt && arrayLength > 0) { altArray[arrayLength] = altArray[arrayLength-1]; } else { if ( constantAlt && !dataPresent->alt) altArray[arrayLength] = 0; else altArray[arrayLength] = (context->athlete->useMetricUnits ? p1->alt : p1->alt * FEET_PER_METER); } } // Unpack: double power = max(0, p1->watts); double v = p1->kph/vfactor; double headwind = v; if( dataPresent->headwind ) { headwind = p1->headwind/vfactor; } double f = 0.0; double a = 0.0; // Use km data instead of formula for file with a stop (gap). //d += v * dt; //distanceArray[arrayLength] = d/1000; distanceArray[arrayLength] = p1->km; if( v > small_number ) { f = power/v; a = ( v*v - vlast*vlast ) / ( 2.0 * dt * v ); } else { a = ( v - vlast ) / dt; } f *= eta; // adjust for drivetrain efficiency if using a crank-based meter double s = slope( f, a, m, crr, cda, rho, headwind ); double de = s * v * dt * (context->athlete->useMetricUnits ? 1 : FEET_PER_METER); e += de; t += dt; veArray[arrayLength] = e; vlast = v; ++arrayLength; } } else { veCurve->setVisible(false); altCurve->setVisible(false); } recalc(new_zoom); adjustEoffset(); } else { //setTitle("no data"); } } void Aerolab::setAxisTitle(int axis, QString label) { // setup the default fonts QFont stGiles; // hoho - Chart Font St. Giles ... ok you have to be British to get this joke stGiles.fromString(appsettings->value(this, GC_FONT_CHARTLABELS, QFont().toString()).toString()); stGiles.setPointSize(appsettings->value(NULL, GC_FONT_CHARTLABELS_SIZE, 8).toInt()); QwtText title(label); title.setFont(stGiles); QwtPlot::setAxisFont(axis, stGiles); QwtPlot::setAxisTitle(axis, title); } void Aerolab::adjustEoffset() { if (autoEoffset && !altArray.empty()) { double idx = axisScaleDiv( QwtPlot::xBottom ).lowerBound(); parent->eoffsetSlider->setEnabled(false); if (bydist) { int v = 100*(altArray.at(rideItem->ride()->distanceIndex(idx))-veArray.at(rideItem->ride()->distanceIndex(idx))); parent->eoffsetSlider->setValue(intEoffset()+v); } else { int v = 100*(altArray.at(rideItem->ride()->timeIndex(60*idx))-veArray.at(rideItem->ride()->timeIndex(60*idx))); parent->eoffsetSlider->setValue(intEoffset()+v); } } else parent->eoffsetSlider->setEnabled(true); } struct DataPoint { double time, hr, watts, speed, cad, alt; DataPoint(double t, double h, double w, double s, double c, double a) : time(t), hr(h), watts(w), speed(s), cad(c), alt(a) {} }; void Aerolab::recalc( bool new_zoom ) { if (timeArray.empty()) return; int rideTimeSecs = (int) ceil(timeArray[arrayLength - 1]); int totalRideDistance = (int ) ceil(distanceArray[arrayLength - 1]); // If the ride is really long, then avoid it like the plague. if (rideTimeSecs > 7*24*60*60) { QVector data; if (!veArray.empty()){ veCurve->setSamples(data, data); } if( !altArray.empty()) { altCurve->setSamples(data, data); } return; } QVector &xaxis = (bydist?distanceArray:timeArray); int startingIndex = 0; int totalPoints = arrayLength - startingIndex; // set curves if (!veArray.empty()) { veCurve->setSamples(xaxis.data() + startingIndex, veArray.data() + startingIndex, totalPoints); } if (!altArray.empty()){ altCurve->setSamples(xaxis.data() + startingIndex, altArray.data() + startingIndex, totalPoints); } if( new_zoom ) setAxisScale(xBottom, 0.0, (bydist?totalRideDistance:rideTimeSecs)); setYMax(new_zoom ); refreshIntervalMarkers(); replot(); } void Aerolab::setYMax(bool new_zoom) { if (veCurve->isVisible()) { if ( context->athlete->useMetricUnits ) { setAxisTitle( yLeft, tr("Elevation (m)") ); } else { setAxisTitle( yLeft, tr("Elevation (')") ); } double minY = 0.0; double maxY = 0.0; //************ //if (veCurve->isVisible()) { // setAxisTitle(yLeft, tr("Elevation")); if ( !altArray.empty() ) { // setAxisScale(yLeft, // min( veCurve->minYValue(), altCurve->minYValue() ) - 10, // 10.0 + max( veCurve->maxYValue(), altCurve->maxYValue() ) ); minY = min( veCurve->minYValue(), altCurve->minYValue() ) - 10; maxY = 10.0 + max( veCurve->maxYValue(), altCurve->maxYValue() ); } else { //setAxisScale(yLeft, // veCurve->minYValue() , // 1.05 * veCurve->maxYValue() ); if ( new_zoom ) { minY = veCurve->minYValue(); maxY = veCurve->maxYValue(); } else { minY = parent->getCanvasTop(); maxY = parent->getCanvasBottom(); } //adjust eooffset // TODO } setAxisScale( yLeft, minY, maxY ); setAxisLabelRotation(yLeft,270); setAxisLabelAlignment(yLeft,Qt::AlignVCenter); } enableAxis(yLeft, veCurve->isVisible()); } void Aerolab::setXTitle() { if (bydist) setAxisTitle(xBottom, tr("Distance ")+QString(context->athlete->useMetricUnits?"(km)":"(miles)")); else setAxisTitle(xBottom, tr("Time (minutes)")); } void Aerolab::setAutoEoffset(int value) { autoEoffset = value; adjustEoffset(); } void Aerolab::setConstantAlt(int value) { constantAlt = value; } void Aerolab::setByDistance(int value) { bydist = value; setXTitle(); recalc(true); } double Aerolab::slope( double f, double a, double m, double crr, double cda, double rho, double v ) { double g = 9.80665; // Small angle version of slope calculation: double s = f/(m*g) - crr - cda*rho*v*v/(2.0*m*g) - a/g; return s; } // At slider 1000, we want to get max Crr=0.1000 // At slider 1 , we want to get min Crr=0.0001 void Aerolab::setIntCrr( int value ) { crr = (double) value / 1000000.0; recalc(false); } // At slider 1000, we want to get max CdA=1.000 // At slider 1 , we want to get min CdA=0.001 void Aerolab::setIntCda( int value ) { cda = (double) value / 10000.0; recalc(false); } // At slider 1000, we want to get max CdA=1.000 // At slider 1 , we want to get min CdA=0.001 void Aerolab::setIntTotalMass( int value ) { totalMass = (double) value / 100.0; recalc(false); } // At slider 1000, we want to get max CdA=1.000 // At slider 1 , we want to get min CdA=0.001 void Aerolab::setIntRho( int value ) { rho = (double) value / 10000.0; recalc(false); } // At slider 1000, we want to get max CdA=1.000 // At slider 1 , we want to get min CdA=0.001 void Aerolab::setIntEta( int value ) { eta = (double) value / 10000.0; recalc(false); } // At slider 1000, we want to get max CdA=1.000 // At slider 1 , we want to get min CdA=0.001 void Aerolab::setIntEoffset( int value ) { eoffset = (double) value / 100.0; recalc(false); } void Aerolab::pointHover (QwtPlotCurve *curve, int index) { if ( index >= 0 && curve != intervalHighlighterCurve ) { double x_value = curve->sample( index ).x(); double y_value = curve->sample( index ).y(); // output the tooltip QString text = QString( "%1 %2 %3 %4 %5" ) . arg( this->axisTitle( curve->xAxis() ).text() ) . arg( x_value, 0, 'f', 3 ) . arg( "\n" ) . arg( this->axisTitle( curve->yAxis() ).text() ) . arg( y_value, 0, 'f', 3 ); // set that text up tooltip->setText( text ); } else { // no point tooltip->setText( "" ); } } void Aerolab::refreshIntervalMarkers() { foreach( QwtPlotMarker *mrk, d_mrk ) { mrk->detach(); delete mrk; } d_mrk.clear(); QRegExp wkoAuto("^(Peak *[0-9]*(s|min)|Entire workout|Find #[0-9]*) *\\([^)]*\\)$"); if ( rideItem->ride() ) { foreach(const RideFileInterval &interval, rideItem->ride()->intervals()) { // skip WKO autogenerated peak intervals if (wkoAuto.exactMatch(interval.name)) continue; QwtPlotMarker *mrk = new QwtPlotMarker; d_mrk.append(mrk); mrk->attach(this); mrk->setLineStyle(QwtPlotMarker::VLine); mrk->setLabelAlignment(Qt::AlignRight | Qt::AlignTop); mrk->setLinePen(QPen(GColor(CPLOTMARKER), 0, Qt::DashDotLine)); QwtText text(interval.name); text.setFont(QFont("Helvetica", 10, QFont::Bold)); text.setColor(GColor(CPLOTMARKER)); if (!bydist) mrk->setValue(interval.start / 60.0, 0.0); else mrk->setValue((context->athlete->useMetricUnits ? 1 : MILES_PER_KM) * rideItem->ride()->timeToDistance(interval.start), 0.0); mrk->setLabel(text); } } } /* * Estimate CdA and Crr usign energy balance in segments defined by * non-zero altitude. * Returns an explanatory error message ff it fails to do the estimation, * otherwise it updates cda and crr and returns an empty error message. * Author: Alejandro Martinez * Date: 23-aug-2012 */ QString Aerolab::estimateCdACrr(RideItem *rideItem) { // HARD-CODED DATA: p1->kph const double vfactor = 3.600; const double g = 9.80665; RideFile *ride = rideItem->ride(); QString errMsg; if(ride) { const RideFileDataPresent *dataPresent = ride->areDataPresent(); if(( dataPresent->alt || constantAlt ) && dataPresent->watts) { double dt = ride->recIntSecs(); int npoints = ride->dataPoints().size(); double X1[npoints], X2[npoints], Egain[npoints]; int nSeg = -1; double altInit = 0, vInit = 0; /* For each segment, defined between points with alt != 0, * this loop computes X1, X2 and Egain to verify: * Aero-Loss + RR-Loss = Egain * where * Aero-Loss = X1[nSgeg] * CdA * RR-Loss = X2[nSgeg] * Crr * are the aero and rr components of the energy loss with * X1[nSeg] = sum(0.5 * rho * headwind*headwind * distance) * X2[nSeg] = sum(totalMass * g * distance) * and the energy gain sums power in the segment with * potential and kinetic variations: * Egain = sum(eta * power * dt) + * totalMass * (g * (altInit - alt) + * 0.5 * (vInit*vInit - v*v)) */ foreach(const RideFilePoint *p1, ride->dataPoints()) { // Unpack: double power = max(0, p1->watts); double v = p1->kph/vfactor; double distance = v * dt; double headwind = v; if( dataPresent->headwind ) { headwind = p1->headwind/vfactor; } double alt = p1->alt; // start initial segment if (nSeg < 0 && alt != 0) { nSeg = 0; X1[nSeg] = X2[nSeg] = Egain[nSeg] = 0.0; altInit = alt; vInit = v; } // accumulate segment data if (nSeg >= 0) { // X1[nSgeg] * CdA == Aero-Loss X1[nSeg] += 0.5 * rho * headwind*headwind * distance; // X2[nSgeg] * Crr == RR-Loss X2[nSeg] += totalMass * g * distance; // Energy supplied Egain[nSeg] += eta * power * dt; } // close current segment and start a new one if (nSeg >= 0 && alt != 0) { // Add change in potential and kinetic energy Egain[nSeg] += totalMass * (g * (altInit - alt) + 0.5 * (vInit*vInit - v*v)); // Start a new segment nSeg++; X1[nSeg] = X2[nSeg] = Egain[nSeg] = 0.0; altInit = alt; vInit = v; } } /* At least two segmentes needed to approximate: * X1 * CdA + X2 * Crr = Egain * which, in matrix form, is: * X * [ CdA ; Crr ] = Egain * and pre-multiplying by X transpose (X'): * X'* X [ CdA ; Crr ] = X' * Egain * which is a system with two equations and two unknowns * A * [ CdA ; Crr ] = B */ if (nSeg >= 2) { // compute the normal equation: A * [ CdA ; Crr ] = B // A = X'*X // B = X'*Egain double A11 = 0, A12 = 0, A21 = 0, A22 = 0, B1 = 0, B2 = 0; for (int i = 0; i < nSeg; i++) { A11 += X1[i] * X1[i]; A12 += X1[i] * X2[i]; A21 += X2[i] * X1[i]; A22 += X2[i] * X2[i]; B1 += X1[i] * Egain[i]; B2 += X2[i] * Egain[i]; } // Solve the normal equation // A11 * CdA + A12 * Crr = B1 // A21 * CdA + A22 * Crr = B2 double det = A11 * A22 - A12 * A21; if (fabs(det) > 0.00) { // round and update if the values are in Aerolab's range double cda = floor(10000 * (A22 * B1 - A12 * B2) / det + 0.5) / 10000; double crr = floor(1000000 * (A11 * B2 - A21 * B1) / det + 0.5) / 1000000; if (cda >= 0.001 and cda <= 1.0 and crr >= 0.0001 and crr <= 0.1) { this->cda = cda; this->crr = crr; errMsg = ""; // No error } else { errMsg = tr("Estimates out-of-range"); } } else { errMsg = tr("At least two segments must be independent"); } } else { errMsg = tr("At least two segments must be defined"); } } else { errMsg = tr("Altitude and Power data must be present"); } } else { errMsg = tr("No activity selected"); } return errMsg; }