exact_derivatives-Hreal.c

Go to the documentation of this file.

#define Power(A,B) pow(A,B)
#define Sqrt(A) sqrt(A)
#define Log(A) log(A)
const double csi1  =  1.0 + (1.-fabs(1.-tortoise)) * (csi - 1.0);
const double csi2  =  1.0 + (0.5-copysign(0.5, 1.5-tortoise)) * (csi - 1.0);
const double prT  =  csi2*(nx*p->data[0] + ny*p->data[1] + nz*p->data[2]);
const double prTtimesoneminuscsi1inv  =  (1. - 1./csi1)*prT;
const double tmpP0  =  p->data[0] - nx*prTtimesoneminuscsi1inv;
const double tmpP1  =  p->data[1] - ny*prTtimesoneminuscsi1inv;
const double tmpP2  =  p->data[2] - nz*prTtimesoneminuscsi1inv;
const double pxir  =  r*(tmpP0*xiUSCOREx + tmpP1*xiUSCOREy + tmpP2*xiUSCOREz);
const double pvr  =  r*(tmpP0*vx + tmpP1*vy + tmpP2*vz);
const double pvrsq  =  ((pvr)*(pvr));
const double pn  =  nx*tmpP0 + ny*tmpP1 + nz*tmpP2;
const double pnsq  =  ((pn)*(pn));
const double pr  =  pn;
const double prsq  =  ((pr)*(pr));
const double pf  =  pxir;
const double pxirsq  =  ((pxir)*(pxir));
const double ptheta2  =  invxi2*pvrsq;
const double prT4  =  ((prT)*(prT)*(prT)*(prT));
const double Hnspt7  =  deltaR*invrho2;
const double Hnspt6  =  invxi2*rho2invLambda;
const double Hnspt5  =  qq*u2;
const double Hnspt4  =  1. + Hnspt6*((pf)*(pf)) + Hnspt7*prsq + Hnspt5*prT4 + invrho2*ptheta2;
const double Hnspt3  =  deltaT*Hnspt4;
const double Hnspt2  =  Hnspt3*rho2;
const double Hnspt1  =  pf*ww;
const double Hns  =  Hnspt1*invLambda + Sqrt(Hnspt2*invLambda);
const double Qpt3  =  deltaR*invrho2;
const double Qpt2  =  invxi2*rho2invLambda;
const double Qpt1  =  invrho2*invxi2;
const double Q  =  1. + pvrsq*Qpt1 + pxirsq*Qpt2 + pnsq*Qpt3;
const double pn2  =  deltaR*invrho2*prsq;
const double pp  =  -1. + Q;
const double sKerrmultfact  =  -8. - 36.*pn2*r + 3.*pp*r;
const double sStarmultfact  =  14. - 30.*pn2*r + 4.*pp*r;
const double deltaSigmaStarUSCOREx1  =  etaover12r*(sKerrmultfact*sKerrUSCOREx + sStarmultfact*sStarUSCOREx);
const double deltaSigmaStarUSCOREy1  =  etaover12r*(sKerrmultfact*sKerrUSCOREy + sStarmultfact*sStarUSCOREy);
const double deltaSigmaStarUSCOREz1  =  etaover12r*(sKerrmultfact*sKerrUSCOREz + sStarmultfact*sStarUSCOREz);
const double pn2pp  =  pn2*pp;
const double pp2  =  ((pp)*(pp));
const double pn2u2  =  pn2*u2;
const double ppu2  =  pp*u2;
const double pn2ppu2  =  pn2pp*u2;
const double sMultiplier1pt6  =  -360.*((pn2)*(pn2)) + 126.*pn2pp + 3.*pp2;
const double sMultiplier1pt5  =  -96.*pn2pp + 23.*pp2;
const double sMultiplier1pt4  =  324.*pn2 - 120.*pp + r*sMultiplier1pt6;
const double sMultiplier1pt3  =  -282.*pn2 + 206.*pp + r*sMultiplier1pt5;
const double sMultiplier1pt2  =  54. + r*sMultiplier1pt4;
const double sMultiplier1pt1  =  -706. + eta*sMultiplier1pt2 + r*sMultiplier1pt3;
const double sMultiplier1  =  -0.013888888888888888*eta*sMultiplier1pt1*u2;
const double sMultiplier2pt6  =  -1.625*pn2ppu2 + 5.625*pn2*pn2u2;
const double sMultiplier2pt5  =  0.25*pn2ppu2 - 0.3125*pp2*u2;
const double sMultiplier2pt4  =  -6.125*pn2u2 + 1.4166666666666665*ppu2 + r*sMultiplier2pt6;
const double sMultiplier2pt3  =  -0.6666666666666666*pn2u2 - 3.0277777777777777*ppu2 + r*sMultiplier2pt5;
const double sMultiplier2pt2  =  r*sMultiplier2pt4 - 2.333333333333333*u2;
const double sMultiplier2pt1  =  eta*sMultiplier2pt2 + r*sMultiplier2pt3 - 6.222222222222221*u2;
const double sMultiplier2  =  eta*sMultiplier2pt1;
const double deltaSigmaStarUSCOREx2  =  deltaSigmaStarUSCOREx1 + sigmaStar->data[0]*sMultiplier1 + sigmaKerr->data[0]*sMultiplier2;
const double deltaSigmaStarUSCOREy2  =  deltaSigmaStarUSCOREy1 + sigmaStar->data[1]*sMultiplier1 + sigmaKerr->data[1]*sMultiplier2;
const double deltaSigmaStarUSCOREz2  =  deltaSigmaStarUSCOREz1 + sigmaStar->data[2]*sMultiplier1 + sigmaKerr->data[2]*sMultiplier2;
const double deltaSigmaStarUSCOREx3  =  deltaSigmaStarUSCOREx2 + coeffs->d1*etau3*sigmaStar->data[0];
const double deltaSigmaStarUSCOREy3  =  deltaSigmaStarUSCOREy2 + coeffs->d1*etau3*sigmaStar->data[1];
const double deltaSigmaStarUSCOREz3  =  deltaSigmaStarUSCOREz2 + coeffs->d1*etau3*sigmaStar->data[2];
const double deltaSigmaStarUSCOREx  =  deltaSigmaStarUSCOREx3 + coeffs->d1v2*etau3*sigmaKerr->data[0];
const double deltaSigmaStarUSCOREy  =  deltaSigmaStarUSCOREy3 + coeffs->d1v2*etau3*sigmaKerr->data[1];
const double deltaSigmaStarUSCOREz  =  deltaSigmaStarUSCOREz3 + coeffs->d1v2*etau3*sigmaKerr->data[2];
const double sx  =  deltaSigmaStarUSCOREx + sStarUSCOREx;
const double sy  =  deltaSigmaStarUSCOREy + sStarUSCOREy;
const double sz  =  deltaSigmaStarUSCOREz + sStarUSCOREz;
const double sxi  =  sx*xiUSCOREx + sy*xiUSCOREy + sz*xiUSCOREz;
const double sv  =  sx*vx + sy*vy + sz*vz;
const double sn  =  nx*sx + ny*sy + nz*sz;
const double s3  =  e3USCOREx*sx + e3USCOREy*sy + e3USCOREz*sz;
const double sqrtQ  =  Sqrt(Q);
const double oneplus2sqrtQ  =  1. + 2.*sqrtQ;
const double oneplus1sqrtQ  =  oneplus2sqrtQ - sqrtQ;
const double twoB1psqrtQsqrtQ  =  2.*B*oneplus1sqrtQ*sqrtQ;
const double invtwoB1psqrtQsqrtQ  =  1./twoB1psqrtQsqrtQ;
const double expMUsqsqrtQplusQ  =  expMUsq*(Q + sqrtQ);
const double Hwrpt4a  =  pxirsq*sv;
const double Hwrpt4  =  expMUsqexpnusq*Hwrpt4a;
const double Hwrpt3c  =  pxir*sxi;
const double Hwrpt3b  =  Hwrpt3c*pvr;
const double Hwrpt3a  =  expMUexpnu*Hwrpt3b;
const double Hwrpt3  =  B*Hwrpt3a;
const double Hwrpt2g  =  deltaR*sv;
const double Hwrpt2f  =  sn*sqrtdeltaR;
const double Hwrpt2e  =  Hwrpt2f*pvr;
const double Hwrpt2d  =  Hwrpt2g*pnsq;
const double Hwrpt2c  =  Hwrpt2e*pn;
const double Hwrpt2b  =  expMUsqsqrtQplusQ*sv;
const double Hwrpt2a  =  (Hwrpt2b + Hwrpt2c - Hwrpt2d)*xi2;
const double Hwrpt2  =  deltaT*Hwrpt2a;
const double Hwrpt1b  =  invtwoB1psqrtQsqrtQ*invxi2;
const double Hwrpt1a  =  Hwrpt1b*sqrtdeltaR;
const double Hwrpt1  =  Hwrpt1a*invexpMUcubinvexpnu;
const double Hwr  =  Hwrpt1*(Hwrpt2 - Hwrpt3 + Hwrpt4);
const double Hwcospt9  =  pxir*sxi;
const double Hwcospt8  =  pvr*sv;
const double Hwcospt7  =  B*Hwcospt8 - expMUexpnu*Hwcospt9;
const double Hwcospt6  =  Hwcospt7*sqrtdeltaR;
const double Hwcospt5  =  pvrsq - expMUsqsqrtQplusQ*xi2;
const double Hwcospt4  =  Hwcospt6*pn;
const double Hwcospt3  =  deltaT*Hwcospt5 - expMUsqexpnusq*pxirsq;
const double Hwcospt2  =  -(B*Hwcospt4) + Hwcospt3*sn;
const double Hwcospt1  =  Hwcospt2*invexpMUcubinvexpnu;
const double Hwcos  =  Hwcospt1*invtwoB1psqrtQsqrtQ;
const double deltaTsqrtQ  =  deltaT*sqrtQ;
const double invdeltatTsqrtQ  =  1./deltaTsqrtQ;
const double HSOLpt5  =  (-B + expMUexpnu)*pxir;
const double HSOLpt4  =  HSOLpt5*invexpMU;
const double HSOLpt3  =  expnusq*HSOLpt4;
const double HSOLpt2  =  HSOLpt3*s3;
const double HSOLpt1  =  HSOLpt2*invxi2;
const double HSOL  =  HSOLpt1*invdeltatTsqrtQ;
const double deltaTsqrtQplusQ  =  deltaT*(Q + sqrtQ);
const double invdeltaTsqrtQplusQ  =  1./deltaTsqrtQplusQ;
const double HSONLmult2  =  invdeltaTsqrtQplusQ*invxi2;
const double HSONLmult  =  expnuinvexpMU2*HSONLmult2;
const double HSONLpt1b  =  pn*xi2;
const double HSONLpt1a  =  HSONLpt1b*(-mucos + nucos) + (mur - nur)*pvr;
const double HSONLpt1  =  -(HSONLpt1b*mucos) + mur*pvr + HSONLpt1a*sqrtQ;
const double HSONLpt2d  =  nur*pxir;
const double HSONLpt2c  =  HSONLpt2d*oneplus2sqrtQ;
const double HSONLpt2b  =  B*sxi;
const double HSONLpt2a  =  expMUexpnu*HSONLpt2c;
const double HSONLpt2  =  HSONLpt1*HSONLpt2b + HSONLpt2a*sv;
const double HSONLpt3c  =  pxir*sv;
const double HSONLpt3b  =  HSONLpt3c*oneplus1sqrtQ;
const double HSONLpt3a  =  expMUexpnu*HSONLpt3b;
const double HSONLpt3  =  B*HSONLpt2 - BR*HSONLpt3a;
const double HSONLpt4e  =  sn*xi2;
const double HSONLpt4d  =  HSONLpt4e*oneplus2sqrtQ;
const double HSONLpt4c  =  HSONLpt4d*pxir;
const double HSONLpt4b  =  HSONLpt4c*nucos;
const double HSONLpt4a  =  expMUexpnu*HSONLpt4b;
const double HSONLpt4  =  -(B*HSONLpt4a) + HSONLpt3*sqrtdeltaR;
const double HSONL  =  HSONLmult*HSONLpt4;
const double Hs  =  HSOL + HSONL + s3*w + Hwcos*wcos + Hwr*wr;
const double Hsspt1  =  -0.5*(-3.*((sn)*(sn)) + ((sx)*(sx)) + ((sy)*(sy)) + ((sz)*(sz)));
const double Hss  =  Hsspt1*u3;
const double sKerrdotsStar  =  sKerrUSCOREx*sStarUSCOREx + sKerrUSCOREy*sStarUSCOREy + sKerrUSCOREz*sStarUSCOREz;
const double Hpt1  =  etau4*(s1dots1 + s2dots2);
const double H  =  Hns + coeffs->dheffSSv2*Hpt1 + Hs + Hss + coeffs->dheffSS*Hpt1*sKerrdotsStar;
const double Hreal  =  Sqrt(1. + 2.*eta*(-1. + H));
const double invHreal  =  1.0/Hreal;