52 REAL,
PARAMETER ::
pi = 3.14159265358979323846
53 REAL,
PARAMETER ::
sqrt_two = 1.41421356237309504880
59 REAL,
DIMENSION(MAX_PL_COEFF),
PARAMETER &
64 INTERFACE legendrepolynomial
99 ELEMENTAL FUNCTION asinh(x)
RESULT(fx)
101 REAL,
INTENT(IN) :: x
103 fx = log(x+sqrt(x*x+1.0))
109 ELEMENTAL FUNCTION acosh(x)
RESULT(fx)
111 REAL,
INTENT(IN) :: x
113 fx = log(x+sqrt(x*x-1.0))
119 ELEMENTAL FUNCTION atanh(x)
RESULT(fx)
121 REAL,
INTENT(IN) :: x
123 fx = 0.5*log((1.0+x)/(1.0-x))
133 REAL,
INTENT(IN) :: x,a
135 fx = 0.5*(sign(1.0,x-a)+1.0)
144 REAL,
INTENT(IN) :: x,a,b
154 REAL,
INTENT(IN) :: u1,u2
155 REAL,
INTENT(OUT) :: x,y
159 r = sqrt(-2.0*log(u1))
194 REAL,
INTENT(IN) :: k
195 REAL,
INTENT(OUT) :: eell,kell
201 IF (abs(k).LT.1.0)
THEN
209 IF (abs(cn).GT.0.5*
eps_agm*abs(an))
THEN
214#if !defined(NECSXAURORA)
221 kell =
pi / (an + bn + tiny(kell))
223 kell = sqrt(-1.0*abs(k))
225 eell = (1.0-0.5*tmp)*kell
240 REAL,
INTENT(IN) :: k
246 IF (abs(k).LT.1.0)
THEN
251 IF (abs(an-bn).GT.
eps_agm*abs(an))
THEN
255#if !defined(NECSXAURORA)
263 kell =
pi / ( an + bn + tiny(kell) )
265 kell = sqrt(-1.0*abs(k))
277 REAL,
INTENT(IN) :: k
291 ELEMENTAL FUNCTION rf(x,y,z)
RESULT(res)
296 REAL,
PARAMETER :: errtol = .0025
297 REAL,
PARAMETER :: third = 1./3.
298 REAL,
PARAMETER :: tin = 1.5e-38
299 REAL,
PARAMETER :: big = 3.e37
300 REAL,
PARAMETER :: c1 = 1./24.
301 REAL,
PARAMETER :: c2 = 0.1
302 REAL,
PARAMETER :: c3 = 3./44.
303 REAL,
PARAMETER :: c4 = 1./14.
304 REAL :: alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,&
309 IF((min(x,y,z).LT.0.).OR.&
310 (min(x+y,x+z,y+z).LT.tin).OR.&
311 (max(x,y,z).GT.big))
THEN
312 res = sqrt(-1.*abs(x))
321 alamb = sqrtx * (sqrty+sqrtz) + sqrty*sqrtz
325 ave = third*(xt+yt+zt)
329 IF(max(abs(delx),abs(dely),abs(delz)).LE.errtol) &
333 e2 = delx*dely-delz**2
335 res = (1. + (c1*e2-c2-c3*e3)*e2 + c4*e3)/sqrt(ave)
346 ELEMENTAL FUNCTION rd(x,y,z)
RESULT(res)
351 REAL,
PARAMETER :: errtol = .0015
352 REAL,
PARAMETER :: third = 1./3.
353 REAL,
PARAMETER :: tin = 1.5e-25
354 REAL,
PARAMETER :: big = 4.5e+21
355 REAL,
PARAMETER :: c1 = 3./14.
356 REAL,
PARAMETER :: c2 = 1./6.
357 REAL,
PARAMETER :: c3 = 9./22.
358 REAL,
PARAMETER :: c4 = 3./26.
359 REAL,
PARAMETER :: c5 = .25*c3
360 REAL,
PARAMETER :: c6 = 1.5*c4
361 REAL :: alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,&
362 sqrty,sqrtz,summ,xt,yt,zt
366 IF((min(x,y).LT.0.).OR.&
367 (min(x+y,z).LT.tin).OR.&
368 (max(x,y,z).GT.big))
THEN
369 res = sqrt(-1.*abs(x))
380 alamb = sqrtx * (sqrty+sqrtz) + sqrty*sqrtz
381 summ = summ + fac/(sqrtz*(zt+alamb))
386 ave = .2*(xt+yt+3.*zt)
390 IF(max(abs(delx),abs(dely),abs(delz)).LE.errtol) &
398 res = 3. * summ + fac*(1.+ed*(-c1+c5*ed-c6*delz*ee) &
399 + delz*(c2*ee+delz*(-c3*ec+delz*c4*ea)))/(ave*sqrt(ave))
417 res = s*
rf(cos(phi)**2,(1.-s*ak)*(1.+s*ak),1.)
436 q = (1.-s*ak)*(1.+s*ak)
437 res = s*(
rf(cc,q,1.)-((s*ak)**2)*
rd(cc,q,1.)/3.)
446 REAL,
DIMENSION(0:l) :: p
471 p(i) = (1.0+p(i))*x*p(i-1) - p(i)*p(i-2)
482 REAL :: x,plminus1,plminus2,pl
486 INTENT(IN) :: l,x,plminus1,plminus2
504 pl = (1.0+c)*x*plminus1 - c*plminus2
517 REAL :: plminus1,plminus2
553 IF (abs(an-bn).GT.
eps_agm*abs(an))
THEN
557#if !defined(NECSXAURORA)
566 q = sqrt(-1.0*abs(x))
584 REAL,
DIMENSION(7),
PARAMETER &
585 :: p = (/ 1.0d0, 3.5156229d0, 3.0899424d0, &
586 1.2067492d0, 0.2659732d0, 0.360768d-1, &
588 REAL,
DIMENSION(9),
PARAMETER &
589 :: q = (/ 0.39894228d0, 0.1328592d-1, 0.225319d-2, &
590 -0.157565d-2, 0.916281d-2, -0.2057706d-1, &
591 0.2635537d-1, -0.1647633d-1, 0.392377d-2 /)
595 IF(abs(x).LT.3.75)
THEN
597 i0 = p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7))))))
601 i0 = (exp(absx)/sqrt(absx)) &
602 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4) &
603 + t*(q(5)+t*(q(6)+t*(q(7)+t*(q(8)+t*q(9)))))))))
622 REAL,
DIMENSION(7),
PARAMETER &
623 :: p = (/ 0.5d0, 0.87890594d0, 0.51498869d0, &
624 0.15084934d0, 0.2658733d-1, 0.301532d-2, &
626 REAL,
DIMENSION(9),
PARAMETER &
627 :: q = (/ 0.39894228d0, -0.3988024d-1, -0.362018d-2, &
628 0.163801d-2, -0.1031555d-1, 0.2282967d-1, &
629 -0.2895312d-1, 0.1787654d-1, -0.420059d-2 /)
633 IF(abs(x).LT.3.75)
THEN
635 i1 = x*(p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
639 i1 = (exp(absx)/sqrt(absx)) &
640 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4) &
641 + t*(q(5)+t*(q(6)+t*(q(7)+t*(q(8)+t*q(9)))))))))
664 REAL,
DIMENSION(7),
PARAMETER &
665 :: p = (/ -0.57721566d0, 0.42278420d0, 0.23069756d0, &
666 0.3488590d-1, 0.262698d-2, 0.10750d-3, &
668 REAL,
DIMENSION(7),
PARAMETER &
669 :: q = (/ 1.25331414d0, -0.7832358d-1, 0.2189568d-1, &
670 -0.1062446d-1, 0.587872d-2, -0.251540d-2, &
678 + (p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
681 k0 = (exp(-x)/sqrt(x)) &
682 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4)+t*(q(5)+t*(q(6)+t*q(7)))))))
699 REAL,
DIMENSION(7),
PARAMETER &
700 :: q = (/ 1.25331414, -0.07832358, 0.02189568,&
701 -0.01062446, 0.00587872, -0.00251540,&
710 k0e = (1.0/sqrt(x)) &
711 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4)+t*(q(5)+t*(q(6)+t*q(7)))))))
728 REAL,
DIMENSION(7),
PARAMETER &
729 :: p = (/ 1.0d0, 0.15443144d0, -0.67278579d0, &
730 -0.18156897d0, -0.1919402d-1, -0.110404d-2, &
732 REAL,
DIMENSION(7),
PARAMETER &
733 :: q = (/ 1.25331414d0, 0.23498619d0, -0.3655620d-1, &
734 0.1504268d-1, -0.780353d-2, 0.325614d-2, &
743 * (p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
746 k1 = (exp(-x)/sqrt(x)) &
747 * (q(1)+t*(q(2)+t*(q(3)+t*(q(4)+t*(q(5)+t*(q(6)+t*q(7)))))))
757 ELEMENTAL FUNCTION ei(x)
RESULT(res)
763 REAL,
PARAMETER :: eps=epsilon(x), euler=.577215664902, fpmin=tiny(x)/eps
767 ELSE IF(x.EQ.0.)
THEN
769 ELSE IF(x.LT.6.8)
THEN
771 ELSE IF(x.LT.50.)
THEN
785 REAL :: am1,a0,bm1,b0,a,b,ap1,bp1,eps
794 ap1 = b* a0 + a * am1
795 bp1 = b* b0 + a * bm1
796 eps = 10. * epsilon(x)
799 DO WHILE(abs(ap1 * b0 - a0 * bp1).GT.eps*abs(a0*bp1))
800 IF(abs(bp1).GT.1.)
THEN
813 ap1 = b * a0 + a * am1
814 bp1 = b * b0 + a * bm1
827 REAL :: xn,sn,sm1,hsum,g,y,factorial,eps
833 g = 0.5772156649015328606065121
839 DO WHILE(abs(sn-sm1).GT.eps*abs(sm1))
843 factorial = factorial * y
845 sn = sn + hsum * xn / factorial
847 res = g + log(abs(x)) - exp(x) * sn
856 REAL,
DIMENSION(44),
PARAMETER :: eic = &
857 (/1.915047433355013959531e2, 4.403798995348382689974e2, &
858 1.037878290717089587658e3, 2.492228976241877759138e3, &
859 6.071406374098611507965e3, 1.495953266639752885229e4, &
860 3.719768849068903560439e4, 9.319251363396537129882e4, &
861 2.349558524907683035782e5, 5.955609986708370018502e5, &
862 1.516637894042516884433e6, 3.877904330597443502996e6, &
863 9.950907251046844760026e6, 2.561565266405658882048e7, &
864 6.612718635548492136250e7, 1.711446713003636684975e8, &
865 4.439663698302712208698e8, 1.154115391849182948287e9, &
866 3.005950906525548689841e9, 7.842940991898186370453e9, &
867 2.049649711988081236484e10, 5.364511859231469415605e10,&
868 1.405991957584069047340e11, 3.689732094072741970640e11,&
869 9.694555759683939661662e11, 2.550043566357786926147e12,&
870 6.714640184076497558707e12, 1.769803724411626854310e13,&
871 4.669055014466159544500e13, 1.232852079912097685431e14,&
872 3.257988998672263996790e14, 8.616388199965786544948e14,&
873 2.280446200301902595341e15, 6.039718263611241578359e15,&
874 1.600664914324504111070e16, 4.244796092136850759368e16,&
875 1.126348290166966760275e17, 2.990444718632336675058e17,&
876 7.943916035704453771510e17, 2.111342388647824195000e18,&
877 5.614329680810343111535e18, 1.493630213112993142255e19,&
878 3.975442747903744836007e19, 1.058563689713169096306e20 /)
880 REAL :: xx,dx,xxj,edx,sm,sn,term,factorial,dxj,eps
889 sn = (edx - 1.) / xxj
895 DO WHILE(abs(term).GT.eps*abs(sn))
897 factorial = factorial * j
900 sm = sm + dxj / factorial
901 term = (factorial * (edx * sm - 1.)) / xxj
904 res = eic(k-6) + sn * exp(xx)
921 REAL,
INTENT(IN) :: x
924 REAL,
PARAMETER :: xerf = 0.46875
925 REAL,
PARAMETER :: xerfc = 4.0
928 IF (abs(x) .LT. tiny(1.) )
THEN
930 ELSE IF(x.GT.0.0)
THEN
931 IF (x .LT. xerf)
THEN
933 ELSE IF (x .LT. xerfc)
THEN
940 IF (-x .LT. xerf)
THEN
942 ELSE IF (-x .LT. xerfc)
THEN
955 REAL,
INTENT(IN) :: x
958 REAL,
PARAMETER :: xerf = 0.46875
959 REAL,
PARAMETER :: xerfc = 4.0
965 ELSE IF (x .LT. xerfc)
THEN
972 IF (-x .LT. xerf)
THEN
974 ELSE IF (-x .LT. xerfc)
THEN
987 REAL,
INTENT(IN) :: x
990 REAL :: denom, nom, x2, x4, x6
991 REAL,
DIMENSION(4) :: p_denom
992 REAL,
DIMENSION(5) :: p_nom
994 p_nom = (/ 3.20937758913846947e03 , 3.77485237685302021e02 ,&
995 1.13864154151050156e02 , 3.16112374387056560 ,&
996 1.85777706184603153e-1 /)
997 p_denom = (/ 2.84423683343917062e03 , 1.28261652607737228e03 ,&
998 2.44024637934444173e02 , 2.36012909523441209e01 /)
1004 nom = p_nom(1) + x2 * (p_nom(2) + p_nom(3)* x2 + p_nom(4) * x4 &
1006 denom = p_denom(1)+ x2 * (p_denom(2) + p_denom(3)* x2 + p_denom(4) * x4 &
1009 erf = x * nom / denom
1017 REAL,
INTENT(IN) :: x
1022 REAL,
DIMENSION(9) :: p_nom , p_denom
1024 p_nom = (/ 1.23033935479799725e03 , 2.05107837782607147e03 ,&
1025 1.71204761263407058e03 , 8.81952221241769090e02 ,&
1026 2.98635138197400131e02 , 6.61191906371416295e01 ,&
1027 8.88314979438837594e0 , 5.64188496988670089e-1 ,&
1028 2.15311535474403846e-8 /)
1029 p_denom = (/ 1.23033935480374942e03 , 3.43936767414372164e03 ,&
1030 4.36261909014324716e03 , 3.29079923573345963e03 ,&
1031 1.62138957456669019e03 , 5.37181101862009858e02 ,&
1032 1.17693950891312499e02 , 1.57449261107098347e01 ,&
1039 nom = (nom + p_nom(9-i)) * x
1040 denom = (denom + p_denom(9-i)) * x
1043 nom = nom + p_nom(1)
1044 denom = denom + p_denom(1)
1046 cerf = exp(-x*x)* nom / denom
1054 REAL,
INTENT(IN) :: x
1057 REAL :: denom, nom, x2, x4, x6, x8
1058 REAL,
DIMENSION(5) :: p_denom
1059 REAL,
DIMENSION(6) :: p_nom
1061 p_nom = (/ 6.58749161529837803e-4 , 1.60837851487422766e-2 ,&
1062 1.25781726111229246e-1 , 3.60344899949804439e-1 ,&
1063 3.05326634961232344e-1 , 1.63153871373020978e-2 /)
1064 p_denom = (/ 2.33520497626869185e-3 , 6.05183413124413191e-2 ,&
1065 5.27905102951428412e-1 , 1.87295284992346047e00 ,&
1066 2.56852019228982242e00 /)
1073 nom = -x2 * p_nom(1) - p_nom(2) - p_nom(3)/x2 - p_nom(4)/(x4) - &
1074 p_nom(5)/(x6) - p_nom(6)/(x8)
1075 denom = x2 * p_denom(1) + p_denom(2) + p_denom(3)/x2 + p_denom(4)/(x4) + &
1076 p_denom(5)/(x6) + 1./(x8)
1078 cerf = exp(-x*x)/x * (1./sqrt(
pi) +1./x2 * nom/denom)
1093 REAL,
INTENT(IN) :: x
1097 REAL,
PARAMETER :: error = 1e-14
1098 REAL :: q,y ,sqrtpi,sqrtpi_0_5,fy
1103 sqrtpi_0_5 = sqrtpi*0.5
1106 y = sqrtpi*(0.5*x +
pi/24.*x**3 + 7.*
pi**2/960.*x**5)
1107 q = sqrtpi_0_5 * exp(y*y)*(
errorfunc(y)-x)
1110 DO WHILE (abs(fy) .GT. error)
1113 q = sqrtpi_0_5 * exp(y*y)*(
errorfunc(y)-x)
1135 REAL,
INTENT(IN) :: x
1138 REAL,
DIMENSION(26),
PARAMETER :: g = &
1139 (/ 1.0e0, 0.5772156649015329e0,&
1140 -0.6558780715202538e0, -0.420026350340952e-1, &
1141 0.1665386113822915e0, -0.421977345555443e-1, &
1142 -0.96219715278770e-2, 0.72189432466630e-2, &
1143 -0.11651675918591e-2, -0.2152416741149e-3, &
1144 0.1280502823882e-3, -0.201348547807e-4, &
1145 -0.12504934821e-5, 0.11330272320e-5, &
1146 -0.2056338417e-6, 0.61160950e-8, &
1147 0.50020075e-8, -0.11812746e-8, &
1148 0.1043427e-9, 0.77823e-11, &
1149 -0.36968e-11, 0.51e-12, &
1150 -0.206e-13, -0.54e-14, &
1156 IF (x.EQ.int(x))
THEN
1167 IF (abs(x).GT.1.)
THEN
1183 IF (abs(x).GT.1.)
THEN
1186 res=-
pi/(x*res*sin(
pi*x))
1199 DOUBLE PRECISION :: x
1201 DOUBLE PRECISION :: alpha,beta,ag,con
1206 DOUBLE PRECISION,
DIMENSION(15),
PARAMETER &
1208 (/ .99999999999999709182d0, 57.156235665862923517d0,&
1209 -59.597960355475491248d0, 14.136097974741747174d0,&
1210 -.49191381609762019978d0, .33994649984811888699d-4,&
1211 .46523628927048575665d-4, -.98374475304879564677d-4,&
1212 .15808870322491248884d-3, -.21026444172410488319d-3,&
1213 .21743961811521264320d-3, -.16431810653676389022d-3,&
1214 .84418223983852743293d-4, -.26190838401581408670d-4,&
1215 .36899182659531622704d-5 /)
1218 alpha = 2.506628274631000502415d0
1219 beta = x + 5.24218750d0
1220 con = log(alpha)+log(beta)*(x+0.5d0) - beta
1227 ag = ag + p(i)/(x+i-1.d0)
elemental subroutine, public ellipticintegrals(k, Kell, Eell)
compute the complete elliptic integrals of first and second kind using the AGM method (arithmetic-geo...
elemental real function argument_addition_series_ei(x)
elemental real function, public lngamma(x)
computes the logarithm of the gammafunction with Lanczos approximation found in Numerical Recipes for...
elemental real function, public acosh(x)
inverse hyperbolic cosine function
elemental real function codycerf2approx(x)
Cody approximation for complementary error function for |x| > 4.
elemental real function, public incompleteellipticintegral_f(phi, ak)
elemental real function, public barrier(x, a, b)
barrier function
elemental real function rf(x, y, z)
Computes Carlson's elliptic integral of the first kind R_F(x,y,z), y,y,z > 0. One can be zero....
integer, parameter max_pl_coeff
elemental real function, public errorfunc(x)
returns the error function of argument x calculation of error fuction and complementary errorfunction...
integer, parameter max_agm_iterations
elemental real function, public bessel_k0e(x)
Compute the exponential scaled modified Bessel function of the second kind e.g. K0e = EXP(x) * K0(x) ...
elemental real function rd(x, y, z)
Computes Carlson's elliptic integral of the second kind R_D(x,y,z), y,y > 0. One can be zero....
elemental real function, public heaviside(x, a)
step function
elemental real function, public comperrorfunc(x)
elemental real function, public bessel_i1(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
elemental real function, public ei(x)
Computes the exponential integral Ei(x) for x != 0 Ei(x) := int_-oo^x exp(t)/t dt Implementation simi...
elemental real function, public ellipticintegral_k(k)
compute the complete elliptic integral of the first kind using the AGM method (arithmetic-geometric-m...
elemental real function legendrepolynomial_one(l, x, Plminus1, Plminus2)
elemental real function, public inverf(x)
Inverse of the error function Argument of this function has to be in ]-1,1[, returns huge(1....
elemental real function, public atanh(x)
inverse hyperbolic tangens function
pure subroutine, public legendrepolynomials(l, x, P)
elemental real function, public bessel_k1(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
elemental real function, public bessel_k0(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
elemental subroutine, public normalrand(u1, u2, x, y)
Box-Muller alg. for gaussian distributed random variables input are uniform distributed random variab...
elemental real function codycerf1approx(x)
Cody approximation for complementary error function for 0.46875 < |x| < 4.
elemental real function, public ellipticintegral_e(k)
compute the complete elliptic integral of the second kind using the AGM method (arithmetic-geometric-...
elemental real function, public gamma(x)
Compute the Gamma function for arguments != 0,-1,-2,..
elemental real function power_series_ei(x)
elemental real function legendrepolynomial_all(l, x)
elemental real function, public incompleteellipticintegral_e(phi, ak)
Legendre elliptic integral of the second kind E(phi,k), evaluated using Carlson's function R_D and R_...
elemental real function codyerfapprox(x)
elemental real function, public asinh(x)
inverse hyperbolic sine function
real, dimension(max_pl_coeff), parameter pl_coeff
elemental real function continued_fraction_ei(x)
elemental real function, public bessel_i0(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
elemental real function, public legendrefunction_qminhalf(x)
Computation of the order 0 and -1/2 degree Legendre function of the second kind Q_{-1/2} using the AG...