Merge branch '157-cyclotron-trim-coil-has-discontinuous-derivative' into 'master'

Resolve "Cyclotron trim coil has discontinuous derivative"

See merge request !35

src/BasicActions/Option.cpp

@@ -61,10 +61,10 @@ namespace {
         REBINFREQ,
         SCSOLVEFREQ,
         MTSSUBSTEPS,
-	REMOTEPARTDEL,
+        REMOTEPARTDEL,
         RHODUMP,
         EBDUMP,
-	CSRDUMP,
+        CSRDUMP,
         AUTOPHASE,
         PPDEBUG,
         SURFDUMPFREQ,
@@ -76,7 +76,7 @@ namespace {
         ENABLEHDF5,
         ASCIIDUMP,
         BOUNDPDESTROYFQ,
-	BEAMHALOBOUNDARY,
+        BEAMHALOBOUNDARY,
         CLOTUNEONLY,
         IDEALIZED,
         LOGBENDTRAJECTORY,

src/Classic/AbsBeamline/Cyclotron.cpp

@@ -19,12 +19,14 @@
 // ------------------------------------------------------------------------
 
 #include "AbsBeamline/Cyclotron.h"
-#include "Algorithms/PartBunchBase.h"
+
 #include "AbsBeamline/BeamlineVisitor.h"
+#include "Algorithms/PartBunchBase.h"
+#include "Fields/Fieldmap.h"
 #include "Physics/Physics.h"
 #include "Structure/LossDataSink.h"
+#include "TrimCoils/TrimCoil.h"
 #include "Utilities/Options.h"
-#include "Fields/Fieldmap.h"
 #include "Utilities/GeneralClassicException.h"
 #include <fstream>
 
@@ -62,10 +64,7 @@ Cyclotron::Cyclotron(const Cyclotron &right):
     type_m(right.type_m),
     harm_m(right.harm_m),
     bscale_m(right.bscale_m),
-    tcr1V_m(right.tcr1V_m),
-    tcr2V_m(right.tcr2V_m),
-    mbtcV_m(right.mbtcV_m),
-    slptcV_m(right.slptcV_m),
+    trimcoils_m(right.trimcoils_m),
     minr_m(right.minr_m),
     maxr_m(right.maxr_m),
     minz_m(right.minz_m),
@@ -88,69 +87,9 @@ Cyclotron::~Cyclotron() {
 
 
 void Cyclotron::applyTrimCoil(const double r, const double z, double *br, double *bz) {
-    /// update bz and br with trim coil contributions
-    // for some discussion on the formulas see
-    // http://doi.org/10.1103/PhysRevSTAB.14.054402
-    // https://gitlab.psi.ch/OPAL/src/issues/157
-    // https://gitlab.psi.ch/OPAL/src/issues/110
-
-    // unitless constants
-    const double Amax1 = 1;
-    const double Amax2 = 3;
-    const double Amin  = -2;
-    const double x01   = 4;
-    const double x02   = 8;
-    const double h1    = 0.03;
-    const double h2    = 0.2;
-    const double justAnotherFudgeFactor = 1 / 2.78;
-    // helper variables
-    const double log10  = std::log(10);
-    const double const3 = -(Amax1 - Amin) * h1 * log10;
-    const double const4 =  (Amax2 - Amin) * h2 * log10;
-
-    for (unsigned int idx = 0; idx < slptcV_m.size(); ++ idx) {
-        const double &tcr1 = tcr1V_m[idx];
-        const double &tcr2 = tcr2V_m[idx];
-        const double &slope = slptcV_m[idx];
-        const double &magnitude = mbtcV_m[idx];
-
-        double part1;
-        double part2;
-
-        if (2 * r < (tcr2 + tcr1)) {
-            part1 = std::pow(10.0,  ((r - tcr1) * slope - x01) * h1);
-            part2 = std::pow(10.0, -((r - tcr1) * slope - x02) * h2);
-        } else {
-            part1 = std::pow(10.0,  ((tcr2 - r) * slope - x01) * h1);
-            part2 = std::pow(10.0, -((tcr2 - r) * slope - x02) * h2);
-        }
-
-        const double part1plusinv = 1 / (1 + part1);
-        const double part2plusinv = 1 / (1 + part2);
-
-        double part3 = const3 * slope * part1 * part1plusinv * part1plusinv;
-        double part4 = const4 * slope * part2 * part2plusinv * part2plusinv;
-
-        const double constmag = magnitude * justAnotherFudgeFactor;
-
-        double dr  = constmag * (part3 + part4);
-        double btr = constmag * (Amin - 1 +
-                                 (Amax1 - Amin) * part1plusinv +
-                                 (Amax2 - Amin) * part2plusinv);
-
-        if (std::isnan(dr)  || std::isinf(dr) ||
-            std::isnan(btr) || std::isinf(btr)) {
-            ERRORMSG("r = " << r << " m, tcr1 = " << tcr1 << " m, tcr2 = " << tcr2 << " m\n");
-            ERRORMSG("slope = " << slope << ", magnitude = " << magnitude << " kG\n");
-            ERRORMSG("part1 = " << part1 << ", part2 = " << part2 << "\n");
-            ERRORMSG("part3 = " << part3 << ", part4 = " << part4 << endl);
-            throw GeneralClassicException("Cyclotron::applyTrimCoil",
-                                          "dr or btr yield results that are either nan or inf");
-        }
-
-        *bz -= btr;
-        *br -= dr * z;
-    }
+     for (auto trimcoil : trimcoils_m) {
+         trimcoil->applyField(r,z,br,bz);
+     }
 }
 
 void Cyclotron::accept(BeamlineVisitor &visitor) const {
@@ -283,7 +222,7 @@ void Cyclotron::setEScale(vector<double> s) {
 }
 
 unsigned int Cyclotron::getNumberOfTrimcoils() const {
-  return tcr1V_m.size();
+  return trimcoils_m.size();
 }
 
 double Cyclotron::getCyclHarm() const {
@@ -336,20 +275,8 @@ double Cyclotron::getMaxZ() const {
     return maxz_m;
 }
 
-void Cyclotron::setTCr1V(const vector<double> & tcr1) {
-    tcr1V_m = tcr1;
-}
-
-void Cyclotron::setTCr2V(const vector<double> & tcr2) {
-    tcr2V_m = tcr2;
-}
-
-void Cyclotron::setMBtcV(const vector<double> & mbtc) {
-    mbtcV_m = mbtc;
-}
-
-void Cyclotron::setSLPtcV(const vector<double> & slptc) {
-    slptcV_m = slptc;
+void Cyclotron::setTrimCoils(const std::vector<TrimCoil*> &trimcoils) {
+    trimcoils_m = trimcoils;
 }
 
 void Cyclotron::setFMLowE(double e) { fmLowE_m = e;}

src/Classic/AbsBeamline/Cyclotron.h

@@ -29,6 +29,7 @@
 
 class Fieldmap;
 class LossDataSink;
+class TrimCoil;
 
 enum BFieldType {PSIBF,CARBONBF,ANSYSBF,AVFEQBF,FFAGBF,BANDRF,SYNCHRO};
 
@@ -158,10 +159,7 @@ public:
 
     void   setEScale(std::vector<double> bs);
 
-    void setTCr1V(const std::vector<double> &tcr1);
-    void setTCr2V(const std::vector<double> &tcr2);
-    void setMBtcV(const std::vector<double> &mbtc);
-    void setSLPtcV(const std::vector<double> &slptc);
+    void setTrimCoils(const std::vector<TrimCoil*> &trimcoils);
 
     void setSuperpose(std::vector<bool> flag);
     //    virtual bool getSuperpose() const;
@@ -243,12 +241,8 @@ private:
 
     double bscale_m; // a scale factor for the B-field
 
-    ///@{ Trim coil variable
-    std::vector<double> tcr1V_m;
-    std::vector<double> tcr2V_m;
-    std::vector<double> mbtcV_m;
-    std::vector<double> slptcV_m;
-    ///@}
+    /// Trim coils
+    std::vector<TrimCoil*> trimcoils_m;
 
     double minr_m;
     double maxr_m;

src/Classic/CMakeLists.txt

@@ -18,6 +18,7 @@ add_subdirectory(Physics)
 add_subdirectory(Solvers)
 add_subdirectory(Structure)
 add_subdirectory(Utilities)
+add_subdirectory(TrimCoils)
 
 install (FILES Utilities/ClassicException.h Utilities/ClassicRandom.h Utilities/Options.h Utilities/OptionTypes.h 
     DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Utilities"

src/Classic/TrimCoils/CMakeLists.txt

+set (_SRCS
+    OpalTrimCoil.cpp
+    TrimCoil.cpp
+    TrimCoilFit.cpp
+    TrimCoilMirrored.cpp
+)
+
+include_directories (
+  ${CMAKE_CURRENT_SOURCE_DIR}
+  )
+
+ADD_OPAL_SOURCES(${_SRCS})
+
+set (HDRS
+    OpalTrimCoil.h
+    TrimCoil.h
+    TrimCoilFit.h
+    TrimCoilMirrored.h
+)
+
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/TrimCoils")
+install (FILES ${_SRCS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/TrimCoils")

src/Classic/TrimCoils/OpalTrimCoil.cpp

+#include "TrimCoils/OpalTrimCoil.h"
+
+#include "AbstractObjects/OpalData.h"
+#include "Attributes/Attributes.h"
+#include "TrimCoils/TrimCoilFit.h"
+#include "TrimCoils/TrimCoilMirrored.h"
+#include "Utilities/OpalException.h"
+#include "Utilities/Util.h"
+#include "Utility/IpplInfo.h"
+
+extern Inform *gmsg;
+
+
+// Class OpalTrimCoil
+// ------------------------------------------------------------------------
+
+// The attributes of class OpalTrimCoil.
+namespace {
+    enum {
+        TYPE,       // The type of trim coil
+        COEFNUM,    //
+        COEFDENOM,  //
+        BMAX,       //
+        RMIN,       //
+        RMAX,       //
+        SLPTC,
+        SIZE
+    };
+}
+
+OpalTrimCoil::OpalTrimCoil():
+    Definition(SIZE, "TRIMCOIL",
+               "The \"TRIMCOIL\" statement defines a trim coil."),
+    trimcoil_m(nullptr) {
+    itsAttr[TYPE]      = Attributes::makeString
+                         ("TYPE", "Specifies the type of trim coil: PSI-RING, PSI-RING-OLD");
+
+    itsAttr[COEFNUM]   = Attributes::makeRealArray
+                         ("COEFNUM", "List of polynomial coefficients for the numerator");
+
+    itsAttr[COEFDENOM] = Attributes::makeRealArray
+                         ("COEFDENOM", "List of polynomial coefficients for the denominator");
+
+    itsAttr[BMAX]      = Attributes::makeReal
+                         ("BMAX", "Maximum magnetic field in Tesla.");
+
+    itsAttr[RMIN]      = Attributes::makeReal
+                         ("RMIN", "Minimum radius in millimeters.");
+
+    itsAttr[RMAX]      = Attributes::makeReal
+                         ("RMAX", "Maximum radius in millimeters.");
+
+    itsAttr[SLPTC]      = Attributes::makeReal
+                         ("SLPTC", "Slopes of the rising edge [1/mm] (for PSI-RING-OLD)");
+    
+
+    registerOwnership(AttributeHandler::STATEMENT);
+
+    OpalTrimCoil *defTrimCoil = clone("UNNAMED_TRIMCOIL");
+    defTrimCoil->builtin = true;
+
+    try {
+        defTrimCoil->update();
+        OpalData::getInstance()->define(defTrimCoil);
+    } catch(...) {
+        delete defTrimCoil;
+    }
+}
+
+
+OpalTrimCoil::OpalTrimCoil(const std::string &name, OpalTrimCoil *parent):
+    Definition(name, parent),
+    trimcoil_m(nullptr)
+{}
+
+
+OpalTrimCoil::~OpalTrimCoil() {
+    
+}
+
+
+bool OpalTrimCoil::canReplaceBy(Object *object) {
+    // Can replace only by another trim coil.
+    return dynamic_cast<OpalTrimCoil *>(object) != nullptr;
+}
+
+
+OpalTrimCoil *OpalTrimCoil::clone(const std::string &name) {
+    return new OpalTrimCoil(name, this);
+}
+
+
+void OpalTrimCoil::execute() {
+    update();
+}
+
+
+OpalTrimCoil *OpalTrimCoil::find(const std::string &name) {
+    OpalTrimCoil *trimcoil = dynamic_cast<OpalTrimCoil *>(OpalData::getInstance()->find(name));
+
+    if (trimcoil == nullptr) {
+        throw OpalException("OpalTrimCoil::find()", "OpalTrimCoil \"" + name + "\" not found.");
+    }
+    return trimcoil;
+}
+
+
+void OpalTrimCoil::update() {
+    // Set default name.
+    if (getOpalName().empty()) setOpalName("UNNAMED_TRIMCOIL");
+}
+
+
+void OpalTrimCoil::initOpalTrimCoil() {
+    if (trimcoil_m == nullptr) {
+        
+        std::string type = Util::toUpper(Attributes::getString(itsAttr[TYPE]));
+        
+        double bmax = Attributes::getReal(itsAttr[BMAX]);
+        double rmin = Attributes::getReal(itsAttr[RMIN]);
+        double rmax = Attributes::getReal(itsAttr[RMAX]);
+        
+        if (type == "PSI-RING") {
+            std::vector<double> coefnum   = Attributes::getRealArray(itsAttr[COEFNUM]);
+            std::vector<double> coefdenom = Attributes::getRealArray(itsAttr[COEFDENOM]);
+            trimcoil_m = std::unique_ptr<TrimCoilFit>      (new TrimCoilFit(bmax, rmin, rmax, coefnum, coefdenom));
+        } else if (type == "PSI-RING-OLD") {
+            double slope = Attributes::getReal(itsAttr[SLPTC]);
+            trimcoil_m = std::unique_ptr<TrimCoilMirrored> (new TrimCoilMirrored(bmax, rmin, rmax, slope));
+        } else {
+            WARNMSG(type << " is not a valid trim coil type" << endl);
+        }
+        
+        *gmsg << level3 << *this << endl;
+    }
+}
+
+Inform& OpalTrimCoil::print(Inform &os) const {
+    os << "* ******************************** T R I M C O I L ********************************\n"
+       << "* TRIMCOIL       " << getOpalName() << '\n'
+       << "* TYPE           " << Attributes::getString(itsAttr[TYPE]) << '\n';
+       
+    if (Util::toUpper(Attributes::getString(itsAttr[TYPE])) == "PSI-RING") {
+        std::vector<double> coefnum = Attributes::getRealArray(itsAttr[COEFNUM]);
+        std::stringstream ss;
+        for (std::size_t i = 0; i < coefnum.size(); ++i) {
+            ss << ((i > 0) ? "+ " : "") << coefnum[i]
+               << ((i > 0) ? (" * x^" + std::to_string(i)) : "") << ' ';
+        }
+        os << "* POLYNOM NUM    " << ss.str() << '\n';
+    
+        std::vector<double> coefdenom = Attributes::getRealArray(itsAttr[COEFDENOM]);
+        ss.str("");
+        for (std::size_t i = 0; i < coefdenom.size(); ++i) {
+            ss << ((i > 0) ? "+ " : "") << coefdenom[i]
+               << ((i > 0) ? (" * x^" + std::to_string(i)) : "") << ' ';
+        }
+        os << "* POLYNOM DENOM  " << ss.str() << '\n';
+    }
+
+    os << "* BMAX           " << Attributes::getReal(itsAttr[BMAX]) << '\n'
+       << "* RMIN           " << Attributes::getReal(itsAttr[RMIN]) << '\n'
+       << "* RMAX           " << Attributes::getReal(itsAttr[RMAX]) << '\n';
+ 
+    if (Util::toUpper(Attributes::getString(itsAttr[TYPE])) == "PSI-RING-OLD") {
+        os << "* SLPTC          " << Attributes::getReal(itsAttr[SLPTC]) << '\n';
+    }
+    os << "* *********************************************************************************" << endl;
+    return os;
+}

src/Classic/TrimCoils/OpalTrimCoil.h

+#ifndef OPAL_TRIM_COIL_H
+#define OPAL_TRIM_COIL_H
+
+#include <string>
+#include "AbstractObjects/Definition.h"
+
+class TrimCoil;
+
+// Class OpalTrimCoil
+// ------------------------------------------------------------------------
+/// The TRIMCOIL definition.
+//  A TRIMCOIL definition is used to define a trim coil which can be applied
+//  to a Cyclotron
+
+class OpalTrimCoil: public Definition {
+
+public:
+
+    /// Exemplar constructor.
+    OpalTrimCoil();
+
+    virtual ~OpalTrimCoil();
+
+    /// Test if replacement is allowed.
+    //  Can replace only by another OpalTrimCoil
+    virtual bool canReplaceBy(Object *object);
+
+    /// Make clone.
+    virtual OpalTrimCoil *clone(const std::string &name);
+
+    /// Check the OpalTrimCoil data.
+    virtual void execute();
+
+    /// Find named trim coil.
+    static OpalTrimCoil *find(const std::string &name);
+
+    /// Update the OpalTrimCoil data.
+    virtual void update();
+
+    /// Print method, called at initialisation
+    Inform& print(Inform& os) const;
+
+    /// Initialise implementation
+    void initOpalTrimCoil();
+
+    /// Actual implementation
+    std::unique_ptr<TrimCoil> trimcoil_m;
+    
+private:
+
+    ///@{ Not implemented.
+    OpalTrimCoil  (const OpalTrimCoil &) = delete;
+    void operator=(const OpalTrimCoil &) = delete;
+    ///@}
+    /// Private copy constructor, called by clone
+    OpalTrimCoil(const std::string &name, OpalTrimCoil *parent);
+
+};
+
+inline Inform &operator<<(Inform &os, const OpalTrimCoil &b) {
+    return b.print(os);
+}
+
+#endif // OPAL_TRIM_COIL_H

src/Classic/TrimCoils/TrimCoil.cpp

+#include "TrimCoil.h"
+
+TrimCoil::TrimCoil(){}
+
+void TrimCoil::applyField(const double r, const double z, double *br, double *bz)
+{
+    doApplyField(r,z,br,bz);
+}
\ No newline at end of file

src/Classic/TrimCoils/TrimCoil.h

+#ifndef TRIM_COIL_H
+#define TRIM_COIL_H
+
+/// Abstract TrimCoil class
+
+class TrimCoil {
+
+public:
+
+    TrimCoil();
+    /// Apply the trim coil at position r and z to Bfields br and bz
+    /// Calls virtual doApplyField
+    void applyField(const double r, const double z, double *br, double *bz);
+
+    virtual ~TrimCoil() { };
+
+private:
+
+    /// virtual implementation of applyField
+    virtual void doApplyField(const double r, const double z, double *br, double *bz) = 0;
+};
+
+#endif //TRIM_COIL_H

src/Classic/TrimCoils/TrimCoilFit.cpp

+#include "TrimCoils/TrimCoilFit.h"
+
+#include <cmath>
+
+TrimCoilFit::TrimCoilFit(double bmax,
+                         double rmin,
+                         double rmax,
+                         const std::vector<double>& coefnum,
+                         const std::vector<double>& coefdenom):
+    TrimCoil(),
+    coefnum_m(coefnum),
+    coefdenom_m(coefdenom)
+{
+    // convert to m
+    const double mm2m = 0.001;
+    rmin_m = rmin * mm2m;
+    rmax_m = rmax * mm2m;
+    // convert to kG
+    bmax_m   = bmax * 10.0;
+
+    // normal polynom if no denominator coefficients (denominator = 1)
+    if (coefdenom_m.empty())
+      coefdenom_m.push_back(1.0);
+}
+
+void TrimCoilFit::doApplyField(const double r, const double z, double *br, double *bz)
+{
+    if (std::abs(bmax_m) < 1e-20) return;
+    // check range
+    if (r < rmin_m || r > rmax_m) return;
+
+    double num     = 0.0; // numerator
+    double dnum_dr = 0.0; // derivative of numerator
+    double powr    = 1.0; // power of r
+
+    // add constant
+    num += coefnum_m[0];
+    for (std::size_t i = 1; i < coefnum_m.size(); ++i) {
+        dnum_dr += coefnum_m[i] * powr * i;
+        powr    *= r;
+        num     += coefnum_m[i] * powr;
+    }
+
+    double denom     = 0.0; // denominator
+    double ddenom_dr = 0.0; // derivative of denominator
+    powr             = 1.0; // power of r
+
+    // add constant
+    denom += coefdenom_m[0];
+    for (std::size_t i = 1; i < coefdenom_m.size(); ++i) {
+        ddenom_dr += coefdenom_m[i] * powr * i;
+        powr      *= r;
+        denom     += coefdenom_m[i] * powr;
+    }
+
+    double btr = num / denom;
+
+    // derivative of dr with quotient rule
+    double dr = (dnum_dr * denom - ddenom_dr * num) / (denom*denom);
+
+    //std::cout << "r " << r << " dr " <<  dr << std::endl;
+
+    *bz += bmax_m * btr;
+    *br += bmax_m * dr * z;
+}
\ No newline at end of file

src/Classic/TrimCoils/TrimCoilFit.h

+#ifndef TRIM_COILFIT_H
+#define TRIM_COILFIT_H
+
+#include "TrimCoils/TrimCoil.h"
+
+#include <vector>
+
+/// TrimCoilFit class
+/// General rational function fit
+/// https://gitlab.psi.ch/OPAL/src/issues/157
+
+class TrimCoilFit : public TrimCoil {
+
+public:
+    TrimCoilFit(double bmax,
+                double rmin,
+                double rmax,
+                const std::vector<double>& coefnum,
+                const std::vector<double>& coefdenom);
+
+    virtual ~TrimCoilFit() { };
+
+private:
+    TrimCoilFit() = delete;
+    /// Maximum B field (kG)
+    double bmax_m;
+    /// Minimum radius (m)
+    double rmin_m;
+    /// Maximum radius (m)
+    double rmax_m;
+
+    std::vector<double> coefnum_m;
+    std::vector<double> coefdenom_m;
+
+    /// @copydoc TrimCoil::doApplyField
+    virtual void doApplyField(const double r, const double z, double *br, double *bz);
+};
+
+#endif //TRIM_COILFIT_H

src/Classic/TrimCoils/TrimCoilMirrored.cpp

+#include "TrimCoils/TrimCoilMirrored.h"
+
+#include "Utility/IpplInfo.h"
+#include "Utilities/GeneralClassicException.h"
+
+#include <cmath>
+
+TrimCoilMirrored::TrimCoilMirrored(double bmax,
+                                   double rmin,
+                                   double rmax,
+                                   double bslope):
+    TrimCoil()
+{
+    // convert to m
+    const double mm2m = 0.001;
+    rmin_m   = rmin   * mm2m;
+    rmax_m   = rmax   * mm2m;
+    bslope_m = bslope / mm2m;
+    // convert to kG
+    bmax_m   = bmax * 10.0;
+}
+
+void TrimCoilMirrored::doApplyField(const double r, const double z, double *br, double *bz)
+{
+    /// update bz and br with trim coil contributions
+    // for some discussion on the formulas see
+    // http://doi.org/10.1103/PhysRevSTAB.14.054402
+    // https://gitlab.psi.ch/OPAL/src/issues/157
+    // https://gitlab.psi.ch/OPAL/src/issues/110
+
+    if (std::abs(bmax_m) < 1e-20) return;
+
+    // unitless constants
+    const double Amax1 = 1;
+    const double Amax2 = 3;
+    const double Amin  = -2;
+    const double x01   = 4;
+    const double x02   = 8;
+    const double h1    = 0.03;
+    const double h2    = 0.2;
+    const double justAnotherFudgeFactor = 1 / 2.78;
+    // helper variables
+    const double log10  = std::log(10);
+    const double const3 = -(Amax1 - Amin) * h1 * log10;
+    const double const4 =  (Amax2 - Amin) * h2 * log10;
+
+    const double &tcr1      = rmin_m;
+    const double &tcr2      = rmax_m;
+    const double &slope     = bslope_m;
+    const double &magnitude = bmax_m;
+
+    double part1;
+    double part2;
+
+    if (2 * r < (tcr2 + tcr1)) {
+        part1 = std::pow(10.0,  ((r - tcr1) * slope - x01) * h1);
+        part2 = std::pow(10.0, -((r - tcr1) * slope - x02) * h2);
+    } else {
+        part1 = std::pow(10.0,  ((tcr2 - r) * slope - x01) * h1);
+        part2 = std::pow(10.0, -((tcr2 - r) * slope - x02) * h2);
+    }
+
+    const double part1plusinv = 1 / (1 + part1);
+    const double part2plusinv = 1 / (1 + part2);
+
+    double part3 = const3 * slope * part1 * part1plusinv * part1plusinv;
+    double part4 = const4 * slope * part2 * part2plusinv * part2plusinv;
+
+    const double constmag = magnitude * justAnotherFudgeFactor;
+
+    double dr  = constmag * (part3 + part4);
+    double btr = constmag * (Amin - 1 +