Helmholtz derivatives are all calculated using parallel code - next step is to remove classes from Helmholtz

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>

dev/codelite/coolprop.project

@@ -195,6 +195,7 @@
         <IncludePath Value="../../include"/>
         <IncludePath Value="../../externals/Eigen"/>
         <Preprocessor Value="NDEBUG"/>
+        <Preprocessor Value="ENABLE_CATCH"/>
       </Compiler>
       <Linker Options="" Required="yes"/>
       <ResourceCompiler Options="" Required="no"/>

include/CoolPropTools.h

@@ -136,12 +136,12 @@
 	double interp1d(std::vector<double> *x, std::vector<double> *y, double x0);
 	double powInt(double x, int y);
     
-    #define POW2(x) x*x
-    #define POW3(x) x*x*x
-    #define POW4(x) x*x*x*x
-    #define POW5(x) x*x*x*x*x
-    #define POW6(x) x*x*x*x*x*x
-    #define POW7(x) x*x*x*x*x*x*x
+    #define POW2(x) (x)*(x)
+    #define POW3(x) (x)*(x)*(x)
+    #define POW4(x) (x)*(x)*(x)*(x)
+    #define POW5(x) (x)*(x)*(x)*(x)*(x)
+    #define POW6(x) (x)*(x)*(x)*(x)*(x)*(x)
+    #define POW7(x) (x)*(x)*(x)*(x)*(x)*(x)*(x)
 
     template<class T> T LinearInterp(T x0, T x1, T y0, T y1, T x)
     {

include/Helmholtz.h

@@ -100,7 +100,171 @@ public:
 
 struct Derivatives
 {
-    long double alphar, dalphar_ddelta, dalphar_dtau, d2alphar_ddelta2, d2alphar_dtau2, d2alphar_ddelta_dtau;
+    long double alphar, dalphar_ddelta, dalphar_dtau, d2alphar_ddelta2, d2alphar_dtau2, d2alphar_ddelta_dtau, 
+                d3alphar_ddelta3, d3alphar_ddelta_dtau2, d3alphar_ddelta2_dtau, d3alphar_dtau3;
+};
+
+struct ResidualHelmholtzGeneralizedExponentialElement
+{
+    /// These variables are for the n*delta^d_i*tau^t_i part
+    long double n,d,t;
+    /// These variables are for the exp(u) part
+    /// u is given by -c*delta^l_i-omega*tau^m_i-eta1*(delta-epsilon1)-eta2*(delta-epsilon2)^2-beta1*(tau-gamma1)-beta2*(tau-gamma2)^2
+    long double c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
+    /// If l_i or m_i are integers, we will store them as integers in order to call pow(double, int) rather than pow(double, double)
+    int l_int, m_int;
+    
+    ResidualHelmholtzGeneralizedExponentialElement()
+    {
+        n = _HUGE; d = _HUGE; t = _HUGE;
+        c = _HUGE; l_double = _HUGE; omega = _HUGE; m_double = _HUGE; 
+        eta1 = _HUGE; epsilon1 = _HUGE; eta2 = _HUGE; epsilon2 = _HUGE;
+        beta1 = _HUGE; gamma1 = _HUGE; beta2 = _HUGE; gamma2 = _HUGE;
+        l_int = -100000;
+        m_int = -100000;
+    }
+};
+class ResidualHelmholtzGeneralizedExponential : public BaseHelmholtzTerm{
+    
+public:
+    bool delta_li_in_u, tau_mi_in_u, eta1_in_u, eta2_in_u, beta1_in_u, beta2_in_u;
+    std::vector<long double> s;
+    std::size_t N;
+    std::vector<ResidualHelmholtzGeneralizedExponentialElement> elements;
+    // Default Constructor
+    ResidualHelmholtzGeneralizedExponential(){N = 0; 
+                                              delta_li_in_u = false; 
+                                              tau_mi_in_u = false;
+                                              eta1_in_u = false;
+                                              eta2_in_u = false;
+                                              beta1_in_u = false;
+                                              beta2_in_u = false;
+                                              };
+    
+    void add_Power(const std::vector<long double> &n, const std::vector<long double> &d, 
+                   const std::vector<long double> &t, const std::vector<long double> &l)
+    {
+        for (std::size_t i = 0; i < n.size(); ++i)
+        {
+            ResidualHelmholtzGeneralizedExponentialElement el;
+            el.n = n[i];
+            el.d = d[i];
+            el.t = t[i];
+            el.l_double = l[i];
+            el.l_int = (int)el.l_double;
+            if (el.l_double > 0)
+                el.c = 1.0;
+            else
+                el.c = 0.0;
+            elements.push_back(el);
+        }
+        delta_li_in_u = true;
+    };
+    void add_Exponential(const std::vector<long double> &n, const std::vector<long double> &d, 
+                         const std::vector<long double> &t, const std::vector<long double> &g, 
+                         const std::vector<long double> &l)
+    {
+        for (std::size_t i = 0; i < n.size(); ++i)
+        {
+            ResidualHelmholtzGeneralizedExponentialElement el;
+            el.n = n[i];
+            el.d = d[i];
+            el.t = t[i];
+            el.c = g[i];
+            el.l_double = l[i];
+            el.l_int = (int)el.l_double;
+            elements.push_back(el);
+        }
+        delta_li_in_u = true;
+    }   
+    void add_Gaussian(const std::vector<long double> &n, 
+                      const std::vector<long double> &d, 
+                      const std::vector<long double> &t, 
+                      const std::vector<long double> &eta, 
+                      const std::vector<long double> &epsilon,
+                      const std::vector<long double> &beta,
+                      const std::vector<long double> &gamma
+                      )
+    { 
+        for (std::size_t i = 0; i < n.size(); ++i)
+        {
+            ResidualHelmholtzGeneralizedExponentialElement el;
+            el.n = n[i];
+            el.d = d[i];
+            el.t = t[i];
+            el.eta2 = eta[i];
+            el.epsilon2 = epsilon[i];
+            el.beta2 = beta[i];
+            el.gamma2 = gamma[i];
+            elements.push_back(el);
+        }
+        eta2_in_u = true;
+        beta2_in_u = true;
+    };
+    void add_GERG2008Gaussian(const std::vector<long double> &n, 
+                              const std::vector<long double> &d, 
+                              const std::vector<long double> &t, 
+                              const std::vector<long double> &eta, 
+                              const std::vector<long double> &epsilon,
+                              const std::vector<long double> &beta,
+                              const std::vector<long double> &gamma)
+    { 
+        for (std::size_t i = 0; i < n.size(); ++i)
+        {
+            ResidualHelmholtzGeneralizedExponentialElement el;
+            el.n = n[i];
+            el.d = d[i];
+            el.t = t[i];
+            el.eta2 = eta[i];
+            el.epsilon2 = epsilon[i];
+            el.eta1 = beta[i];
+            el.epsilon1 = gamma[i];
+            elements.push_back(el);
+        }
+        eta2_in_u = true;
+        eta1_in_u = true;
+    };
+    void add_Lemmon2005(const std::vector<long double> &n, 
+                        const std::vector<long double> &d, 
+                        const std::vector<long double> &t, 
+                        const std::vector<long double> &l, 
+                        const std::vector<long double> &m)
+    {
+        for (std::size_t i = 0; i < n.size(); ++i)
+        {
+            ResidualHelmholtzGeneralizedExponentialElement el;
+            el.n = n[i];
+            el.d = d[i];
+            el.t = t[i];
+            el.c = 1.0;
+            el.omega = 1.0;
+            el.l_double = l[i];
+            el.m_double = m[i];
+            el.l_int = (int)el.l_double;
+            el.m_int = (int)el.m_double;
+            elements.push_back(el);
+        }
+        delta_li_in_u = true;
+        tau_mi_in_u = true;
+    };
+
+    ///< Destructor for the class.  No implementation
+    ~ResidualHelmholtzGeneralizedExponential(){};
+
+    void to_json(rapidjson::Value &el, rapidjson::Document &doc);
+    
+    long double base(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dDelta(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dTau(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dDelta2(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dDelta_dTau(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dTau2(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dDelta3(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dDelta2_dTau(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dDelta_dTau2(const long double &tau, const long double &delta) throw(){return 0;};
+    long double dTau3(const long double &tau, const long double &delta) throw(){return 0;};
+    
+    void all(const long double &tau, const long double &delta, Derivatives &derivs) throw();
 };
 
 struct ResidualHelmholtzPowerElement
@@ -108,6 +272,8 @@ struct ResidualHelmholtzPowerElement
     long double n,d,t,ld;
     int l;
 };
+
+
 /// Power term
 /*!
 Term of the form 
@@ -144,7 +310,7 @@ public:
             elements.push_back(el);
         }
     };
-
+    
     ///< Destructor for the alphar_power class.  No implementation
     ~ResidualHelmholtzPower(){};
 
@@ -160,8 +326,6 @@ public:
     long double dDelta2_dTau(const long double &tau, const long double &delta) throw();
     long double dDelta_dTau2(const long double &tau, const long double &delta) throw();
     long double dTau3(const long double &tau, const long double &delta) throw();
-    
-    void all(const long double &tau, const long double &delta, Derivatives &derivs) throw();
 };
 
 struct ResidualHelmholtzExponentialElement
@@ -497,79 +661,73 @@ class ResidualHelmholtzContainer
 {
     
 public:
-    ResidualHelmholtzPower Power;
-    ResidualHelmholtzExponential Exponential;
-    ResidualHelmholtzGaussian Gaussian;
-    ResidualHelmholtzLemmon2005 Lemmon2005;
     ResidualHelmholtzNonAnalytic NonAnalytic;
     ResidualHelmholtzSAFTAssociating SAFT;
+    ResidualHelmholtzGeneralizedExponential GenExp; 
 
     long double base(long double tau, long double delta)
     {
-        return (Power.base(tau, delta) + Exponential.base(tau, delta)
-                +Gaussian.base(tau, delta) + Lemmon2005.base(tau, delta)
-                +NonAnalytic.base(tau, delta) + SAFT.base(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.alphar + NonAnalytic.base(tau, delta) + SAFT.base(tau,delta));
     };
-
     long double dDelta(long double tau, long double delta)
     {
-        return (Power.dDelta(tau, delta) + Exponential.dDelta(tau, delta) 
-                + Gaussian.dDelta(tau, delta) + Lemmon2005.dDelta(tau, delta) 
-                + NonAnalytic.dDelta(tau, delta) + SAFT.dDelta(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.dalphar_ddelta + NonAnalytic.dDelta(tau, delta) + SAFT.dDelta(tau,delta));
     };
     long double dTau(long double tau, long double delta)
     {
-        return (Power.dTau(tau, delta) + Exponential.dTau(tau, delta) 
-                + Gaussian.dTau(tau, delta) + Lemmon2005.dTau(tau, delta) 
-                + NonAnalytic.dTau(tau, delta) + SAFT.dTau(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.dalphar_dtau + NonAnalytic.dTau(tau, delta) + SAFT.dTau(tau,delta));
     };
-
     long double dDelta2(long double tau, long double delta)
     {
-        return (Power.dDelta2(tau, delta) + Exponential.dDelta2(tau, delta)
-                +Gaussian.dDelta2(tau, delta) + Lemmon2005.dDelta2(tau, delta)
-                +NonAnalytic.dDelta2(tau, delta) + SAFT.dDelta2(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.d2alphar_ddelta2 + NonAnalytic.dDelta2(tau, delta) + SAFT.dDelta2(tau,delta));
     };
     long double dDelta_dTau(long double tau, long double delta)
     {
-        return (Power.dDelta_dTau(tau, delta) + Exponential.dDelta_dTau(tau, delta) 
-                + Gaussian.dDelta_dTau(tau, delta) + Lemmon2005.dDelta_dTau(tau, delta) 
-                + NonAnalytic.dDelta_dTau(tau, delta) + SAFT.dDelta_dTau(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.d2alphar_ddelta_dtau + NonAnalytic.dDelta_dTau(tau, delta) + SAFT.dDelta_dTau(tau,delta));
     };
     long double dTau2(long double tau, long double delta)
     {
-        return (Power.dTau2(tau, delta) + Exponential.dTau2(tau, delta) 
-                + Gaussian.dTau2(tau, delta) + Lemmon2005.dTau2(tau, delta) 
-                + NonAnalytic.dTau2(tau, delta) + SAFT.dTau2(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.d2alphar_dtau2 + NonAnalytic.dTau2(tau, delta) + SAFT.dTau2(tau,delta));
     };
     
     long double dDelta3(long double tau, long double delta) 
     {
-        return (Power.dDelta3(tau, delta) + Exponential.dDelta3(tau, delta)
-                +Gaussian.dDelta3(tau, delta) + Lemmon2005.dDelta3(tau, delta)
-                +NonAnalytic.dDelta3(tau, delta) + SAFT.dDelta3(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.d3alphar_ddelta3 + NonAnalytic.dDelta3(tau, delta) + SAFT.dDelta3(tau,delta));
     };
     long double dDelta2_dTau(long double tau, long double delta)
     {
-        return (Power.dDelta2_dTau(tau, delta) + Exponential.dDelta2_dTau(tau, delta) 
-                + Gaussian.dDelta2_dTau(tau, delta) + Lemmon2005.dDelta2_dTau(tau, delta) 
-                + NonAnalytic.dDelta2_dTau(tau, delta) + SAFT.dDelta2_dTau(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.d3alphar_ddelta2_dtau + NonAnalytic.dDelta2_dTau(tau, delta) + SAFT.dDelta2_dTau(tau,delta));
     };
     long double dDelta_dTau2(long double tau, long double delta)
     {
-        return (Power.dDelta_dTau2(tau, delta) + Exponential.dDelta_dTau2(tau, delta) 
-                + Gaussian.dDelta_dTau2(tau, delta) + Lemmon2005.dDelta_dTau2(tau, delta) 
-                + NonAnalytic.dDelta_dTau2(tau, delta) + SAFT.dDelta_dTau2(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.d3alphar_ddelta_dtau2 + NonAnalytic.dDelta_dTau2(tau, delta) + SAFT.dDelta_dTau2(tau,delta));
     };
     long double dTau3(long double tau, long double delta)
     {
-        return (Power.dTau3(tau, delta) + Exponential.dTau3(tau, delta)
-                +Gaussian.dTau3(tau, delta) + Lemmon2005.dTau3(tau, delta)
-                +NonAnalytic.dTau3(tau, delta) + SAFT.dTau3(tau,delta));
+        Derivatives derivs;
+        GenExp.all(tau, delta, derivs);
+        return (derivs.d3alphar_dtau3 + NonAnalytic.dTau3(tau, delta) + SAFT.dTau3(tau,delta));
     };
 };
 
-
 // #############################################################################
 // #############################################################################
 // #############################################################################

src/Backends/Helmholtz/Fluids/FluidLibrary.h

@@ -42,7 +42,6 @@ protected:
 
             if (!type.compare("ResidualHelmholtzPower"))
             {
-                if (EOS.alphar.Power.N > 0){throw ValueError("Cannot add ");}
                 std::vector<long double> n = cpjson::get_long_double_array(contribution["n"]);
                 std::vector<long double> d = cpjson::get_long_double_array(contribution["d"]);
                 std::vector<long double> t = cpjson::get_long_double_array(contribution["t"]);
@@ -50,11 +49,10 @@ protected:
                 assert(n.size() == d.size());
                 assert(n.size() == t.size());
                 assert(n.size() == l.size());
-                EOS.alphar.Power = ResidualHelmholtzPower(n,d,t,l);
+                EOS.alphar.GenExp.add_Power(n,d,t,l);
             }
             else if (!type.compare("ResidualHelmholtzGaussian"))
             {
-                if (EOS.alphar.Gaussian.N > 0){throw ValueError("Cannot add ");}
                 std::vector<long double> n = cpjson::get_long_double_array(contribution["n"]);
                 std::vector<long double> d = cpjson::get_long_double_array(contribution["d"]);
                 std::vector<long double> t = cpjson::get_long_double_array(contribution["t"]);
@@ -68,7 +66,7 @@ protected:
                 assert(n.size() == epsilon.size());
                 assert(n.size() == beta.size());
                 assert(n.size() == gamma.size());
-                EOS.alphar.Gaussian = ResidualHelmholtzGaussian(n,d,t,eta,epsilon,beta,gamma);
+                EOS.alphar.GenExp.add_Gaussian(n,d,t,eta,epsilon,beta,gamma);
             }
             else if (!type.compare("ResidualHelmholtzNonAnalytic"))
             {
@@ -92,7 +90,6 @@ protected:
             }
             else if (!type.compare("ResidualHelmholtzLemmon2005"))
             {
-                if (EOS.alphar.Lemmon2005.N > 0){throw ValueError("Cannot add ");}
                 std::vector<long double> n = cpjson::get_long_double_array(contribution["n"]);
                 std::vector<long double> d = cpjson::get_long_double_array(contribution["d"]);
                 std::vector<long double> t = cpjson::get_long_double_array(contribution["t"]);
@@ -102,11 +99,10 @@ protected:
                 assert(n.size() == t.size());
                 assert(n.size() == l.size());
                 assert(n.size() == m.size());
-                EOS.alphar.Lemmon2005 = ResidualHelmholtzLemmon2005(n,d,t,l,m);
+                EOS.alphar.GenExp.add_Lemmon2005(n,d,t,l,m);
             }
             else if (!type.compare("ResidualHelmholtzExponential"))
             {
-                if (EOS.alphar.Exponential.N > 0){throw ValueError("Cannot add ");}
                 std::vector<long double> n = cpjson::get_long_double_array(contribution["n"]);
                 std::vector<long double> d = cpjson::get_long_double_array(contribution["d"]);
                 std::vector<long double> t = cpjson::get_long_double_array(contribution["t"]);
@@ -116,7 +112,7 @@ protected:
                 assert(n.size() == t.size());
                 assert(n.size() == g.size());
                 assert(n.size() == l.size());
-                EOS.alphar.Exponential = ResidualHelmholtzExponential(n,d,t,g,l);
+                EOS.alphar.GenExp.add_Exponential(n,d,t,g,l);
             }
             else if (!type.compare("ResidualHelmholtzAssociating"))
             {

src/Helmholtz.cpp

@@ -35,56 +35,143 @@ long double kahanSum(std::vector<long double> &x)
 }
 bool wayToSort(long double i, long double j) { return std::abs(i) > std::abs(j); }
 
-void ResidualHelmholtzPower::all(const long double &tau, const long double &delta, Derivatives &derivs) throw()
+void ResidualHelmholtzGeneralizedExponential::all(const long double &tau, const long double &delta, Derivatives &derivs) throw()
 {
-    long double log_tau = log(tau), log_delta = log(delta), u, du_ddelta, du_dtau, d2u_ddelta2, d2u_dtau2, ndteu, 
+    long double log_tau = log(tau), log_delta = log(delta), ndteu, 
+                u, du_ddelta, du_dtau, d2u_ddelta2, d2u_dtau2, d3u_ddelta3, d3u_dtau3,
                 one_over_delta = 1/delta, one_over_tau = 1/tau, // division is much slower than multiplication, so do one division here
-                B_delta, B_tau, B_delta2, B_tau2;
+                B_delta, B_tau, B_delta2, B_tau2, B_delta3, B_tau3,
+                u_increment, du_ddelta_increment, du_dtau_increment, 
+                d2u_ddelta2_increment, d3u_ddelta3_increment, d2u_dtau2_increment, d3u_dtau3_increment;
+                
     derivs.alphar = 0.0;
     derivs.dalphar_ddelta = 0.0;
     derivs.dalphar_dtau = 0.0;
     derivs.d2alphar_ddelta2 = 0.0;
     derivs.d2alphar_dtau2 = 0.0;
     derivs.d2alphar_ddelta_dtau = 0.0;
+    derivs.d3alphar_dtau3 = 0.0;
+    derivs.d3alphar_ddelta_dtau2 = 0.0;
+    derivs.d3alphar_ddelta2_dtau = 0.0;
+    derivs.d3alphar_ddelta3 = 0.0;
+    
+    const std::size_t N = elements.size();
+    
+    // We define them locally as constants so that the compiler 
+    // can more intelligently branch the for loop
+    const bool delta_li_in_u = this->delta_li_in_u;
+    const bool tau_mi_in_u = this->tau_mi_in_u;
+    const bool eta2_in_u = this->eta2_in_u;
+    const bool beta2_in_u = this->beta2_in_u;
+    const bool eta1_in_u = this->eta1_in_u;
+    const bool beta1_in_u = this->beta1_in_u;
     
     for (std::size_t i = 0; i < N; ++i)
     {
-        ResidualHelmholtzPowerElement &el = elements[i];
+        ResidualHelmholtzGeneralizedExponentialElement &el = elements[i];
         long double ni = el.n, di = el.d, ti = el.t;
-        int li = el.l;
-        if (li > 0){
-            u = -pow(delta, li);
-            du_ddelta = li*u*one_over_delta;
-            du_dtau = 0;
-            d2u_ddelta2 = (li-1)*du_ddelta*one_over_delta;
-            d2u_dtau2 = 0;
+        
+        // Set the u part of exp(u) to zero
+        u = 0;
+        du_ddelta = 0;
+        du_dtau = 0;
+        d2u_ddelta2 = 0;
+        d2u_dtau2 = 0;
+        d3u_ddelta3 = 0;
+        d3u_dtau3 = 0;
+        
+        if (delta_li_in_u){
+            long double  ci = el.c, l_double = el.l_double;
+            int l_int = el.l_int;
+            if (ValidNumber(l_double) && l_int > 0){
+                u_increment = -ci*pow(delta, l_int);
+                du_ddelta_increment = l_double*u_increment*one_over_delta;
+                d2u_ddelta2_increment = (l_double-1)*du_ddelta_increment*one_over_delta;
+                d3u_ddelta3_increment = (l_double-2)*d2u_ddelta2_increment*one_over_delta;
+                u += u_increment;
+                du_ddelta += du_ddelta_increment;
+                d2u_ddelta2 += d2u_ddelta2_increment;
+                d3u_ddelta3 += d3u_ddelta3_increment;
+            }
         }
-        else{
-            u = 0;
-            du_ddelta = 0;
-            du_dtau = 0;
-            d2u_ddelta2 = 0;
-            d2u_dtau2 = 0;
+        if (tau_mi_in_u){
+            long double omegai = el.omega, m_double = el.m_double;
+            if (std::abs(m_double) > 0){
+                u_increment = -omegai*pow(tau, m_double);
+                du_dtau_increment = m_double*u_increment*one_over_tau;
+                d2u_dtau2_increment = (m_double-1)*du_dtau_increment*one_over_tau;
+                d3u_dtau3_increment = (m_double-2)*d2u_dtau2_increment*one_over_tau;
+                u += u_increment;
+                du_dtau += du_dtau_increment;
+                d2u_dtau2 += d2u_dtau2_increment;
+                d3u_dtau3 += d3u_dtau3_increment;
+            }
         }
+        if (eta1_in_u){
+            long double eta1 = el.eta1, epsilon1 = el.epsilon1;
+            if (ValidNumber(eta1)){
+                u += -eta1*(delta-epsilon1);
+                du_ddelta += -eta1;
+            }
+        }
+        if (eta2_in_u){
+            long double eta2 = el.eta2, epsilon2 = el.epsilon2;
+            if (ValidNumber(eta2)){
+                u += -eta2*POW2(delta-epsilon2);
+                du_ddelta += -2*eta2*(delta-epsilon2);
+                d2u_ddelta2 += -2*eta2;
+            }
+        }
+        if (beta1_in_u){
+            long double beta1 = el.beta1, gamma1 = el.gamma1;
+            if (ValidNumber(beta1)){
+                u += -beta1*(tau-gamma1);
+                du_dtau += -beta1;
+            }
+        }
+        if (beta2_in_u){
+            long double beta2 = el.beta2, gamma2 = el.gamma2;
+            if (ValidNumber(beta2)){
+                u += -beta2*POW2(tau-gamma2);
+                du_dtau += -2*beta2*(tau-gamma2);
+                d2u_dtau2 += -2*beta2;
+            }
+        }
+        
         ndteu = ni*exp(ti*log_tau+di*log_delta+u);
         
         B_delta = (delta*du_ddelta + di);
         B_tau = (tau*du_dtau + ti);
         B_delta2 = (POW2(delta)*(d2u_ddelta2 + POW2(du_ddelta)) + 2*di*delta*du_ddelta + di*(di-1));
         B_tau2 = (POW2(tau)*(d2u_dtau2 + POW2(du_dtau)) + 2*ti*tau*du_dtau + ti*(ti-1));
+        B_delta3 = POW3(delta)*d3u_ddelta3 + 3*di*POW2(delta)*d2u_ddelta2+3*POW3(delta)*d2u_ddelta2*du_ddelta+3*di*POW2(delta*du_ddelta)+3*di*(di-1)*delta*du_ddelta+di*(di-1)*(di-2)+POW3(delta*du_ddelta);
+        B_tau3 = POW3(tau)*d3u_dtau3 + 3*ti*POW2(tau)*d2u_dtau2+3*POW3(tau)*d2u_dtau2*du_dtau+3*ti*POW2(tau*du_dtau)+3*ti*(ti-1)*tau*du_dtau+ti*(ti-1)*(ti-2)+POW3(tau*du_dtau);
         
         derivs.alphar += ndteu;
+        
         derivs.dalphar_ddelta += ndteu*B_delta;
         derivs.dalphar_dtau += ndteu*B_tau;
+        
         derivs.d2alphar_ddelta2 += ndteu*B_delta2;
-        derivs.d2alphar_dtau2 += ndteu*B_tau2;
         derivs.d2alphar_ddelta_dtau += ndteu*B_delta*B_tau;
+        derivs.d2alphar_dtau2 += ndteu*B_tau2;
+        
+        derivs.d3alphar_ddelta3 += ndteu*B_delta3;
+        derivs.d3alphar_ddelta2_dtau += ndteu*B_delta2*B_tau;
+        derivs.d3alphar_ddelta_dtau2 += ndteu*B_delta*B_tau2;
+        derivs.d3alphar_dtau3 += ndteu*B_tau3;
     }
-    derivs.dalphar_ddelta       *= one_over_delta;
-    derivs.dalphar_dtau         *= one_over_tau;
-    derivs.d2alphar_ddelta2     *= POW2(one_over_delta);
-    derivs.d2alphar_dtau2       *= POW2(one_over_tau);
-    derivs.d2alphar_ddelta_dtau *= one_over_delta*one_over_tau;
+    derivs.dalphar_ddelta        *= one_over_delta;
+    derivs.dalphar_dtau          *= one_over_tau;
+    derivs.d2alphar_ddelta2      *= POW2(one_over_delta);
+    derivs.d2alphar_dtau2        *= POW2(one_over_tau);
+    derivs.d2alphar_ddelta_dtau  *= one_over_delta*one_over_tau;
+    
+    derivs.d3alphar_ddelta3      *= POW3(one_over_delta);
+    derivs.d3alphar_dtau3        *= POW3(one_over_tau);
+    derivs.d3alphar_ddelta2_dtau *= POW2(one_over_delta)*one_over_tau;
+    derivs.d3alphar_ddelta_dtau2 *= one_over_delta*POW2(one_over_tau);
+    
     return;
 };