functions.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: functions.f90 #
5!# #
6!# Copyright (C) 2006-2013 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9!# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
31!----------------------------------------------------------------------------!
47!----------------------------------------------------------------------------!
48MODULE functions
49 IMPLICIT NONE
50 !--------------------------------------------------------------------------!
51 PRIVATE
52 REAL, PARAMETER :: pi = 3.14159265358979323846
53 REAL, PARAMETER :: sqrt_two = 1.41421356237309504880
54 REAL, PARAMETER :: eps_agm = epsilon(eps_agm) ! precision of the AGM
55 INTEGER, PARAMETER :: max_agm_iterations = 20 ! limit iteration steps
56 ! store the first 20 coefficients used in computing the Legendre polynomials
57 INTEGER :: iii
58 INTEGER, PARAMETER :: max_pl_coeff = 20
59 REAL, DIMENSION(MAX_PL_COEFF), PARAMETER &
60 :: pl_coeff = (/ (iii/(iii+1.0), iii=1,max_pl_coeff) /)
61 !--------------------------------------------------------------------------!
62 ! exclude interface block from doxygen processing
64 INTERFACE legendrepolynomial
65 MODULE PROCEDURE legendrepolynomial_one, legendrepolynomial_all
66 END INTERFACE
68 !--------------------------------------------------------------------------!
69 PUBLIC :: &
70 asinh, acosh, atanh, &
71 heaviside, &
72 barrier, &
73 normalrand, &
74 ellipticintegrals, &
75 ellipticintegral_k, &
76 ellipticintegral_e, &
77 incompleteellipticintegral_f, &
78 incompleteellipticintegral_e, &
79 errorfunc,&
80 comperrorfunc,&
81 inverf,&
82 legendrepolynomials, &
83 legendrepolynomial, &
84 legendrefunction_qminhalf, &
85 bessel_i0, &
86 bessel_i1, &
87 bessel_k0, &
88 bessel_k0e, &
89 bessel_k1, &
90 lngamma, &
91 ei, &
92 gamma
93 !--------------------------------------------------------------------------!
94
95CONTAINS
96
98 ! not included in FORTRAN before 2008 standard
99 ELEMENTAL FUNCTION asinh(x) RESULT(fx)
100 IMPLICIT NONE
101 REAL, INTENT(IN) :: x
102 REAL :: fx
103 fx = log(x+sqrt(x*x+1.0))
104 END FUNCTION asinh
105
106
108 ! not included in FORTRAN before 2008 standard
109 ELEMENTAL FUNCTION acosh(x) RESULT(fx)
110 IMPLICIT NONE
111 REAL, INTENT(IN) :: x
112 REAL :: fx
113 fx = log(x+sqrt(x*x-1.0))
114 END FUNCTION acosh
115
116
118 ! not included in FORTRAN before 2008 standard
119 ELEMENTAL FUNCTION atanh(x) RESULT(fx)
120 IMPLICIT NONE
121 REAL, INTENT(IN) :: x
122 REAL :: fx
123 fx = 0.5*log((1.0+x)/(1.0-x))
124 END FUNCTION atanh
125
126
128 ! returns: 0 for x<a
129 ! 1/2 for x=a
130 ! 1 for x>a
131 ELEMENTAL FUNCTION heaviside(x,a) RESULT(fx)
132 IMPLICIT NONE
133 REAL, INTENT(IN) :: x,a
134 REAL :: fx
135 fx = 0.5*(sign(1.0,x-a)+1.0)
136 END FUNCTION heaviside
137
139 ! returns 0 for x<a and x>b
140 ! 1/2 for x=a and x=b
141 ! 1 for a<x<b
142 ELEMENTAL FUNCTION barrier(x,a,b) RESULT(fx)
143 IMPLICIT NONE
144 REAL, INTENT(IN) :: x,a,b
145 REAL :: fx
146 fx = heaviside(x,a)-heaviside(x,b)
147 END FUNCTION barrier
148
151 ELEMENTAL SUBROUTINE normalrand(u1,u2,x,y)
152 IMPLICIT NONE
153 !------------------------------------------------------------------------!
154 REAL, INTENT(IN) :: u1,u2
155 REAL, INTENT(OUT) :: x,y
156 !------------------------------------------------------------------------!
157 REAL :: r,t
158 !------------------------------------------------------------------------!
159 r = sqrt(-2.0*log(u1))
160 t = 2.0*pi*u2
161 x = r*cos(t)
162 y = r*sin(t)
163 END SUBROUTINE normalrand
164
165
191 ELEMENTAL SUBROUTINE ellipticintegrals(k,Kell,Eell)
192 IMPLICIT NONE
193 !------------------------------------------------------------------------!
194 REAL, INTENT(IN) :: k ! elliptic modulus
195 REAL, INTENT(OUT) :: eell,kell
196 !------------------------------------------------------------------------!
197 INTEGER :: i,n
198 REAL :: an,bn,cn,tmp
199 !------------------------------------------------------------------------!
200 tmp = 0.0
201 IF (abs(k).LT.1.0) THEN
202 i = 1
203 an = 1.0
204 bn = sqrt(1.0-k*k)
205 cn = k
206!NEC$ UNROLL(20)
207 DO n=1,max_agm_iterations ! do not loop forever
208 tmp = tmp + cn*cn*i
209 IF (abs(cn).GT.0.5*eps_agm*abs(an)) THEN
210 cn = 0.5*(an-bn)
211 bn = sqrt(an*bn)
212 an = an-cn ! = 0.5*(an+bn_old)
213 i = ishft(i,1) ! = 2*i
214#if !defined(NECSXAURORA)
215 ELSE
216 EXIT ! exit prohibits vectorization
217#endif
218 END IF
219 END DO
220 ! Kell = 0.5*PI / an better: Kell = 0.5*PI/0.5*(an+bn)
221 kell = pi / (an + bn + tiny(kell)) ! avoid division by 0
222 ELSE
223 kell = sqrt(-1.0*abs(k)) ! return NaN
224 END IF
225 eell = (1.0-0.5*tmp)*kell
226 END SUBROUTINE ellipticintegrals
227
228
237 ELEMENTAL FUNCTION ellipticintegral_k(k) RESULT(Kell)
238 IMPLICIT NONE
239 !------------------------------------------------------------------------!
240 REAL, INTENT(IN) :: k ! elliptic modulus
241 REAL :: kell
242 !------------------------------------------------------------------------!
243 REAL :: an,bn,tmp
244 INTEGER :: i
245 !------------------------------------------------------------------------!
246 IF (abs(k).LT.1.0) THEN
247 an = 1.0
248 bn = sqrt(1.0-k*k) ! = NaN if |k| > 1
249!NEC$ UNROLL(20)
250 DO i=1,max_agm_iterations ! do not loop forever
251 IF (abs(an-bn).GT.eps_agm*abs(an)) THEN
252 tmp = 0.5*(an + bn)
253 bn = sqrt(an*bn)
254 an = tmp
255#if !defined(NECSXAURORA)
256 ELSE
257 EXIT ! exit prohibits vectorization
258#endif
259 END IF
260 END DO
261 ! Kell = 0.5*PI / an -> next iteration would compute an = 0.5*(an+bn)
262 ! hence return an even better approximation:
263 kell = pi / ( an + bn + tiny(kell) ) ! avoid division by 0
264 ELSE
265 kell = sqrt(-1.0*abs(k)) ! return NaN
266 END IF
267 END FUNCTION ellipticintegral_k
268
269
274 ELEMENTAL FUNCTION ellipticintegral_e(k) RESULT(Eell)
275 IMPLICIT NONE
276 !------------------------------------------------------------------------!
277 REAL, INTENT(IN) :: k ! elliptic modulus
278 REAL :: eell
279 !------------------------------------------------------------------------!
280 REAL :: dummy
281 !------------------------------------------------------------------------!
282 CALL ellipticintegrals(k,dummy,eell)
283 END FUNCTION ellipticintegral_e
284
285
291 ELEMENTAL FUNCTION rf(x,y,z) RESULT(res)
292 IMPLICIT NONE
293 !------------------------------------------------------------------------!
294 REAL :: x,y,z,res
295 !------------------------------------------------------------------------!
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,&
305 yt,zt
306 !------------------------------------------------------------------------!
307 INTENT(IN) :: x,y,z
308 !------------------------------------------------------------------------!
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)) ! return NAN
313 ELSE
314 xt = x
315 yt = y
316 zt = z
317 DO
318 sqrtx = sqrt(xt)
319 sqrty = sqrt(yt)
320 sqrtz = sqrt(zt)
321 alamb = sqrtx * (sqrty+sqrtz) + sqrty*sqrtz
322 xt = .25*(xt+alamb)
323 yt = .25*(yt+alamb)
324 zt = .25*(zt+alamb)
325 ave = third*(xt+yt+zt)
326 delx = (ave-xt)/ave
327 dely = (ave-yt)/ave
328 delz = (ave-zt)/ave
329 IF(max(abs(delx),abs(dely),abs(delz)).LE.errtol) &
330 EXIT
331 END DO
332
333 e2 = delx*dely-delz**2
334 e3 = delx*dely*delz
335 res = (1. + (c1*e2-c2-c3*e3)*e2 + c4*e3)/sqrt(ave)
336 END IF
337 END FUNCTION rf
338
339
346 ELEMENTAL FUNCTION rd(x,y,z) RESULT(res)
347 IMPLICIT NONE
348 !------------------------------------------------------------------------!
349 REAL :: x,y,z,res
350 !------------------------------------------------------------------------!
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
363 !------------------------------------------------------------------------!
364 INTENT(IN) :: x,y,z
365 !------------------------------------------------------------------------!
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)) ! return NAN
370 ELSE
371 xt = x
372 yt = y
373 zt = z
374 summ = 0.
375 fac = 1.
376 DO
377 sqrtx = sqrt(xt)
378 sqrty = sqrt(yt)
379 sqrtz = sqrt(zt)
380 alamb = sqrtx * (sqrty+sqrtz) + sqrty*sqrtz
381 summ = summ + fac/(sqrtz*(zt+alamb))
382 fac = .25*fac
383 xt = .25*(xt+alamb)
384 yt = .25*(yt+alamb)
385 zt = .25*(zt+alamb)
386 ave = .2*(xt+yt+3.*zt)
387 delx = (ave-xt)/ave
388 dely = (ave-yt)/ave
389 delz = (ave-zt)/ave
390 IF(max(abs(delx),abs(dely),abs(delz)).LE.errtol) &
391 EXIT
392 END DO
393 ea = delx*dely
394 eb = delz*delz
395 ec = ea - eb
396 ed = ea - 6. * eb
397 ee = ed + ec + ec
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))
400 END IF
401 END FUNCTION rd
402
403 ! Legendre elliptic integral of the first kind F(phi,k), evaluated using
404 ! Carlson's function R_F. The argument ranges are 0 <= phi <= pi/2,
405 ! 0 <= k*sin(phi) <= 1
406 ! similar to Numerical recipes 2nd edition
407 ELEMENTAL FUNCTION incompleteellipticintegral_f(phi,ak) RESULT(res)
408 IMPLICIT NONE
409 !------------------------------------------------------------------------!
410 REAL :: phi,ak,res
411 !------------------------------------------------------------------------!
412 REAL :: s
413 !------------------------------------------------------------------------!
414 INTENT(IN) :: phi,ak
415 !------------------------------------------------------------------------!
416 s = sin(phi)
417 res = s*rf(cos(phi)**2,(1.-s*ak)*(1.+s*ak),1.)
418 END FUNCTION incompleteellipticintegral_f
419
420
425 ELEMENTAL FUNCTION incompleteellipticintegral_e(phi,ak) RESULT(res)
426 IMPLICIT NONE
427 !------------------------------------------------------------------------!
428 REAL :: phi,ak,res
429 !------------------------------------------------------------------------!
430 REAL :: cc,q,s
431 !------------------------------------------------------------------------!
432 INTENT(IN) :: phi,ak
433 !------------------------------------------------------------------------!
434 s = sin(phi)
435 cc = cos(phi)**2
436 q = (1.-s*ak)*(1.+s*ak)
437 res = s*(rf(cc,q,1.)-((s*ak)**2)*rd(cc,q,1.)/3.)
438 END FUNCTION incompleteellipticintegral_e
439
440
441 PURE SUBROUTINE legendrepolynomials(l,x,P)
442 IMPLICIT NONE
443 !------------------------------------------------------------------------!
444 INTEGER :: l
445 REAL :: x
446 REAL, DIMENSION(0:l) :: p
447 !------------------------------------------------------------------------!
448 INTEGER :: i
449 !------------------------------------------------------------------------!
450 INTENT(IN) :: l,x
451 INTENT(OUT) :: p
452 !------------------------------------------------------------------------!
453 IF(l.LT.0) THEN
454 p(:) = sqrt(1.0*l) ! return NaN
455 ELSE
456 p(0) = 1.0
457 IF (l.GE.1) THEN
458 p(1) = x
459 ! use the constant coefficients for the first
460 ! MAX_PL_COEFF Legendre Polynomials
461 DO i=2,min(l,max_pl_coeff)
462 p(i) = (1.0+pl_coeff(i))*x*p(i-1) - pl_coeff(i)*p(i-2)
463 END DO
464 ! from MAX_PL_COEFF+1 to l compute the coefficients i/(i+1);
465 ! this is probably slower, because of the division
466 DO i=max_pl_coeff+1,l
467 ! recurrence formula for the Legendre polynomials:
468 ! P(i) = ( (2*i+1)*x*P(i-1) - i*P(i-2) ) / (i+1)
469 ! = (1+i/(i+1))*x*P(i-1) - i/(i+1)*P(i-2)
470 p(i) = i/(i+1.0) ! temporary (is a number close to 1)
471 p(i) = (1.0+p(i))*x*p(i-1) - p(i)*p(i-2)
472 END DO
473 END IF
474 END IF
475 END SUBROUTINE legendrepolynomials
476
477
478 ELEMENTAL FUNCTION legendrepolynomial_one(l,x,Plminus1,Plminus2) RESULT (Pl)
479 IMPLICIT NONE
480 !------------------------------------------------------------------------!
481 INTEGER :: l
482 REAL :: x,plminus1,plminus2,pl
483 !------------------------------------------------------------------------!
484 REAL :: c
485 !------------------------------------------------------------------------!
486 INTENT(IN) :: l,x,plminus1,plminus2
487 !------------------------------------------------------------------------!
488 IF(l.LT.0) THEN
489 pl = sqrt(1.0*l) ! return NaN
490 ELSE
491 SELECT CASE(l)
492 CASE(0)
493 pl = 1.0
494 CASE(1)
495 pl = x
496 CASE DEFAULT
497 ! speed up things a little with predefined coefficients
498 IF (l.LE.max_pl_coeff) THEN
499 c = pl_coeff(l)
500 ELSE
501 c = l/(l+1.0) ! temporary storage
502 END IF
503 ! do recursion step
504 pl = (1.0+c)*x*plminus1 - c*plminus2
505 END SELECT
506 END IF
507 END FUNCTION legendrepolynomial_one
508
509
510 ELEMENTAL FUNCTION legendrepolynomial_all(l,x) RESULT (Pl)
511 IMPLICIT NONE
512 !------------------------------------------------------------------------!
513 INTEGER :: l
514 REAL :: x,pl
515 !------------------------------------------------------------------------!
516 INTEGER :: i
517 REAL :: plminus1,plminus2
518 !------------------------------------------------------------------------!
519 INTENT(IN) :: l,x
520 !------------------------------------------------------------------------!
521 IF(l.LT.0) THEN
522 pl = sqrt(1.0*l) ! return NaN
523 ELSE
524!NEC$ UNROLL(20)
525 DO i=0,l
526 pl = legendrepolynomial_one(i,x,plminus1,plminus2)
527 plminus2 = plminus1
528 plminus1 = pl
529 END DO
530 END IF
531 END FUNCTION legendrepolynomial_all
532
538 ELEMENTAL FUNCTION legendrefunction_qminhalf(x) RESULT (q)
539 IMPLICIT NONE
540 !------------------------------------------------------------------------!
541 REAL :: x,q
542 !------------------------------------------------------------------------!
543 INTEGER :: i
544 REAL :: an,bn,tmp
545 !------------------------------------------------------------------------!
546 INTENT(IN) :: x
547 !------------------------------------------------------------------------!
548 IF (x.GE.1.0) THEN
549 an = sqrt(x-1.0)
550 bn = sqrt(x+1.0) ! if x < 1 bn = NaN, hence q = NaN
551!NEC$ UNROLL(20)
552 DO i=1,max_agm_iterations ! do not loop forever
553 IF (abs(an-bn).GT.eps_agm*abs(an)) THEN
554 tmp = 0.5*(an + bn)
555 bn = sqrt(an*bn)
556 an = tmp
557#if !defined(NECSXAURORA)
558 ELSE
559 EXIT ! exit prohibits vectorization
560#endif
561 END IF
562 END DO
563! q = 0.5*PI*SQRT_TWO / an
564 q = pi*sqrt_two / (an+bn+tiny(q))
565 ELSE
566 q = sqrt(-1.0*abs(x)) ! return NaN
567 END IF
568 END FUNCTION legendrefunction_qminhalf
569
570
571
576 ELEMENTAL FUNCTION bessel_i0(x) RESULT(I0)
577 IMPLICIT NONE
578 !------------------------------------------------------------------------!
579 REAL :: x
580 REAL :: i0
581 !------------------------------------------------------------------------!
582 INTENT(IN) :: x
583 !------------------------------------------------------------------------!
584 REAL, DIMENSION(7), PARAMETER &
585 :: p = (/ 1.0d0, 3.5156229d0, 3.0899424d0, &
586 1.2067492d0, 0.2659732d0, 0.360768d-1, &
587 0.45813d-2 /)
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 /)
592 REAL :: t, absx
593 !------------------------------------------------------------------------!
594
595 IF(abs(x).LT.3.75) THEN
596 t = (x/3.75)**2
597 i0 = p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7))))))
598 ELSE
599 absx = abs(x)
600 t = 3.75/absx
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)))))))))
604 ENDIF
605
606 END FUNCTION bessel_i0
607
608
609
614 ELEMENTAL FUNCTION bessel_i1(x) RESULT(I1)
615 IMPLICIT NONE
616 !------------------------------------------------------------------------!
617 REAL :: x
618 REAL :: i1
619 !------------------------------------------------------------------------!
620 INTENT(IN) :: x
621 !------------------------------------------------------------------------!
622 REAL, DIMENSION(7), PARAMETER &
623 :: p = (/ 0.5d0, 0.87890594d0, 0.51498869d0, &
624 0.15084934d0, 0.2658733d-1, 0.301532d-2, &
625 0.32411d-3 /)
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 /)
630 REAL :: t, absx
631 !------------------------------------------------------------------------!
632
633 IF(abs(x).LT.3.75) THEN
634 t = (x/3.75)**2
635 i1 = x*(p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
636 ELSE
637 absx = abs(x)
638 t = 3.75/absx
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)))))))))
642 IF(x.LT.0.) THEN
643 i1 = -i1
644 END IF
645
646 ENDIF
647
648 END FUNCTION bessel_i1
649
650
651
656 ELEMENTAL FUNCTION bessel_k0(x) RESULT(K0)
657 IMPLICIT NONE
658 !------------------------------------------------------------------------!
659 REAL :: x
660 REAL :: k0
661 !------------------------------------------------------------------------!
662 INTENT(IN) :: x
663 !------------------------------------------------------------------------!
664 REAL, DIMENSION(7), PARAMETER &
665 :: p = (/ -0.57721566d0, 0.42278420d0, 0.23069756d0, &
666 0.3488590d-1, 0.262698d-2, 0.10750d-3, &
667 0.74d-5 /)
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, &
671 0.53208d-3 /)
672 REAL :: t
673 !------------------------------------------------------------------------!
674
675 IF(x.LE.2.0) THEN
676 t = x*x / 4.0
677 k0 = (-log(x/2.0)*bessel_i0(x)) &
678 + (p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
679 ELSE
680 t = (2.0/x)
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)))))))
683 ENDIF
684
685 END FUNCTION bessel_k0
686
691 ELEMENTAL FUNCTION bessel_k0e(x) RESULT(K0e)
692 IMPLICIT NONE
693 !------------------------------------------------------------------------!
694 REAL :: x
695 REAL :: k0e
696 !------------------------------------------------------------------------!
697 INTENT(IN) :: x
698 !------------------------------------------------------------------------!
699 REAL, DIMENSION(7), PARAMETER &
700 :: q = (/ 1.25331414, -0.07832358, 0.02189568,&
701 -0.01062446, 0.00587872, -0.00251540,&
702 0.00053208 /)
703 REAL :: t
704 !------------------------------------------------------------------------!
705
706 IF(x.LT.2.0) THEN
707 k0e = exp(x) * bessel_k0(x)
708 ELSE
709 t = (2.0/x)
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)))))))
712 ENDIF
713
714 END FUNCTION bessel_k0e
715
720 ELEMENTAL FUNCTION bessel_k1(x) RESULT(K1)
721 IMPLICIT NONE
722 !------------------------------------------------------------------------!
723 REAL :: x
724 REAL :: k1
725 !------------------------------------------------------------------------!
726 INTENT(IN) :: x
727 !------------------------------------------------------------------------!
728 REAL, DIMENSION(7), PARAMETER &
729 :: p = (/ 1.0d0, 0.15443144d0, -0.67278579d0, &
730 -0.18156897d0, -0.1919402d-1, -0.110404d-2, &
731 -0.4686d-4 /)
732 REAL, DIMENSION(7), PARAMETER &
733 :: q = (/ 1.25331414d0, 0.23498619d0, -0.3655620d-1, &
734 0.1504268d-1, -0.780353d-2, 0.325614d-2, &
735 -0.68245d-3 /)
736 REAL :: t
737 !------------------------------------------------------------------------!
738
739 IF(x.LE.2.0) THEN
740 t = x*x / 4.0
741 k1 = (log(x/2.0)*bessel_i1(x)) &
742 + (1.0/x) &
743 * (p(1)+t*(p(2)+t*(p(3)+t*(p(4)+t*(p(5)+t*(p(6)+t*p(7)))))))
744 ELSE
745 t = (2.0/x)
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)))))))
748 ENDIF
749
750 END FUNCTION bessel_k1
751
752
757 ELEMENTAL FUNCTION ei(x) RESULT(res)
758 IMPLICIT NONE
759 !------------------------------------------------------------------------!
760 REAL, INTENT(IN):: x
761 REAL :: res
762 !------------------------------------------------------------------------!
763 REAL, PARAMETER :: eps=epsilon(x), euler=.577215664902, fpmin=tiny(x)/eps
764 !------------------------------------------------------------------------!
765 IF(x.LT.-5.) THEN
766 res = continued_fraction_ei(x)
767 ELSE IF(x.EQ.0.) THEN
768 res = -huge(res)
769 ELSE IF(x.LT.6.8) THEN
770 res = power_series_ei(x)
771 ELSE IF(x.LT.50.) THEN
772 res = argument_addition_series_ei(x)
773 ELSE
774 res = continued_fraction_ei(x)
775 END IF
776 END FUNCTION ei
777
778
779 ELEMENTAL FUNCTION continued_fraction_ei(x) RESULT(res)
780 IMPLICIT NONE
781 !------------------------------------------------------------------------!
782 REAL, INTENT(IN):: x
783 REAL :: res
784 !------------------------------------------------------------------------!
785 REAL :: am1,a0,bm1,b0,a,b,ap1,bp1,eps
786 INTEGER :: j
787 !------------------------------------------------------------------------!
788 am1 = 1.
789 a0 = 0.
790 bm1 = 0.
791 b0 = 1.
792 a = exp(x)
793 b = -x + 1.
794 ap1 = b* a0 + a * am1
795 bp1 = b* b0 + a * bm1
796 eps = 10. * epsilon(x)
797 j = 1
798 a = 1.
799 DO WHILE(abs(ap1 * b0 - a0 * bp1).GT.eps*abs(a0*bp1))
800 IF(abs(bp1).GT.1.) THEN
801 am1 = a0 / bp1
802 a0 = ap1 / bp1
803 bm1 = b0 / bp1
804 b0 = 1.
805 ELSE
806 am1 = a0
807 a0 = ap1
808 bm1 = b0
809 b0 = bp1
810 END IF
811 a = -j*j
812 b = b + 2.
813 ap1 = b * a0 + a * am1
814 bp1 = b * b0 + a * bm1
815 j = j + 1
816 END DO
817 res = -ap1 / bp1
818 END FUNCTION continued_fraction_ei
819
820
821 ELEMENTAL FUNCTION power_series_ei(x) RESULT(res)
822 IMPLICIT NONE
823 !------------------------------------------------------------------------!
824 REAL, INTENT(IN):: x
825 REAL :: res
826 !------------------------------------------------------------------------!
827 REAL :: xn,sn,sm1,hsum,g,y,factorial,eps
828 !------------------------------------------------------------------------!
829 xn = -x
830 sn = -x
831 sm1 = 0.
832 hsum = 1.
833 g = 0.5772156649015328606065121
834 y = 1.
835 factorial = 1.
836 eps = 10.*epsilon(x)
837 IF(x.EQ.0.) &
838 res = -huge(res)
839 DO WHILE(abs(sn-sm1).GT.eps*abs(sm1))
840 sm1 = sn
841 y = y + 1.
842 xn = xn * (-x)
843 factorial = factorial * y
844 hsum = hsum + 1. / y
845 sn = sn + hsum * xn / factorial
846 END DO
847 res = g + log(abs(x)) - exp(x) * sn
848 END FUNCTION power_series_ei
849
850 ELEMENTAL FUNCTION argument_addition_series_ei(x) RESULT(res)
851 IMPLICIT NONE
852 !------------------------------------------------------------------------!
853 REAL, INTENT(IN):: x
854 REAL :: res
855 !------------------------------------------------------------------------!
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 /)
879 INTEGER :: k, j
880 REAL :: xx,dx,xxj,edx,sm,sn,term,factorial,dxj,eps
881 !------------------------------------------------------------------------!
882 k = nint(x + 0.5)
883 j = 0
884 xx = k
885 dx = x - xx
886 xxj = xx
887 edx = exp(dx)
888 sm = 1.
889 sn = (edx - 1.) / xxj
890 term = huge(x)
891 factorial = 1.
892 dxj = 1.
893 eps = 10.*epsilon(x)
894
895 DO WHILE(abs(term).GT.eps*abs(sn))
896 j = j + 1
897 factorial = factorial * j
898 xxj = xxj * xx
899 dxj = dxj * (-dx)
900 sm = sm + dxj / factorial
901 term = (factorial * (edx * sm - 1.)) / xxj
902 sn = sn + term
903 END DO
904 res = eic(k-6) + sn * exp(xx)
905 END FUNCTION argument_addition_series_ei
906
907
918 ELEMENTAL FUNCTION errorfunc(x) RESULT(erf)
919 IMPLICIT NONE
920 !------------------------------------------------------------------------!
921 REAL, INTENT(IN) :: x
922 REAL :: erf
923 !------------------------------------------------------------------------!
924 REAL, PARAMETER :: xerf = 0.46875
925 REAL, PARAMETER :: xerfc = 4.0
926 !------------------------------------------------------------------------!
927 !Cody error function approximation
928 IF (abs(x) .LT. tiny(1.) ) THEN !!0 besser abfangen
929 erf = 0.0
930 ELSE IF(x.GT.0.0)THEN
931 IF (x .LT. xerf) THEN
932 erf = codyerfapprox(x)
933 ELSE IF (x .LT. xerfc) THEN
934 erf = 1. - codycerf1approx(x)
935 ELSE
936 erf = 1. - codycerf2approx(x)
937 END IF
938 !use erf,erfc symmetry for negative arguments
939 ELSE
940 IF (-x .LT. xerf) THEN
941 erf = - codyerfapprox(-x)
942 ELSE IF (-x .LT. xerfc) THEN
943 erf = codycerf1approx(-x) - 1.
944 ELSE
945 erf = codycerf2approx(-x) - 1.
946 END IF
947 END IF
948
949 END FUNCTION errorfunc
950
951 ! returns the complementary error function
952 ELEMENTAL FUNCTION comperrorfunc(x) RESULT(cerf)
953 IMPLICIT NONE
954 !------------------------------------------------------------------------!
955 REAL, INTENT(IN) :: x
956 REAL :: cerf
957 !------------------------------------------------------------------------!
958 REAL, PARAMETER :: xerf = 0.46875
959 REAL, PARAMETER :: xerfc = 4.0
960 !------------------------------------------------------------------------!
961 !Cody complementary error function approximation
962 IF(x.GE.0.0) THEN
963 IF (x .LT.xerf) THEN
964 cerf = 1. - codyerfapprox(x)
965 ELSE IF (x .LT. xerfc) THEN
966 cerf = codycerf1approx(x)
967 ELSE
968 cerf = codycerf2approx(x)
969 END IF
970 !use erf,erfc symmetry for negative arguments
971 ELSE
972 IF (-x .LT. xerf) THEN
973 cerf = 1. + codyerfapprox(-x)
974 ELSE IF (-x .LT. xerfc) THEN
975 cerf = 2. - codycerf1approx(-x)
976 ELSE
977 cerf = 2. - codycerf2approx(-x)
978 END IF
979 END IF
980
981 END FUNCTION comperrorfunc
982
983
984 ELEMENTAL FUNCTION codyerfapprox(x) RESULT(erf)
985 IMPLICIT NONE
986 !------------------------------------------------------------------------!
987 REAL, INTENT(IN) :: x
988 REAL :: erf
989 !------------------------------------------------------------------------!
990 REAL :: denom, nom, x2, x4, x6
991 REAL,DIMENSION(4) :: p_denom
992 REAL,DIMENSION(5) :: p_nom
993 !------------------------------------------------------------------------!
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 /)
999
1000 x2 = x * x
1001 x4 = x2 * x2
1002 x6 = x4 * x2
1003
1004 nom = p_nom(1) + x2 * (p_nom(2) + p_nom(3)* x2 + p_nom(4) * x4 &
1005 + p_nom(5)*x6)
1006 denom = p_denom(1)+ x2 * (p_denom(2) + p_denom(3)* x2 + p_denom(4) * x4 &
1007 + x6)
1008
1009 erf = x * nom / denom
1010
1011 END FUNCTION codyerfapprox
1012
1014 ELEMENTAL FUNCTION codycerf1approx(x) RESULT(cerf)
1015 IMPLICIT NONE
1016 !------------------------------------------------------------------------!
1017 REAL, INTENT(IN) :: x
1018 REAL :: cerf
1019 !------------------------------------------------------------------------!
1020 INTEGER :: i
1021 REAL :: denom, nom
1022 REAL,DIMENSION(9) :: p_nom , p_denom
1023 !------------------------------------------------------------------------!
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 ,&
1033 1.0 /)
1034
1035 nom = 0.0
1036 denom = 0.0
1037
1038 DO i= 0,7
1039 nom = (nom + p_nom(9-i)) * x
1040 denom = (denom + p_denom(9-i)) * x
1041 END DO
1042
1043 nom = nom + p_nom(1)
1044 denom = denom + p_denom(1)
1045
1046 cerf = exp(-x*x)* nom / denom
1047
1048 END FUNCTION codycerf1approx
1049
1051 ELEMENTAL FUNCTION codycerf2approx(x) RESULT(cerf)
1052 IMPLICIT NONE
1053 !------------------------------------------------------------------------!
1054 REAL, INTENT(IN) :: x
1055 REAL :: cerf
1056 !------------------------------------------------------------------------!
1057 REAL :: denom, nom, x2, x4, x6, x8
1058 REAL,DIMENSION(5) :: p_denom
1059 REAL,DIMENSION(6) :: p_nom
1060 !------------------------------------------------------------------------!
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 /)
1067
1068 x2 = x * x
1069 x4 = x2 * x2
1070 x6 = x4 * x2
1071 x8 = x4 * x4
1072
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)
1077
1078 cerf = exp(-x*x)/x * (1./sqrt(pi) +1./x2 * nom/denom)
1079
1080 END FUNCTION codycerf2approx
1081
1090 ELEMENTAL FUNCTION inverf(x) RESULT(ierf)
1091 IMPLICIT NONE
1092 !------------------------------------------------------------------------!
1093 REAL, INTENT(IN) :: x
1094 REAL :: ierf
1095 !------------------------------------------------------------------------!
1096! FIXME : find global error variable, depending on precision(single,double,..)
1097 REAL, PARAMETER :: error = 1e-14
1098 REAL :: q,y ,sqrtpi,sqrtpi_0_5,fy
1099 !------------------------------------------------------------------------!
1100
1101 IF(x.LT.1.0)THEN
1102 sqrtpi = sqrt(pi)
1103 sqrtpi_0_5 = sqrtpi*0.5
1104
1105 ! initial guess:
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) !f(y)/f'(y)
1108 fy = errorfunc(y) - x
1109
1110 DO WHILE (abs(fy) .GT. error)
1111 y = y - q/(1.+q*y)
1112 fy = errorfunc(y) - x
1113 q = sqrtpi_0_5 * exp(y*y)*(errorfunc(y)-x)
1114 END DO
1115 ierf = y
1116 ELSE
1117 ! infinity
1118 ierf = huge(1.)
1119 END IF
1120
1121 END FUNCTION inverf
1122
1132 ELEMENTAL FUNCTION gamma(x) RESULT(res)
1133 IMPLICIT NONE
1134 !------------------------------------------------------------------------!
1135 REAL, INTENT(IN) :: x
1136 REAL :: res
1137 !------------------------------------------------------------------------!
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, &
1151 0.14e-14, 0.1e-15/)
1152 INTEGER :: m1,k,m
1153 REAL :: z,r,gr
1154 !------------------------------------------------------------------------!
1155 r=1 ! this surpresses the 'maybe-uninitialized' compiler warning
1156 IF (x.EQ.int(x)) THEN
1157 IF (x.GT.0.) THEN
1158 res=1.0
1159 m1=x-1
1160 DO k=2,m1
1161 res=res*k
1162 END DO
1163 ELSE
1164 res=1.0e+300
1165 END IF
1166 ELSE
1167 IF (abs(x).GT.1.) THEN
1168 z=abs(x)
1169 m=int(z)
1170 !R=1.
1171 DO k=1,m
1172 r=r*(z-k)
1173 END DO
1174 z=z-m
1175 ELSE
1176 z=x
1177 ENDIF
1178 gr=g(26)
1179 DO k=25,1,-1
1180 gr=gr*z+g(k)
1181 END DO
1182 res=1.0/(gr*z)
1183 IF (abs(x).GT.1.) THEN
1184 res=res*r
1185 IF (x.LT.0.) &
1186 res=-pi/(x*res*sin(pi*x))
1187 ENDIF
1188 ENDIF
1189 END FUNCTION gamma
1190
1196 ELEMENTAL FUNCTION lngamma(x) RESULT(lg)
1197 IMPLICIT NONE
1198 !------------------------------------------------------------------------!
1199 DOUBLE PRECISION :: x
1200 REAL :: lg
1201 DOUBLE PRECISION :: alpha,beta,ag,con
1202 INTEGER :: i
1203 !------------------------------------------------------------------------!
1204 INTENT(IN) :: x
1205 !------------------------------------------------------------------------!
1206 DOUBLE PRECISION, DIMENSION(15), PARAMETER &
1207 :: p = &
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 /)
1216
1217 ! prefactor
1218 alpha = 2.506628274631000502415d0
1219 beta = x + 5.24218750d0
1220 con = log(alpha)+log(beta)*(x+0.5d0) - beta
1221 lg = 0
1222
1223 ! sum
1224 IF (x.GE.0.0) THEN
1225 ag = p(1)
1226 DO i=2,15
1227 ag = ag + p(i)/(x+i-1.d0)
1228 END DO
1229 lg = con + log(ag)
1230 ELSE
1231! CALL Error(this,"LnGamma","bad value for for gammafunction")
1232 END IF
1233 END FUNCTION lngamma
1234END MODULE functions
functions
Mathematical functions.
Definition: functions.f90:48
functions::ellipticintegrals
elemental subroutine, public ellipticintegrals(k, Kell, Eell)
compute the complete elliptic integrals of first and second kind using the AGM method (arithmetic-geo...
Definition: functions.f90:192
functions::argument_addition_series_ei
elemental real function argument_addition_series_ei(x)
Definition: functions.f90:851
functions::lngamma
elemental real function, public lngamma(x)
computes the logarithm of the gammafunction with Lanczos approximation found in Numerical Recipes for...
Definition: functions.f90:1197
functions::acosh
elemental real function, public acosh(x)
inverse hyperbolic cosine function
Definition: functions.f90:110
functions::codycerf2approx
elemental real function codycerf2approx(x)
Cody approximation for complementary error function for |x| > 4.
Definition: functions.f90:1052
functions::incompleteellipticintegral_f
elemental real function, public incompleteellipticintegral_f(phi, ak)
Definition: functions.f90:408
functions::barrier
elemental real function, public barrier(x, a, b)
barrier function
Definition: functions.f90:143
functions::iii
integer iii
Definition: functions.f90:57
functions::rf
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....
Definition: functions.f90:292
functions::max_pl_coeff
integer, parameter max_pl_coeff
Definition: functions.f90:58
functions::errorfunc
elemental real function, public errorfunc(x)
returns the error function of argument x calculation of error fuction and complementary errorfunction...
Definition: functions.f90:919
functions::max_agm_iterations
integer, parameter max_agm_iterations
Definition: functions.f90:55
functions::bessel_k0e
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) ...
Definition: functions.f90:692
functions::rd
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....
Definition: functions.f90:347
functions::heaviside
elemental real function, public heaviside(x, a)
step function
Definition: functions.f90:132
functions::comperrorfunc
elemental real function, public comperrorfunc(x)
Definition: functions.f90:953
functions::bessel_i1
elemental real function, public bessel_i1(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
Definition: functions.f90:615
functions::ei
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...
Definition: functions.f90:758
functions::ellipticintegral_k
elemental real function, public ellipticintegral_k(k)
compute the complete elliptic integral of the first kind using the AGM method (arithmetic-geometric-m...
Definition: functions.f90:238
functions::legendrepolynomial_one
elemental real function legendrepolynomial_one(l, x, Plminus1, Plminus2)
Definition: functions.f90:479
functions::inverf
elemental real function, public inverf(x)
Inverse of the error function Argument of this function has to be in ]-1,1[, returns huge(1....
Definition: functions.f90:1091
functions::atanh
elemental real function, public atanh(x)
inverse hyperbolic tangens function
Definition: functions.f90:120
functions::legendrepolynomials
pure subroutine, public legendrepolynomials(l, x, P)
Definition: functions.f90:442
functions::bessel_k1
elemental real function, public bessel_k1(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
Definition: functions.f90:721
functions::bessel_k0
elemental real function, public bessel_k0(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
Definition: functions.f90:657
functions::normalrand
elemental subroutine, public normalrand(u1, u2, x, y)
Box-Muller alg. for gaussian distributed random variables input are uniform distributed random variab...
Definition: functions.f90:152
functions::codycerf1approx
elemental real function codycerf1approx(x)
Cody approximation for complementary error function for 0.46875 < |x| < 4.
Definition: functions.f90:1015
functions::ellipticintegral_e
elemental real function, public ellipticintegral_e(k)
compute the complete elliptic integral of the second kind using the AGM method (arithmetic-geometric-...
Definition: functions.f90:275
functions::gamma
elemental real function, public gamma(x)
Compute the Gamma function for arguments != 0,-1,-2,..
Definition: functions.f90:1133
functions::pi
real, parameter pi
Definition: functions.f90:52
functions::sqrt_two
real, parameter sqrt_two
Definition: functions.f90:53
functions::power_series_ei
elemental real function power_series_ei(x)
Definition: functions.f90:822
functions::legendrepolynomial_all
elemental real function legendrepolynomial_all(l, x)
Definition: functions.f90:511
functions::incompleteellipticintegral_e
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_...
Definition: functions.f90:426
functions::codyerfapprox
elemental real function codyerfapprox(x)
Definition: functions.f90:985
functions::asinh
elemental real function, public asinh(x)
inverse hyperbolic sine function
Definition: functions.f90:100
functions::pl_coeff
real, dimension(max_pl_coeff), parameter pl_coeff
Definition: functions.f90:59
functions::continued_fraction_ei
elemental real function continued_fraction_ei(x)
Definition: functions.f90:780
functions::bessel_i0
elemental real function, public bessel_i0(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
Definition: functions.f90:577
functions::legendrefunction_qminhalf
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...
Definition: functions.f90:539
functions::eps_agm
real, parameter eps_agm
Definition: functions.f90:54
Generated by   doxygen 1.9.4