From e4c6e83728dcfe5a91bc3546e071be23c77063e7 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 10 Jul 2014 15:07:12 +0200 Subject: [PATCH] Fixed Helmholtz dTau2 term for NonAnalytic --- src/Helmholtz.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp index 0bef5e86..e45a7eca 100644 --- a/src/Helmholtz.cpp +++ b/src/Helmholtz.cpp @@ -1097,15 +1097,16 @@ long double ResidualHelmholtzNonAnalytic::dTau2(const long double &tau, const lo ResidualHelmholtzNonAnalyticElement &el = elements[i]; long double ni = el.n, ai = el.a, bi = el.b, betai = el.beta; long double Ai = el.A, Bi = el.B, Ci = el.C, Di = el.D; - long double theta=(1.0-tau)+Ai*pow(pow(delta-1.0,2),1.0/(2.0*betai)); - long double DELTA=pow(theta,2)+Bi*pow(pow(delta-1.0,2),ai); - long double PSI=exp(-Ci*pow(delta-1.0,2)-Di*pow(tau-1.0,2)); - long double dPSI_dTau=-2.0*Di*(tau-1.0)*PSI; - long double dDELTAbi_dTau=-2.0*theta*bi*pow(DELTA,bi-1.0); - long double dPSI2_dTau2=(2.0*Di*pow(tau-1.0,2)-1.0)*2.0*Di*PSI; - long double dDELTAbi2_dTau2=2.0*bi*pow(DELTA,bi-1.0)+4.0*pow(theta,2)*bi*(bi-1.0)*pow(DELTA,bi-2.0); - s[i] = ni*(dDELTAbi2_dTau2*PSI+2.0*dDELTAbi_dTau*dPSI_dTau+pow(DELTA,bi)*dPSI2_dTau2); + long double theta = (1.0-tau)+Ai*pow(pow(delta-1.0,2),1.0/(2.0*betai)); + long double DELTA = pow(theta,2)+Bi*pow(pow(delta-1.0,2),ai); + long double PSI = exp(-Ci*pow(delta-1.0,2)-Di*pow(tau-1.0,2)); + long double dPSI_dTau = -2.0*Di*(tau-1.0)*PSI; + long double dDELTAbi_dTau = -2.0*theta*bi*pow(DELTA,bi-1.0); + long double dPSI2_dTau2 = (2.0*Di*pow(tau-1.0,2)-1.0)*2.0*Di*PSI; + long double dDELTAbi2_dTau2 = 2.0*bi*pow(DELTA,bi-1.0)+4.0*pow(theta,2)*bi*(bi-1.0)*pow(DELTA,bi-2.0); + + s[i] = ni*delta*(dDELTAbi2_dTau2*PSI+2.0*dDELTAbi_dTau*dPSI_dTau+pow(DELTA,bi)*dPSI2_dTau2); } return std::accumulate(s.begin(), s.end(), 0.0); }