← Back to branch summary

~ubuntu-branches/ubuntu/saucy/goldencheetah/saucy

« back to all changes in this revision

Viewing changes to src/Aerolab.cpp

  • Committer: Package Import Robot
  • Author(s): KURASHIKI Satoru
  • Date: 2013-08-18 07:02:45 UTC
  • mfrom: (4.1.8 sid)
  • Revision ID: package-import@ubuntu.com-20130818070245-zgdvb47e1k3mtgil
Tags: 3.0-3
http://bugs.debian.org/719571
debian/control: remove needless dependency. (Closes: #719571)

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Aerolab.cpp
19
19
#include "Aerolab.h"
20
20
#include "AerolabWindow.h"
21
21
#include "MainWindow.h"
 
22
#include "IntervalItem.h"
22
23
#include "RideFile.h"
23
24
#include "RideItem.h"
24
25
#include "Settings.h"
27
28
 
28
29
#include <math.h>
29
30
#include <assert.h>
30
 
#include <qwt_data.h>
 
31
#include <qwt_series_data.h>
31
32
#include <qwt_legend.h>
32
33
#include <qwt_plot_curve.h>
33
34
#include <qwt_plot_grid.h>
43
44
static inline double
44
45
min(double a, double b) { if (a < b) return a; else return b; }
45
46
 
46
 
 
47
47
/*----------------------------------------------------------------------
48
48
 * Interval plotting
49
49
 *--------------------------------------------------------------------*/
50
50
 
51
 
class IntervalAerolabData : public QwtData
 
51
class IntervalAerolabData : public QwtSeriesData<QPointF>
52
52
{
53
53
    public:
54
54
        Aerolab *aerolab;
62
62
        double x( size_t ) const;
63
63
        double y( size_t ) const;
64
64
        size_t size() const;
65
 
        virtual QwtData *copy() const;
 
65
        //virtual QwtData *copy() const;
66
66
 
67
67
        void init();
68
68
 
69
69
        IntervalItem *intervalNum( int ) const;
70
70
 
71
71
        int intervalCount() const;
 
72
 
 
73
        virtual QPointF sample(size_t i) const;
 
74
        virtual QRectF boundingRect() const;
72
75
};
73
76
 
74
77
/*
210
213
    return intervalCount() * 4;
211
214
}
212
215
 
213
 
QwtData *IntervalAerolabData::copy() const
 
216
QPointF IntervalAerolabData::sample(size_t i) const {
 
217
    return QPointF(x(i), y(i));
 
218
}
 
219
 
 
220
QRectF IntervalAerolabData::boundingRect() const
214
221
{
215
 
    return new IntervalAerolabData( aerolab, mainWindow );
 
222
    return QRectF(-5000, 5000, 10000, 10000);
216
223
}
217
224
 
218
 
 
219
225
//**********************************************
220
226
//**        END IntervalAerolabData           **
221
227
//**********************************************
226
232
    MainWindow    *mainWindow
227
233
):
228
234
  QwtPlot(parent),
 
235
  mainWindow(mainWindow),
229
236
  parent(parent),
230
 
  unit(0),
231
237
  rideItem(NULL),
232
238
  smooth(1), bydist(true), autoEoffset(true) {
233
239
 
238
244
  eta       = 1.0;
239
245
  eoffset   = 0.0;
240
246
 
241
 
  boost::shared_ptr<QSettings> settings = GetApplicationSettings();
242
 
  unit = settings->value(GC_UNIT);
243
 
  useMetricUnits = true;
 
247
  useMetricUnits = mainWindow->useMetricUnits;
244
248
 
245
249
  insertLegend(new QwtLegend(), QwtPlot::BottomLegend);
246
250
  setCanvasBackground(Qt::white);
 
251
  canvas()->setFrameStyle(QFrame::NoFrame);
247
252
 
248
253
  setXTitle();
249
254
  setAxisTitle(yLeft, tr("Elevation (m)"));
261
266
  intervalHighlighterCurve = new QwtPlotCurve();
262
267
  intervalHighlighterCurve->setBaseline(-5000);
263
268
  intervalHighlighterCurve->setYAxis( yLeft );
264
 
  intervalHighlighterCurve->setData( IntervalAerolabData( this, mainWindow ) );
 
269
  intervalHighlighterCurve->setData( new IntervalAerolabData( this, mainWindow ) );
265
270
  intervalHighlighterCurve->attach( this );
266
271
  this->legend()->remove( intervalHighlighterCurve ); // don't show in legend
267
272
 
276
281
void
277
282
Aerolab::configChanged()
278
283
{
 
284
  useMetricUnits = mainWindow->useMetricUnits;
 
285
 
279
286
  // set colors
280
287
  setCanvasBackground(GColor(CPLOTBACKGROUND));
281
288
  QPen vePen = QPen(GColor(CAEROVE));
282
 
  vePen.setWidth(1);
 
289
  vePen.setWidth(appsettings->value(this, GC_LINEWIDTH, 2.0).toDouble());
283
290
  veCurve->setPen(vePen);
284
291
  QPen altPen = QPen(GColor(CAEROEL));
285
 
  altPen.setWidth(1);
 
292
  altPen.setWidth(appsettings->value(this, GC_LINEWIDTH, 2.0).toDouble());
286
293
  altCurve->setPen(altPen);
287
294
  QPen gridPen(GColor(CPLOTGRID));
288
295
  gridPen.setStyle(Qt::DotLine);
319
326
  if( ride ) {
320
327
 
321
328
    const RideFileDataPresent *dataPresent = ride->areDataPresent();
322
 
    setTitle(ride->startTime().toString(GC_DATETIME_FORMAT));
 
329
    //setTitle(ride->startTime().toString(GC_DATETIME_FORMAT));
323
330
 
324
331
    if( dataPresent->watts ) {
325
332
 
356
363
      double t = 0.0;
357
364
      double vlast = 0.0;
358
365
      double e     = 0.0;
359
 
      double d = 0;
360
366
      arrayLength = 0;
361
367
      foreach(const RideFilePoint *p1, ride->dataPoints()) {
362
368
      if ( arrayLength == 0 )
413
419
    recalc(new_zoom);
414
420
    adjustEoffset();
415
421
  } else {
416
 
    setTitle("no data");
 
422
    //setTitle("no data");
417
423
 
418
424
  }
419
425
}
420
426
 
421
427
void
 
428
Aerolab::setAxisTitle(int axis, QString label)
 
429
{
 
430
    // setup the default fonts
 
431
    QFont stGiles; // hoho - Chart Font St. Giles ... ok you have to be British to get this joke
 
432
    stGiles.fromString(appsettings->value(this, GC_FONT_CHARTLABELS, QFont().toString()).toString());
 
433
    stGiles.setPointSize(appsettings->value(NULL, GC_FONT_CHARTLABELS_SIZE, 8).toInt());
 
434
 
 
435
    QwtText title(label);
 
436
    title.setFont(stGiles);
 
437
    QwtPlot::setAxisFont(axis, stGiles);
 
438
    QwtPlot::setAxisTitle(axis, title);
 
439
}
 
440
 
 
441
void
422
442
Aerolab::adjustEoffset() {
423
443
 
424
444
    if (autoEoffset && !altArray.empty()) {
454
474
 
455
475
  // If the ride is really long, then avoid it like the plague.
456
476
  if (rideTimeSecs > 7*24*60*60) {
457
 
    QwtArray<double> data;
 
477
    QVector<double> data;
 
478
 
458
479
    if (!veArray.empty()){
459
480
      veCurve->setData(data, data);
460
481
    }
469
490
  int totalPoints   = arrayLength - startingIndex;
470
491
 
471
492
  // set curves
472
 
  if (!veArray.empty())
473
 
    veCurve->setData(xaxis.data() + startingIndex,
474
 
             veArray.data() + startingIndex, totalPoints);
 
493
  if (!veArray.empty()) {
 
494
      veCurve->setData(xaxis.data() + startingIndex, veArray.data() + startingIndex, totalPoints);
 
495
  }
475
496
 
476
497
  if (!altArray.empty()){
477
 
    altCurve->setData(xaxis.data() + startingIndex,
478
 
              altArray.data() + startingIndex, totalPoints);
 
498
      altCurve->setData(xaxis.data() + startingIndex, altArray.data() + startingIndex, totalPoints);
479
499
  }
480
500
 
481
501
  if( new_zoom )
501
521
 
502
522
             {
503
523
 
504
 
                 setAxisTitle( yLeft, "Elevation (m)" );
 
524
                 setAxisTitle( yLeft, tr("Elevation (m)") );
505
525
 
506
526
             }
507
527
 
509
529
 
510
530
             {
511
531
 
512
 
                 setAxisTitle( yLeft, "Elevation (')" );
 
532
                 setAxisTitle( yLeft, tr("Elevation (')") );
513
533
 
514
534
             }
515
535
 
574
594
Aerolab::setXTitle() {
575
595
 
576
596
  if (bydist)
577
 
    setAxisTitle(xBottom, tr("Distance ")+QString(unit.toString() == "Metric"?"(km)":"(miles)"));
 
597
    setAxisTitle(xBottom, tr("Distance ")+QString(useMetricUnits?"(km)":"(miles)"));
578
598
  else
579
599
    setAxisTitle(xBottom, tr("Time (minutes)"));
580
600
}
695
715
{
696
716
    if ( index >= 0 && curve != intervalHighlighterCurve )
697
717
    {
698
 
        double x_value = curve->x( index );
699
 
        double y_value = curve->y( index );
 
718
        double x_value = curve->sample( index ).x();
 
719
        double y_value = curve->sample( index ).y();
700
720
        // output the tooltip
701
721
 
702
722
        QString text = QString( "%1 %2 %3 %4 %5" )
750
770
        }
751
771
    }
752
772
}
 
773
 
 
774
/*
 
775
 * Estimate CdA and Crr usign energy balance in segments defined by
 
776
 * non-zero altitude.
 
777
 * Returns an explanatory error message ff it fails to do the estimation,
 
778
 * otherwise it updates cda and crr and returns an empty error message.
 
779
 * Author: Alejandro Martinez
 
780
 * Date: 23-aug-2012
 
781
 */
 
782
QString Aerolab::estimateCdACrr(RideItem *rideItem)
 
783
{
 
784
    // HARD-CODED DATA: p1->kph
 
785
    const double vfactor = 3.600;
 
786
    const double g = 9.80665;
 
787
    RideFile *ride = rideItem->ride();
 
788
    QString errMsg;
 
789
 
 
790
    if(ride) {
 
791
        const RideFileDataPresent *dataPresent = ride->areDataPresent();
 
792
        if(dataPresent->alt && dataPresent->watts) {
 
793
            double dt = ride->recIntSecs();
 
794
            int npoints = ride->dataPoints().size();
 
795
            double X1[npoints], X2[npoints], Egain[npoints];
 
796
            int nSeg = -1;
 
797
            double altInit = 0, vInit = 0;
 
798
            /* For each segment, defined between points with alt != 0,
 
799
             * this loop computes X1, X2 and Egain to verify:
 
800
             * Aero-Loss + RR-Loss = Egain
 
801
             * where
 
802
             *      Aero-Loss = X1[nSgeg] * CdA
 
803
             *      RR-Loss = X2[nSgeg] * Crr
 
804
             * are the aero and rr components of the energy loss with
 
805
             *      X1[nSeg] = sum(0.5 * rho * headwind*headwind * distance)
 
806
             *      X2[nSeg] = sum(totalMass * g * distance)
 
807
             * and the energy gain sums power in the segment with
 
808
             * potential and kinetic variations:
 
809
             *      Egain = sum(eta * power * dt) +
 
810
             *              totalMass * (g * (altInit - alt) +
 
811
             *              0.5 * (vInit*vInit - v*v))
 
812
             */
 
813
            foreach(const RideFilePoint *p1, ride->dataPoints()) {
 
814
                // Unpack:
 
815
                double power = max(0, p1->watts);
 
816
                double v     = p1->kph/vfactor;
 
817
                double distance = v * dt;
 
818
                double headwind = v;
 
819
                if( dataPresent->headwind ) {
 
820
                    headwind   = p1->headwind/vfactor;
 
821
                }
 
822
                double alt = p1->alt;
 
823
                // start initial segment
 
824
                if (nSeg < 0 && alt != 0) {
 
825
                    nSeg = 0;
 
826
                    X1[nSeg] = X2[nSeg] = Egain[nSeg] = 0.0;
 
827
                    altInit = alt;
 
828
                    vInit = v;
 
829
                }
 
830
                // accumulate segment data
 
831
                if (nSeg >= 0) {
 
832
                    // X1[nSgeg] * CdA == Aero-Loss
 
833
                    X1[nSeg] += 0.5 * rho * headwind*headwind * distance;
 
834
                    // X2[nSgeg] * Crr == RR-Loss
 
835
                    X2[nSeg] += totalMass * g * distance;
 
836
                    // Energy supplied
 
837
                    Egain[nSeg] += eta * power * dt;
 
838
                }
 
839
                // close current segment and start a new one
 
840
                if (nSeg >= 0 && alt != 0) {
 
841
                    // Add change in potential and kinetic energy
 
842
                    Egain[nSeg] += totalMass * (g * (altInit - alt) + 0.5 * (vInit*vInit - v*v));
 
843
                    // Start a new segment
 
844
                    nSeg++;
 
845
                    X1[nSeg] = X2[nSeg] = Egain[nSeg] = 0.0;
 
846
                    altInit = alt;
 
847
                    vInit = v;
 
848
                }
 
849
            }
 
850
            /* At least two segmentes needed to approximate:
 
851
             *     X1 * CdA + X2 * Crr = Egain
 
852
             * which, in matrix form, is:
 
853
             *    X * [ CdA ; Crr ] = Egain
 
854
             * and pre-multiplying by X transpose (X'):
 
855
             *    X'* X [ CdA ; Crr ] = X' * Egain
 
856
             * which is a system with two equations and two unknowns
 
857
             *    A * [ CdA ; Crr ] = B
 
858
             */
 
859
            if (nSeg >= 2) {
 
860
                // compute the normal equation: A * [ CdA ; Crr ] = B
 
861
                // A = X'*X
 
862
                // B = X'*Egain
 
863
                double A11 = 0, A12 = 0, A21 = 0, A22 = 0, B1 = 0, B2 = 0;
 
864
                for (int i = 0; i < nSeg; i++) {
 
865
                    A11 += X1[i] * X1[i];
 
866
                    A12 += X1[i] * X2[i];
 
867
                    A21 += X2[i] * X1[i];
 
868
                    A22 += X2[i] * X2[i];
 
869
                    B1  += X1[i] * Egain[i];
 
870
                    B2  += X2[i] * Egain[i];
 
871
                }
 
872
                // Solve the normal equation
 
873
                // A11 * CdA + A12 * Crr = B1
 
874
                // A21 * CdA + A22 * Crr = B2
 
875
                double det = A11 * A22 - A12 * A21;
 
876
                if (fabs(det) > 0.00) {
 
877
                    // round and update if the values are in Aerolab's range
 
878
                    double cda = floor(10000 * (A22 * B1 - A12 * B2) / det + 0.5) / 10000;
 
879
                    double crr = floor(1000000 * (A11 * B2 - A21 * B1) / det + 0.5) / 1000000;
 
880
                    if (cda >= 0.001 and cda <= 1.0 and crr >= 0.0001 and crr <= 0.1) {
 
881
                        this->cda = cda;
 
882
                        this->crr = crr;
 
883
                        errMsg = ""; // No error
 
884
                    } else {
 
885
                        errMsg = tr("Estimates out-of-range");
 
886
                    }
 
887
                } else {
 
888
                    errMsg = tr("At least two segments must be independent");
 
889
                }
 
890
            } else {
 
891
                errMsg = tr("At least two segments must be defined");
 
892
            }
 
893
        } else {
 
894
            errMsg = tr("Altitude and Power data must be present");
 
895
        }
 
896
    } else {
 
897
        errMsg = tr("No ride selected");
 
898
    }
 
899
    return errMsg;
 
900
}