Changeset 3148 for trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc

Timestamp:

02/14/04 11:48:43 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc (modified) (7 diffs)

trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc

@@ -24 +24 @@
 
 //////////////////////////////////////////////////////////////////////////////
-//                                                                          // 
-//  Cubic Spline Interpolation                                              //
-//                                                                          //
+//
+//  Cubic Spline Interpolation
+//
 //////////////////////////////////////////////////////////////////////////////
-
 #include "MCubicCoeff.h"
 
-#include "TMath.h"
+#include <TMath.h>
 
 #include "MLog.h"
@@ -45 +44 @@
 // where x is the independent variable, 2 and 3 are exponents
 //
-
 MCubicCoeff::MCubicCoeff(Double_t x, Double_t xNext, Double_t y, Double_t yNext, 
                          Double_t a, Double_t b, Double_t c) : 
@@ -51 +49 @@
 {
     fH = fXNext - fX;
-    if(!EvalMinMax())
-    {
-        gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl;
-        fMin = 0;
-        fMax = 0;
-    }
+    if (EvalMinMax())
+        return;
+
+    gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl;
+    fMin = 0;
+    fMax = 0;
 }
 
@@ -66 +64 @@
 Double_t MCubicCoeff::Eval(Double_t x)
 {
-    Double_t dx = x - fX;
-    return (fY+dx*(fC+dx*(fB+dx*fA)));
+    const Double_t dx = x - fX;
+    return fY + dx*(fC + dx*(fB + dx*fA));
 }
 
@@ -79 +77 @@
 Bool_t MCubicCoeff::EvalMinMax()
 {
-    fMin = fMax = fY;
-    fAbMin = fAbMax = fX;
+    fMin = fY;
+    fMax = fY;
+
+    fAbMin = fX;
+    fAbMax = fX;
+
     if (fYNext < fMin)
     {
-        fMin = fYNext;
+        fMin   = fYNext;
         fAbMin = fXNext;
     }
     if (fYNext > fMax)
     {
-        fMax = fYNext;
+        fMax   = fYNext;
         fAbMax = fXNext;
     }
-    Double_t tempMinMax;
-    Double_t delta = 4.0*fB*fB - 12.0*fA*fC;
-        if (delta >= 0.0 && fA != 0.0)
-        {
-            Double_t sqrtDelta = sqrt(delta);
-            Double_t xPlus  = (-2.0*fB + sqrtDelta)/(6.0*fA);
-            Double_t xMinus = (-2.0*fB - sqrtDelta)/(6.0*fA);
-            if (xPlus >= 0.0 && xPlus <= fH)
-            {
-                tempMinMax = this->Eval(fX+xPlus);
-                if (tempMinMax < fMin)
-                {
-                    fMin = tempMinMax;
-                    fAbMin = fX + xPlus;
-                }
-                if (tempMinMax > fMax)
-                {
-                    fMax = tempMinMax;
-                    fAbMax = fX + xPlus;
-                }
-            }
-            if (xMinus >= 0.0 && xMinus <= fH)
-            {
-                tempMinMax = this->Eval(fX+xMinus);
-                if (tempMinMax < fMin)
-                {
-                    fMin = tempMinMax;
-                    fAbMin = fX + xMinus;
-                }
-                if (tempMinMax > fMax)
-                {
-                    fMax = tempMinMax;
-                    fAbMax = fX + xMinus;
-                }
-            }
-        }
-/* if fA is zero then we have only one possible solution */
-        else if (fA == 0.0 && fB != 0.0)
-        {
-            Double_t xSolo = - (fC/(2.0*fB));
-            if (xSolo >= 0.0 && xSolo <= fH)
-            {
-                tempMinMax = this->Eval(fX+xSolo);
-                if (tempMinMax < fMin)
-                {
-                    fMin = tempMinMax;
-                    fAbMin = fX + xSolo;
-                }
-                if (tempMinMax > fMax)
-                {
-                    fMax = tempMinMax;
-                    fAbMax = fX + xSolo;
-                }
-            }
-        }
+
+    const Double_t delta = fB*fB*4 - fA*fC*12;
+    if (delta >= 0 && fA != 0)
+    {
+        const Double_t sqrtDelta = TMath::Sqrt(delta);
+
+        const Double_t xPlus  = (-fB*2 + sqrtDelta)/(fA*6);
+        const Double_t xMinus = (-fB*2 - sqrtDelta)/(fA*6);
+
+        if (xPlus >= 0 && xPlus <= fH)
+        {
+            const Double_t tempMinMax = Eval(fX+xPlus);
+            if (tempMinMax < fMin)
+            {
+                fMin = tempMinMax;
+                fAbMin = fX + xPlus;
+            }
+            if (tempMinMax > fMax)
+            {
+                fMax = tempMinMax;
+                fAbMax = fX + xPlus;
+            }
+        }
+        if (xMinus >= 0 && xMinus <= fH)
+        {
+            const Double_t tempMinMax = Eval(fX+xMinus);
+            if (tempMinMax < fMin)
+            {
+                fMin = tempMinMax;
+                fAbMin = fX + xMinus;
+            }
+            if (tempMinMax > fMax)
+            {
+                fMax = tempMinMax;
+                fAbMax = fX + xMinus;
+            }
+        }
+        return kTRUE;
+    }
+
+    /* if fA is zero then we have only one possible solution */
+    if (fA == 0 && fB != 0)
+    {
+        const Double_t xSolo = -fC/(fB*2);
+
+        if (xSolo < 0 || xSolo > fH)
+            return kTRUE;
+
+        const Double_t tempMinMax = Eval(fX+xSolo);
+        if (tempMinMax < fMin)
+        {
+            fMin = tempMinMax;
+            fAbMin = fX + xSolo;
+        }
+        if (tempMinMax > fMax)
+        {
+            fMax = tempMinMax;
+            fAbMax = fX + xSolo;
+        }
+        return kTRUE;
+    }
+
     return kTRUE;
 }
