diff --git a/src/CPPlot.cpp b/src/CPPlot.cpp
index 0eac9c951..08079e992 100644
--- a/src/CPPlot.cpp
+++ b/src/CPPlot.cpp
@@ -50,6 +50,9 @@
 
 CPPlot::CPPlot(QWidget *parent, Context *context, bool rangemode) : QwtPlot(parent), parent(parent),
 
+    // model
+    model(0), modelVariant(0),
+
     // state
     context(context), rideCache(NULL), bestsCache(NULL), rideSeries(RideFile::watts), isFiltered(false), shadeMode(2),
     shadeIntervals(true), rangemode(rangemode), showPercent(false), showHeat(false), showHeatByDate(false),
@@ -498,13 +501,7 @@ CPPlot::plotModel()
 
             case 4:
             {
-                cp  = tau = t0  = 0;
-                deriveCPParameters();
-
-                // ooopsie no model for us!
-                if (cp == 0 && tau == 0 && t0 == 0) return;
-
-                // the Veloclinic model has the following formulation
+                // the Michael Puchowicz (aka @Veloclinic) model has the following formulation
                 //
                 // p(t) = pc1 + pc2
                 //        Power at time t is the sum of;
@@ -518,37 +515,44 @@ CPPlot::plotModel()
                 //      p1   - pmax - cp as derived from the CP2-20 model
                 //      p2   - cp as derived from the CP2-20 model
                 //      tau1 - W'1 / p1
-                //      tau2 = 15,000
-                //      w2  -  A slow twitch W' derived from p2 * tau2
-                //      alpha - 0.1 
-                //      beta
+                //      tau2 - 15,000
+                //      w2   -  A slow twitch W' derived from p2 * tau2
+                //      alpha- 0.1 
+                //      beta - 1.0
                 //
                 // Fast twitch component is:
                 //      pc1(t) = W'1 / t * (1-exp(-t/tau1)) * ((1-exp(-t/10)) ^ (1/alpha))
                 // 
                 // Slow twitch component has three formulations:
-                //      a) pc2(t) = p2 * tau2 * (1-exp(-t/tau2))
-                //      b) pc2(t) = p2 / (1 + t/tau2)
-                //      c) pc2(t) = p2 / (1 + t/5400) ^ (1/beta)
+                //      sprint capped linear)          pc2(t) = p2 * tau2 * (1-exp(-t/tau2))
+                //      sprint capped regeneration)    pc2(t) = p2 / (1 + t/tau2)
+                //      sprint capped exponential)     pc2(t) = p2 / (1 + t/5400) ^ (1/beta)
                 //
                 // Currently deciding which of the three formulations to use
                 // as the base for GoldenCheetah (we have enough models already !)
                 //
                 // to keep things simple we just use formulation (a) for now.
+                cp  = tau = t0  = 0;
+                deriveCPParameters();
+
+                // ooopsie no model for us!
+                if (cp == 0 && tau == 0 && t0 == 0) return;
+
                 double pmax = cp * (double(1.00f)+tau /(((double(1)/double(60))+t0)));
                 double w1 = cp*tau*60;
                 double p1 = pmax - cp;
                 double p2 = cp;
                 double tau1 = w1 / p1;
-                double tau2 = 15000;
-                double alpha = 0.1;
+                const double tau2 = 15000;
+                const double alpha = 0.1;
+                const double beta = 1.0;
 
                 //double w2 = p2 * tau2;
                 //qDebug()<<"model parameters: pmax="<<pmax<<"w1="<<w1<<"p1="<<p1
                 //            <<"p2="<<p2<<"tau1="<<tau1<<"tau2="<<tau2<<"alpha="<<alpha;
 
                 // populate curve data with a CP curve of 100 points resolution
-                const int points = 3600 * 5; // 5 hours is enough
+                const int points = 3600 * 10; // 10 hours is enough
                 QVector<double> cp_curve_power(points);
                 QVector<double> cp_curve_time(points);
 
@@ -562,8 +566,27 @@ CPPlot::plotModel()
                     } else {
 
                         // two component model
-                        double pc1 = w1 / t * (1.00f - exp(-t/tau1)) * pow(1-exp(-t/10), alpha); // ^ (1/alpha) missing XXX
-                        double pc2 = p2 * tau2 / t * (1-exp(-t/tau2));
+                        double pc1 = w1 / t * (1.00f - exp(-t/tau1)) * pow(1-exp(-t/10), alpha); 
+
+                        // which variant for pc2 ?
+                        double pc2 = 0.0f;
+                        switch (modelVariant) {
+
+                            default:
+                            case 0 : // exponential top and bottom
+                                pc2 = p2 * tau2 / t * (1-exp(-t/tau2));
+                                break;
+
+                            case 1 : // linear feedback
+                                pc2 = p2 / (1+t/tau2);
+                                break;
+
+                            case 2 : // regeneration
+                                pc2 = pow(p2 / (1+t/5400),1/beta);
+                                //pc2 = p2 / pow((1+t/5400),(1/beta));
+                                break;
+
+                        }
 
                         // p(t) as a sum of both component powers
                         cp_curve_power[i] = pc1 + pc2;
@@ -1293,7 +1316,7 @@ CPPlot::setShadeIntervals(int x)
 
 // model parameters!
 void
-CPPlot::setModel(int sanI1, int sanI2, int anI1, int anI2, int aeI1, int aeI2, int laeI1, int laeI2, int model)
+CPPlot::setModel(int sanI1, int sanI2, int anI1, int anI2, int aeI1, int aeI2, int laeI1, int laeI2, int model, int variant)
 {
     this->anI1 = double(anI1) / double(60.00f);
     this->anI2 = double(anI2) / double(60.00f);
@@ -1306,6 +1329,7 @@ CPPlot::setModel(int sanI1, int sanI2, int anI1, int anI2, int aeI1, int aeI2, i
     this->laeI2 = double(laeI2) / double(60.00f);
 
     this->model = model;
+    this->modelVariant = variant;
     clearCurves();
 }
 
diff --git a/src/CPPlot.h b/src/CPPlot.h
index 87f772225..6bf2b2e9a 100644
--- a/src/CPPlot.h
+++ b/src/CPPlot.h
@@ -70,7 +70,7 @@ class CPPlot : public QwtPlot
         void setSeries(CriticalPowerWindow::CriticalSeriesType);
         void setPlotType(int index);
         void setModel(int sanI1, int sanI2, int anI1, int anI2, 
-                      int aeI1, int aeI2, int laeI1, int laeI2, int model);
+                      int aeI1, int aeI2, int laeI1, int laeI2, int model, int variant);
 
         // getters
         QVector<double> getBests();
@@ -119,7 +119,7 @@ class CPPlot : public QwtPlot
         void refreshReferenceLines(RideItem*);
 
         // Models and Extended Models
-        int model;
+        int model, modelVariant;
         double sanI1, sanI2, anI1, anI2, aeI1, aeI2, laeI1, laeI2;
         double cp, tau, t0; // CP model parameters
         ExtendedCriticalPower *ecp;
diff --git a/src/CriticalPowerWindow.cpp b/src/CriticalPowerWindow.cpp
index 2ce0d54ae..66becda4a 100644
--- a/src/CriticalPowerWindow.cpp
+++ b/src/CriticalPowerWindow.cpp
@@ -196,7 +196,6 @@ CriticalPowerWindow::CriticalPowerWindow(const QDir &home, Context *context, boo
     mcl->addRow(new QLabel(tr("CP Model")), modelCombo);
 
     mcl->addRow(new QLabel(tr(" ")));
-
     intervalLabel = new QLabel(tr("Search Interval"));
     secondsLabel = new QLabel(tr("(seconds)"));
     mcl->addRow(intervalLabel, secondsLabel);
@@ -293,6 +292,15 @@ CriticalPowerWindow::CriticalPowerWindow(const QDir &home, Context *context, boo
     laeLayout->addWidget(laeI2SpinBox);
     mcl->addRow(laeLabel, laeLayout);
 
+    mcl->addRow(new QLabel(""), new QLabel(""));
+    velo1 = new QRadioButton(tr("Exponential"));
+    velo1->setChecked(true);
+    mcl->addRow(new QLabel(tr("Variant")), velo1);
+    velo2 = new QRadioButton(tr("Linear feedback"));
+    mcl->addRow(new QLabel(""), velo2);
+    velo3 = new QRadioButton(tr("Regeneration"));
+    mcl->addRow(new QLabel(""), velo3);
+
     // point 2 + 3 -or- point 1 + 2 in a 2 point model
 
     grid = new QwtPlotGrid();
@@ -413,6 +421,9 @@ CriticalPowerWindow::CriticalPowerWindow(const QDir &home, Context *context, boo
     connect(sanI2SpinBox, SIGNAL(valueChanged(double)), this, SLOT(modelParametersChanged()));
     connect(laeI1SpinBox, SIGNAL(valueChanged(double)), this, SLOT(modelParametersChanged()));
     connect(laeI2SpinBox, SIGNAL(valueChanged(double)), this, SLOT(modelParametersChanged()));
+    connect(velo1, SIGNAL(toggled(bool)), this, SLOT(modelParametersChanged()));
+    connect(velo2, SIGNAL(toggled(bool)), this, SLOT(modelParametersChanged()));
+    connect(velo3, SIGNAL(toggled(bool)), this, SLOT(modelParametersChanged()));
 
     // redraw on config change -- this seems the simplest approach
     connect(context, SIGNAL(filterChanged()), this, SLOT(forceReplot()));
@@ -507,6 +518,13 @@ CriticalPowerWindow::modelChanged()
     // for best results in predicting both W' and CP and providing
     // a reasonable fit for durations < 2mins.
     active = true;
+
+    // hide veloclinic's variation for everyone
+    // it will get shown for model 4 below
+    velo1->hide();
+    velo2->hide();
+    velo3->hide();
+
     switch (modelCombo->currentIndex()) {
 
     case 0 : // None
@@ -529,7 +547,14 @@ CriticalPowerWindow::modelChanged()
             // No default values !
             break;
 
-    case 4 : // Veloclinic Model uses 2 parameter classic
+    case 4 : // Veloclinic Model uses 2 parameter classic but 
+             // also lets you select a variation ..
+            velo1->show();
+            velo2->show();
+            velo3->show();
+
+            // and drop through into case 1 below ...
+
     case 1 : // Classic 2 param model 2-20 default (per literature)
 
             intervalLabel->show();
@@ -553,6 +578,7 @@ CriticalPowerWindow::modelChanged()
             anI2SpinBox->setValue(120);
             aeI1SpinBox->setValue(1000);
             aeI2SpinBox->setValue(1200);
+
             break;
 
     case 2 : // 3 param model: 3-30 model
@@ -615,6 +641,27 @@ CriticalPowerWindow::modelChanged()
     modelParametersChanged();
 }
 
+// kind of tedious but return index into a radio button group 
+// for the button that is actually checked.
+int 
+CriticalPowerWindow::variant() const
+{
+    if (velo1->isChecked()) return 0;
+    if (velo2->isChecked()) return 1;
+    if (velo3->isChecked()) return 2;
+    return 0; // default
+}
+
+void
+CriticalPowerWindow::setVariant(int index)
+{
+    switch (index) {
+    case 0 : velo1->setChecked(true); velo2->setChecked(false); velo3->setChecked(false); break;
+    case 1 : velo1->setChecked(false); velo2->setChecked(true); velo3->setChecked(false); break;
+    case 2 : velo1->setChecked(false); velo2->setChecked(false); velo3->setChecked(true); break;
+    }
+}
+
 void
 CriticalPowerWindow::modelParametersChanged()
 {
@@ -636,7 +683,8 @@ CriticalPowerWindow::modelParametersChanged()
                         aeI2SpinBox->value(),
                         laeI1SpinBox->value(),
                         laeI2SpinBox->value(),
-                        modelCombo->currentIndex());
+                        modelCombo->currentIndex(),
+                        variant());
 
     // and apply
     if (amVisible() && myRideItem != NULL) {
diff --git a/src/CriticalPowerWindow.h b/src/CriticalPowerWindow.h
index 6e0214029..280cf40cc 100644
--- a/src/CriticalPowerWindow.h
+++ b/src/CriticalPowerWindow.h
@@ -58,6 +58,7 @@ class CriticalPowerWindow : public GcChartWindow
     Q_PROPERTY(QString season READ season WRITE setSeason USER true)
 
     Q_PROPERTY(int cpmodel READ cpModel WRITE setCPModel USER true)
+    Q_PROPERTY(int variant READ variant WRITE setVariant USER true)
     Q_PROPERTY(int ani1 READ anI1 WRITE setAnI1 USER true)
     Q_PROPERTY(int ani2 READ anI2 WRITE setAnI2 USER true)
     Q_PROPERTY(int aei1 READ aeI1 WRITE setAeI1 USER true)
@@ -102,6 +103,9 @@ class CriticalPowerWindow : public GcChartWindow
         int cpModel() const { return modelCombo->currentIndex(); }
         void setCPModel(int x) { modelCombo->setCurrentIndex(x); }
 
+        int variant() const;
+        void setVariant(int x);
+
 #ifdef GC_HAVE_LUCENE
         // filter
         bool isFiltered() const { return (searchBox->isFiltered() || context->ishomefiltered || context->isfiltered); }
@@ -252,6 +256,7 @@ class CriticalPowerWindow : public GcChartWindow
         QLabel *cpintCPValue;
         QComboBox *seriesCombo;
         QComboBox *modelCombo;
+        QRadioButton *velo1, *velo2, *velo3; // for selecting veloclinic formulation
         QComboBox *cComboSeason;
         QComboBox *ridePlotStyleCombo;
         QComboBox *shadeCombo;