@@ -157 +166 @@
 // There could be three or one real solutions
 //
-
 Short_t MCubicCoeff::FindCardanRoot(Double_t y, Double_t *x)
 {
-
-    Short_t whichRoot = -1;
-
-    Double_t a = fB/fA;
-    Double_t b = fC/fA;
-    Double_t c = (fY - y)/fA;
-
-    Double_t q = (a*a-3.0*b)/9.0;
-    Double_t r = (2.0*a*a*a-9.0*a*b+27.0*c)/54.0;
-
-    Double_t aOver3 = a/3.0;
-    Double_t r2 = r*r;
-    Double_t q3 = q*q*q;
+    const Double_t a = fB/fA;
+    const Double_t b = fC/fA;
+    const Double_t c = (fY - y)/fA;
+
+    const Double_t q = (a*a - b*3)/9;
+    const Double_t r = (a*a*a*2 - a*b*9 + c*27)/54;
+
+    const Double_t aOver3 = a/3;
+    const Double_t r2 = r*r;
+    const Double_t q3 = q*q*q;
     
     if (r2 < q3) //3 real sol
     {
-        Double_t sqrtQ = TMath::Sqrt(q);
-        Double_t min2SqQ = -2.0*sqrtQ;
-        Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ));
-
-        x[0] = min2SqQ * TMath::Cos(theta/3.0) - aOver3;
-        x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3.0) - aOver3;
-        x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3.0) - aOver3;
+        const Double_t sqrtQ = TMath::Sqrt(q);
+        const Double_t min2SqQ = -sqrtQ*2;
+        const Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ));
+
+        x[0] = min2SqQ * TMath::Cos(theta/3) - aOver3;
+        x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3) - aOver3;
+        x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3) - aOver3;
+
         for (Int_t i = 0; i < 3; i++)
-            if (x[i] >= 0.0 && x[i] <= fH)
+            if (x[i] >= 0 && x[i] <= fH)
             {
-                x[i] = x[i] + fX;
-                whichRoot = i;
-                break;
+                x[i] += fX;
+                return i;
             }
-    }
-    else //only 1 real sol
-    {
-        Double_t s;
-        if (r == 0.0)
-            s = 0.0;
-        else if (r < 0.0)
-            s = TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
-        else // r > 0.0
-            s = - TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
-        if (s == 0.0)
-            x[0] = - aOver3;
-        else
-            x[0] = (s + q/s) - aOver3;
-        if (x[0] >= 0.0 && x[0] <= fH)
-        {
-            x[0] = x[0] + fX;
-            whichRoot = 0;
-        }
-    }
-    return whichRoot;
+        return -1;
+    }
+
+    const Double_t s = r==0 ? 0 : -TMath::Sign(TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3), 1./3), r);
+
+    x[0] = s==0 ? - aOver3 : (s + q/s) - aOver3;
+
+    if (x[0] < 0 || x[0] > fH)
+        return -1;
+
+    x[0] += fX;
+    return 0;
 }
 
@@ -220 +216 @@
 Bool_t MCubicCoeff :: IsIn(Double_t x)
 {
-    if (x >= fX && x <= fXNext)
-        return kTRUE;
-    else
-        return kFALSE;
-}
+    return x >= fX && x <= fXNext;
+}