# include "define.inc" ! ! special function wrapper routine written by Tomo Tatsuno (5/7/08) ! ! RN 2018/08/01: Added _SPF200X_ for SPFUNC to use Fortran 2003/2008 standard functions, ! bessel_j0, log_gamma etc. ! RN 2018/07/31: Added partial support for quad precision functions ! define QUAD_PRECISION to use quad precision functions if available ! RN ??: Added Log Gamma function ! RN 2008/07/01: Error function is added ! RN 2008/07/01: Compilers not having intrinsic those special functions ! must choose one of the following ! 1: local [USE_LOCAL_SPFUNC=on] ! 2: NAG Library [USE_NAGLIB=spfunc] ! TT 2008/05/07: Only has Bessel function at the moment ! ! Unfortunately we do not support elemental feature for ifort... ! XL-fortran does not seem to have intrinsic Bessel function ! ! To do: avoid explicit specification of kind and use kind_rs or kind_rd ! support of other compilers such as absoft, lahay etc... ! ! PGI cannot use kind_rs, kind_rd as kind specifications, but it's ok ! because we know kind_rs=4 kind_rd=8 for PGI ! !> FIXME : Add documentation (or tidy existing documentation) module spfunc use constants, only: kind_rs, kind_rd # ifdef QUAD_PRECISION use constants, only: kind_rq # endif # if SPFUNC == _SPNAG_ use constants, only: nag_kind # endif implicit none private public :: j0, j1 public :: erf_ext public :: lgamma_ext # if SPFUNC == _SPNAG_ ! error handling ! ifail=0: terminate the program ! ifail=1 (-1): continues the program w/o (w) error messages integer :: ifail=-1 # endif ! define interface # if SPFUNC == _SPLOCAL_ # elif SPFUNC == _SPNAG_ # elif SPFUNC == _SPF200X_ interface j0 module procedure sj0, dj0 # ifdef QUAD_PRECISION module procedure qj0 # endif end interface j0 interface j1 module procedure sj1, dj1 # ifdef QUAD_PRECISION module procedure qj1 # endif end interface j1 interface erf_ext module procedure serf_ext, derf_ext # ifdef QUAD_PRECISION module procedure qerf_ext # endif end interface erf_ext interface lgamma_ext module procedure slgamma_ext, dlgamma_ext # ifdef QUAD_PRECISION module procedure qlgamma_ext # endif end interface lgamma_ext # else /* if not _SPLOCAL_ and not _SPNAG_ and not _F200X_ */ # if FCOMPILER == _PATHSCALE_ interface j0 module procedure sj0, dj0 end interface interface j1 module procedure sj1, dj1 end interface interface erf_ext module procedure serf_ext, derf_ext end interface # elif ( FCOMPILER == _INTEL_ || FCOMPILER == _GFORTRAN_ || FCOMPILER == _IFX_ ) /* not _PATHSCALE_ */ interface j0 module procedure sj0, dj0 # ifdef QUAD_PRECISION # if FCOMPILER == _GFORTRAN_ module procedure qj0 # endif # endif end interface j0 interface j1 module procedure sj1, dj1 # ifdef QUAD_PRECISION # if FCOMPILER == _GFORTRAN_ module procedure qj1 # endif # endif end interface j1 interface erf_ext module procedure serf_ext, derf_ext # ifdef QUAD_PRECISION # if FCOMPILER == _GFORTRAN_ module procedure qerf_ext # endif # endif end interface erf_ext # elif FCOMPILER == _G95_ /* not _GFORTRAN_, _INTEL_, _PATHSCALE_ */ interface j0 module procedure sj0, dj0 end interface interface j1 module procedure sj1, dj1 end interface interface erf_ext module procedure serf_ext, derf_ext end interface interface lgamma_ext module procedure slgamma_ext, dlgamma_ext end interface # elif FCOMPILER == _PGI_ /* not _G95_, _GFORTRAN_, _INTEL_, _PATHSCALE_*/ interface j0 elemental function besj0(x) real (kind=4), intent (in) :: x real (kind=4) :: besj0 end function besj0 elemental function dbesj0(x) real (kind=8), intent (in) :: x real (kind=8) :: dbesj0 end function dbesj0 end interface j0 interface j1 module procedure sj1, dj1 end interface j1 interface erf_ext elemental function erf(x) real (kind=4), intent (in) :: x real (kind=4) :: erf end function erf elemental function derf(x) real (kind=8), intent (in) :: x real (kind=8) :: derf end function derf end interface erf_ext # elif FCOMPILER == _NEC_ /* not _FUJ_, _PGI_, _G95_, _GFORTRAN_, _INTEL_, _PATHSCALE */ interface j0 module procedure j0s, j0a end interface interface j1 module procedure j1s, j1a end interface # endif /* FCOMPILER */ # endif /* SPFUNC */ contains # if SPFUNC == _SPLOCAL_ !-------------------------------------------------------------------------! ! double precision Bessel functions (dbesj0.f dbesj1.f) ! ! http://www.kurims.kyoto-u.ac.jp/~ooura/bessel.html ! !-------------------------------------------------------------------------! ! Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp). ! ! You may use, copy, modify this code for any purpose and ! ! without fee. You may distribute this ORIGINAL package. ! !-------------------------------------------------------------------------! ! Modified by Ryusuke NUMATA 2008/06/27 ! ! to fit F90 format, and for j1 to give besj1/x ! !-------------------------------------------------------------------------! ! these routines are declared as real, but can be promoted to ! double presicion using compiler option. ! this is so because constatns used are only double precision, ! and cannot be promoted to quad-precision. ! differences from GNU's gfortran besj0 and besj1 are ! |err| < 1.e-15 !> FIXME : Add documentation elemental function j0(x) ! Bessel J_0(x) function in double precision use constants, only: pi implicit none real :: j0 real, intent(in) :: x integer :: k real :: w, t, y, theta, v real, parameter :: a(0:7) = (/ & & -0.0000000000023655394, 0.0000000004708898680, & & -0.0000000678167892231, 0.0000067816840038636, & & -0.0004340277777716935, 0.0156249999999992397, & & -0.2499999999999999638, 0.9999999999999999997 /) real, parameter :: b(0:64) = (/ & & 0.0000000000626681117, -0.0000000022270614428, & & 0.0000000662981656302, -0.0000016268486502196, & & 0.0000321978384111685, -0.0005005237733315830, & & 0.0059060313537449816, -0.0505265323740109701, & & 0.2936432097610503985, -1.0482565081091638637, & & 1.9181123286040428113, -1.1319199475221700100, & & -0.1965480952704682000, 0.0000000000457457332, & & -0.0000000015814772025, 0.0000000455487446311, & & -0.0000010735201286233, 0.0000202015179970014, & & -0.0002942392368203808, 0.0031801987726150648, & & -0.0239875209742846362, 0.1141447698973777641, & & -0.2766726722823530233, 0.1088620480970941648, & & 0.5136514645381999197, -0.2100594022073706033, & & 0.0000000000331366618, -0.0000000011119090229, & & 0.0000000308823040363, -0.0000006956602653104, & & 0.0000123499947481762, -0.0001662951945396180, & & 0.0016048663165678412, -0.0100785479932760966, & & 0.0328996815223415274, -0.0056168761733860688, & & -0.2341096400274429386, 0.2551729256776404262, & & 0.2288438186148935667, 0.0000000000238007203, & & -0.0000000007731046439, 0.0000000206237001152, & & -0.0000004412291442285, 0.0000073107766249655, & & -0.0000891749801028666, 0.0007341654513841350, & & -0.0033303085445352071, 0.0015425853045205717, & & 0.0521100583113136379, -0.1334447768979217815, & & -0.1401330292364750968, 0.2685616168804818919, & & 0.0000000000169355950, -0.0000000005308092192, & & 0.0000000135323005576, -0.0000002726650587978, & & 0.0000041513240141760, -0.0000443353052220157, & & 0.0002815740758993879, -0.0004393235121629007, & & -0.0067573531105799347, 0.0369141914660130814, & & 0.0081673361942996237, -0.2573381285898881860, & & 0.0459580257102978932 /) real, parameter :: c(0:69) = (/ & & -0.00000000003009451757, -0.00000000014958003844, & & 0.00000000506854544776, 0.00000001863564222012, & & -0.00000060304249068078, -0.00000147686259937403, & & 0.00004714331342682714, 0.00006286305481740818, & & -0.00214137170594124344, -0.00089157336676889788, & & 0.04508258728666024989, -0.00490362805828762224, & & -0.27312196367405374426, 0.04193925184293450356, & & -0.00000000000712453560, -0.00000000041170814825, & & 0.00000000138012624364, 0.00000005704447670683, & & -0.00000019026363528842, -0.00000533925032409729, & & 0.00001736064885538091, 0.00030692619152608375, & & -0.00092598938200644367, -0.00917934265960017663, & & 0.02287952522866389076, 0.10545197546252853195, & & -0.16126443075752985095, -0.19392874768742235538, & & 0.00000000002128344556, -0.00000000031053910272, & & -0.00000000334979293158, 0.00000004507232895050, & & 0.00000036437959146427, -0.00000446421436266678, & & -0.00002523429344576552, 0.00027519882931758163, & & 0.00097185076358599358, -0.00898326746345390692, & & -0.01665959196063987584, 0.11456933464891967814, & & 0.07885001422733148815, -0.23664819446234712621, & & 0.00000000003035295055, 0.00000000005486066835, & & -0.00000000501026824811, -0.00000000501246847860, & & 0.00000058012340163034, 0.00000016788922416169, & & -0.00004373270270147275, 0.00001183898532719802, & & 0.00189863342862291449, -0.00113759249561636130, & & -0.03846797195329871681, 0.02389746880951420335, & & 0.22837862066532347461, -0.06765394811166522844, & & 0.00000000001279875977, 0.00000000035925958103, & & -0.00000000228037105967, -0.00000004852770517176, & & 0.00000028696428000189, 0.00000440131125178642, & & -0.00002366617753349105, -0.00024412456252884129, & & 0.00113028178539430542, 0.00708470513919789080, & & -0.02526914792327618386, -0.08006137953480093426, & & 0.16548380461475971846, 0.14688405470042110229 /) real, parameter :: d(0:51) = (/ & & 1.059601355592185731e-14, -2.71150591218550377e-13, & & 8.6514809056201638e-12, -4.6264028554286627e-10, & & 5.0815403835647104e-8, -1.76722552048141208e-5, & & 0.16286750396763997378, 2.949651820598278873e-13, & & -8.818215611676125741e-12, 3.571119876162253451e-10, & & -2.631924120993717060e-8, 4.709502795656698909e-6, & & -5.208333333333283282e-3, 7.18344107717531977e-15, & & -2.51623725588410308e-13, 8.6017784918920604e-12, & & -4.6256876614290359e-10, 5.0815343220437937e-8, & & -1.76722551764941970e-5, 0.16286750396763433767, & & 2.2327570859680094777e-13, -8.464594853517051292e-12, & & 3.563766464349055183e-10, -2.631843986737892965e-8, & & 4.709502342288659410e-6, -5.2083333332278466225e-3, & & 5.15413392842889366e-15, -2.27740238380640162e-13, & & 8.4827767197609014e-12, -4.6224753682737618e-10, & & 5.0814848128929134e-8, -1.76722547638767480e-5, & & 0.16286750396748926663, 1.7316195320192170887e-13, & & -7.971122772293919646e-12, 3.544039469911895749e-10, & & -2.631443902081701081e-8, 4.709498228695400603e-6, & & -5.2083333315143653610e-3, 3.84653681453798517e-15, & & -2.04464520778789011e-13, 8.3089298605177838e-12, & & -4.6155016158412096e-10, 5.0813263696466650e-8, & & -1.76722528311426167e-5, 0.16286750396650065930, & & 1.3797879972460878797e-13, -7.448089381011684812e-12, & & 3.512733797106959780e-10, -2.630500895563592722e-8, & & 4.709483934775839193e-6, -5.2083333227940760113e-3 /) w = abs(x) if (w < 1) then t = w * w y = ((((((a(0) * t + a(1)) * t + & & a(2)) * t + a(3)) * t + a(4)) * t + & & a(5)) * t + a(6)) * t + a(7) else if (w < 8.5) then t = w * w * 0.0625 k = int(t) t = t - (k + 0.5) k = k * 13 y = (((((((((((b(k) * t + b(k + 1)) * t + & & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & & b(k + 11)) * t + b(k + 12) else if (w < 12.5) then k = int(w) t = w - (k + 0.5) k = 14 * (k - 8) y = ((((((((((((c(k) * t + c(k + 1)) * t + & & c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + & & c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + & & c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + & & c(k + 11)) * t + c(k + 12)) * t + c(k + 13) else v = 24. / w t = v * v k = 13 * (int(t)) y = ((((((d(k) * t + d(k + 1)) * t + & & d(k + 2)) * t + d(k + 3)) * t + d(k + 4)) * t + & & d(k + 5)) * t + d(k + 6)) * sqrt(v) theta = (((((d(k + 7) * t + d(k + 8)) * t + & & d(k + 9)) * t + d(k + 10)) * t + d(k + 11)) * t + & & d(k + 12)) * v - .25*pi y = y * cos(w + theta) end if j0 = y end function j0 !> FIXME : Add documentation elemental function j1(x) ! Bessel J_1(x) function devided by x in double precision use constants, only: pi implicit none real :: j1 real, intent(in) :: x integer :: k real :: w, t, y, theta, v real, parameter :: a(0:7) = (/ & & -0.00000000000014810349, 0.00000000003363594618, & & -0.00000000565140051697, 0.00000067816840144764, & & -0.00005425347222188379, 0.00260416666666662438, & & -0.06249999999999999799, 0.49999999999999999998 /) real, parameter :: b(0:64) = (/ & & 0.00000000000243721316, -0.00000000009400554763, & & 0.00000000306053389980, -0.00000008287270492518, & & 0.00000183020515991344, -0.00003219783841164382, & & 0.00043795830161515318, -0.00442952351530868999, & & 0.03157908273375945955, -0.14682160488052520107, & & 0.39309619054093640008, -0.47952808215101070280, & & 0.14148999344027125140, 0.00000000000182119257, & & -0.00000000006862117678, 0.00000000217327908360, & & -0.00000005693592917820, 0.00000120771046483277, & & -0.00002020151799736374, 0.00025745933218048448, & & -0.00238514907946126334, 0.01499220060892984289, & & -0.05707238494868888345, 0.10375225210588234727, & & -0.02721551202427354117, -0.06420643306727498985, & & 0.000000000001352611196, -0.000000000049706947875, & & 0.000000001527944986332, -0.000000038602878823401, & & 0.000000782618036237845, -0.000012349994748451100, & & 0.000145508295194426686, -0.001203649737425854162, & & 0.006299092495799005109, -0.016449840761170764763, & & 0.002106328565019748701, 0.058527410006860734650, & & -0.031896615709705053191, 0.000000000000997982124, & & -0.000000000035702556073, 0.000000001062332772617, & & -0.000000025779624221725, 0.000000496382962683556, & & -0.000007310776625173004, 0.000078028107569541842, & & -0.000550624088538081113, 0.002081442840335570371, & & -0.000771292652260286633, -0.019541271866742634199, & & 0.033361194224480445382, 0.017516628654559387164, & & 0.000000000000731050661, -0.000000000025404499912, & & 0.000000000729360079088, -0.000000016915375004937, & & 0.000000306748319652546, -0.000004151324014331739, & & 0.000038793392054271497, -0.000211180556924525773, & & 0.000274577195102593786, 0.003378676555289966782, & & -0.013842821799754920148, -0.002041834048574905921, & & 0.032167266073736023299 /) real, parameter :: c(0:69) = (/ & & -0.00000000001185964494, 0.00000000039110295657, & & 0.00000000180385519493, -0.00000005575391345723, & & -0.00000018635897017174, 0.00000542738239401869, & & 0.00001181490114244279, -0.00033000319398521070, & & -0.00037717832892725053, 0.01070685852970608288, & & 0.00356629346707622489, -0.13524776185998074716, & & 0.00980725611657523952, 0.27312196367405374425, & & -0.00000000003029591097, 0.00000000009259293559, & & 0.00000000496321971223, -0.00000001518137078639, & & -0.00000057045127595547, 0.00000171237271302072, & & 0.00004271400348035384, -0.00012152454198713258, & & -0.00184155714921474963, 0.00462994691003219055, & & 0.03671737063840232452, -0.06863857568599167175, & & -0.21090395092505707655, 0.16126443075752985095, & & -0.00000000002197602080, -0.00000000027659100729, & & 0.00000000374295124827, 0.00000003684765777023, & & -0.00000045072801091574, -0.00000327941630669276, & & 0.00003571371554516300, 0.00017664005411843533, & & -0.00165119297594774104, -0.00485925381792986774, & & 0.03593306985381680131, 0.04997877588191962563, & & -0.22913866929783936544, -0.07885001422733148814, & & 0.00000000000516292316, -0.00000000039445956763, & & -0.00000000066220021263, 0.00000005511286218639, & & 0.00000005012579400780, -0.00000522111059203425, & & -0.00000134311394455105, 0.00030612891890766805, & & -0.00007103391195326182, -0.00949316714311443491, & & 0.00455036998246516948, 0.11540391585989614784, & & -0.04779493761902840455, -0.22837862066532347460, & & 0.00000000002697817493, -0.00000000016633326949, & & -0.00000000433134860350, 0.00000002508404686362, & & 0.00000048528284780984, -0.00000258267851112118, & & -0.00003521049080466759, 0.00016566324273339952, & & 0.00146474737522491617, -0.00565140892697147306, & & -0.02833882055679300400, 0.07580744376982855057, & & 0.16012275906960187978, -0.16548380461475971845 /) real, parameter :: d(0:51) = (/ & & -1.272346002224188092e-14, 3.370464692346669075e-13, & & -1.144940314335484869e-11, 6.863141561083429745e-10, & & -9.491933932960924159e-8, 5.301676561445687562e-5, & & 0.1628675039676399740, -3.652982212914147794e-13, & & 1.151126750560028914e-11, -5.165585095674343486e-10, & & 4.657991250060549892e-8, -1.186794704692706504e-5, & & 1.562499999999994026e-2, -8.713069680903981555e-15, & & 3.140780373478474935e-13, -1.139089186076256597e-11, & & 6.862299023338785566e-10, -9.491926788274594674e-8, & & 5.301676558106268323e-5, 0.1628675039676466220, & & -2.792555727162752006e-13, 1.108650207651756807e-11, & & -5.156745588549830981e-10, 4.657894859077370979e-8, & & -1.186794650130550256e-5, 1.562499999987299901e-2, & & -6.304859171204770696e-15, 2.857249044208791652e-13, & & -1.124956921556753188e-11, 6.858482894906716661e-10, & & -9.491867953516898460e-8, 5.301676509057781574e-5, & & 0.1628675039678191167, -2.185193490132496053e-13, & & 1.048820673697426074e-11, -5.132819367467680132e-10, & & 4.657409437372994220e-8, -1.186794150862988921e-5, & & 1.562499999779270706e-2, -4.740417209792009850e-15, & & 2.578715253644144182e-13, -1.104148898414138857e-11, & & 6.850134201626289183e-10, -9.491678234174919640e-8, & & 5.301676277588728159e-5, 0.1628675039690033136, & & -1.755122057493842290e-13, 9.848723331445182397e-12, & & -5.094535425482245697e-10, 4.656255982268609304e-8, & & -1.186792402114394891e-5, 1.562499998712198636e-2 /) w = abs(x) if (w < 1) then t = w * w ! y = (((((((a(0) * t + a(1)) * t + & ! & a(2)) * t + a(3)) * t + a(4)) * t + & ! & a(5)) * t + a(6)) * t + a(7)) * w y = ((((((a(0) * t + a(1)) * t + & & a(2)) * t + a(3)) * t + a(4)) * t + & & a(5)) * t + a(6)) * t + a(7) else if (w < 8.5) then t = w * w * 0.0625 k = int(t) t = t - (k + 0.5) k = k * 13 ! y = ((((((((((((b(k) * t + b(k + 1)) * t + & ! & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & ! & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & ! & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & ! & b(k + 11)) * t + b(k + 12)) * w y = (((((((((((b(k) * t + b(k + 1)) * t + & & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & & b(k + 11)) * t + b(k + 12) else if (w < 12.5) then k = int(w) t = w - (k + 0.5) k = 14 * (k - 8) ! y = ((((((((((((c(k) * t + c(k + 1)) * t + & ! & c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + & ! & c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + & ! & c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + & ! & c(k + 11)) * t + c(k + 12)) * t + c(k + 13) y = (((((((((((((c(k) * t + c(k + 1)) * t + & & c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + & & c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + & & c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + & & c(k + 11)) * t + c(k + 12)) * t + c(k + 13)) / w else v = 24. / w t = v * v k = 13 * (int(t)) y = ((((((d(k) * t + d(k + 1)) * t + & & d(k + 2)) * t + d(k + 3)) * t + d(k + 4)) * t + & & d(k + 5)) * t + d(k + 6)) * sqrt(v) theta = (((((d(k + 7) * t + d(k + 8)) * t + & & d(k + 9)) * t + d(k + 10)) * t + d(k + 11)) * t + & & d(k + 12)) * v - .25*pi y = y * sin(w + theta) ! y=y/w end if ! if (x < 0) y = -y j1 = y end function j1 !-------------------------------------------------------------------------! ! double precision error functions (derf.f) ! ! http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html ! !-------------------------------------------------------------------------! ! Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp). ! ! You may use, copy, modify this code for any purpose and ! ! without fee. You may distribute this ORIGINAL package. ! !-------------------------------------------------------------------------! ! Modified by Ryusuke NUMATA 2008/06/27 ! ! to fit F90 format ! !-------------------------------------------------------------------------! ! these routines are declared as real, but can be promoted to ! double presicion using compiler option. ! this is so because constatns used are only double precision, ! and cannot be promoted to quad-precision. ! differences from GNU's gfortran besj0 and besj1 are ! |err| < 1.e-15 !> FIXME : Add documentation elemental function erf_ext(x) implicit none ! error function in double precision real :: erf_ext real, intent(in) :: x integer :: k real :: w, t, y real, parameter :: a(0:64) = (/ & & 0.00000000005958930743, -0.00000000113739022964, & & 0.00000001466005199839, -0.00000016350354461960, & & 0.00000164610044809620, -0.00001492559551950604, & & 0.00012055331122299265, -0.00085483269811296660, & & 0.00522397762482322257, -0.02686617064507733420, & & 0.11283791670954881569, -0.37612638903183748117, & & 1.12837916709551257377, 0.00000000002372510631, & & -0.00000000045493253732, 0.00000000590362766598, & & -0.00000006642090827576, 0.00000067595634268133, & & -0.00000621188515924000, 0.00005103883009709690, & & -0.00037015410692956173, 0.00233307631218880978, & & -0.01254988477182192210, 0.05657061146827041994, & & -0.21379664776456006580, 0.84270079294971486929, & & 0.00000000000949905026, -0.00000000018310229805, & & 0.00000000239463074000, -0.00000002721444369609, & & 0.00000028045522331686, -0.00000261830022482897, & & 0.00002195455056768781, -0.00016358986921372656, & & 0.00107052153564110318, -0.00608284718113590151, & & 0.02986978465246258244, -0.13055593046562267625, & & 0.67493323603965504676, 0.00000000000382722073, & & -0.00000000007421598602, 0.00000000097930574080, & & -0.00000001126008898854, 0.00000011775134830784, & & -0.00000111992758382650, 0.00000962023443095201, & & -0.00007404402135070773, 0.00050689993654144881, & & -0.00307553051439272889, 0.01668977892553165586, & & -0.08548534594781312114, 0.56909076642393639985, & & 0.00000000000155296588, -0.00000000003032205868, & & 0.00000000040424830707, -0.00000000471135111493, & & 0.00000005011915876293, -0.00000048722516178974, & & 0.00000430683284629395, -0.00003445026145385764, & & 0.00024879276133931664, -0.00162940941748079288, & & 0.00988786373932350462, -0.05962426839442303805, & & 0.49766113250947636708 /) real, parameter :: b(0:64) = (/ & & -0.00000000029734388465, 0.00000000269776334046, & & -0.00000000640788827665, -0.00000001667820132100, & & -0.00000021854388148686, 0.00000266246030457984, & & 0.00001612722157047886, -0.00025616361025506629, & & 0.00015380842432375365, 0.00815533022524927908, & & -0.01402283663896319337, -0.19746892495383021487, & & 0.71511720328842845913, -0.00000000001951073787, & & -0.00000000032302692214, 0.00000000522461866919, & & 0.00000000342940918551, -0.00000035772874310272, & & 0.00000019999935792654, 0.00002687044575042908, & & -0.00011843240273775776, -0.00080991728956032271, & & 0.00661062970502241174, 0.00909530922354827295, & & -0.20160072778491013140, 0.51169696718727644908, & & 0.00000000003147682272, -0.00000000048465972408, & & 0.00000000063675740242, 0.00000003377623323271, & & -0.00000015451139637086, -0.00000203340624738438, & & 0.00001947204525295057, 0.00002854147231653228, & & -0.00101565063152200272, 0.00271187003520095655, & & 0.02328095035422810727, -0.16725021123116877197, & & 0.32490054966649436974, 0.00000000002319363370, & & -0.00000000006303206648, -0.00000000264888267434, & & 0.00000002050708040581, 0.00000011371857327578, & & -0.00000211211337219663, 0.00000368797328322935, & & 0.00009823686253424796, -0.00065860243990455368, & & -0.00075285814895230877, 0.02585434424202960464, & & -0.11637092784486193258, 0.18267336775296612024, & & -0.00000000000367789363, 0.00000000020876046746, & & -0.00000000193319027226, -0.00000000435953392472, & & 0.00000018006992266137, -0.00000078441223763969, & & -0.00000675407647949153, 0.00008428418334440096, & & -0.00017604388937031815, -0.00239729611435071610, & & 0.02064129023876022970, -0.06905562880005864105, & & 0.09084526782065478489 /) w = abs(x) if (w < 2.2) then t = w * w k = int(t) t = t - k k = k * 13 y = ((((((((((((a(k) * t + a(k + 1)) * t + & & a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + & & a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + & & a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + & & a(k + 11)) * t + a(k + 12)) * w else if (w < 6.9) then k = int(w) t = w - k k = 13 * (k - 2) y = (((((((((((b(k) * t + b(k + 1)) * t + & & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & & b(k + 11)) * t + b(k + 12) y = y * y y = y * y y = y * y y = 1. - y * y else y = 1. end if if (x < 0.) y = -y erf_ext = y end function erf_ext !-------------------------------------------------------------------------! # elif SPFUNC == _SPNAG_ /* if SPFUNC != _SPLOCAL_ */ !> FIXME : Add documentation function j0 (x) implicit none real :: j0 real, intent(in) :: x # if NAG_PREC == _NAGDBLE_ real (kind=nag_kind), external :: s17aef ! bessel J0 double # elif NAG_PREC == _NAGSNGL_ real (kind=nag_kind), external :: s17aee ! bessel J0 single # endif /* NAG_PREC */ real (kind=nag_kind) :: xarg xarg=x # if NAG_PREC == _NAGDBLE_ j0=s17aef(xarg,ifail) # elif NAG_PREC == _NAGSNGL_ j0=s17aee(xarg,ifail) # endif /* NAG_PREC */ end function j0 !> FIXME : Add documentation function j1 (x) real :: j1 real, intent(in) :: x # if NAG_PREC == _NAGDBLE_ real (kind=nag_kind), external :: s17aff ! bessel J1 double # elif NAG_PREC == _NAGSNGL_ real (kind=nag_kind), external :: s17afe ! bessel J1 single # endif /* NAG_PREC */ real (kind=nag_kind) :: xarg xarg=x # if NAG_PREC == _NAGDBLE_ j1=s17aff(xarg,ifail) # elif NAG_PREC == _NAGSNGL_ j1=s17afe(xarg,ifail) # endif /* NAG_PREC */ if(x == 0.) then j1=.5 else j1=j1/x endif end function j1 !> FIXME : Add documentation function erf_ext (x) implicit none real :: erf_ext real, intent(in) :: x # if NAG_PREC == _NAGDBLE_ real (kind=nag_kind), external :: s15aef ! error function double # elif NAG_PREC == _NAGSNGL_ real (kind=nag_kind), external :: s15aee ! error function single # endif /* NAG_PREC */ real (kind=nag_kind) :: xarg xarg=x # if NAG_PREC == _NAGDBLE_ erf_ext=s15aef(xarg,ifail) # elif NAG_PREC == _NAGSNGL_ erf_ext=s15aee(xarg,ifail) # endif /* NAG_PREC */ end function erf_ext # elif SPFUNC == _SPF200X_ !> FIXME : Add documentation elemental function sj0 (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj0 sj0 = bessel_j0(x) end function sj0 !> FIXME : Add documentation elemental function dj0 (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj0 dj0 = bessel_j0(x) end function dj0 # ifdef QUAD_PRECISION !> FIXME : Add documentation elemental function qj0 (x) implicit none real (kind=kind_rq), intent (in) :: x real (kind=kind_rq) :: qj0 qj0 = bessel_j0(x) end function qj0 # endif !> FIXME : Add documentation elemental function sj1 (x) use warning_helpers, only: is_zero implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj1 if (is_zero(x)) then sj1 = 0.5_kind_rs else sj1 = bessel_j1(x) / x end if end function sj1 !> FIXME : Add documentation elemental function dj1 (x) use warning_helpers, only: is_zero implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj1 if (is_zero(x)) then dj1 = 0.5_kind_rd else dj1 = bessel_j1(x) / x end if end function dj1 # ifdef QUAD_PRECISION !> FIXME : Add documentation elemental function qj1 (x) implicit none real (kind=kind_rq), intent (in) :: x real (kind=kind_rq) :: qj1 if (x == 0.0_kind_rq) then qj1 = 0.5_kind_rq else qj1 = bessel_j1(x) / x end if end function qj1 # endif !> FIXME : Add documentation elemental function serf_ext (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: serf_ext serf_ext = erf(x) end function serf_ext !> FIXME : Add documentation elemental function derf_ext (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: derf_ext derf_ext = erf(x) end function derf_ext # ifdef QUAD_PRECISION !> FIXME : Add documentation elemental function qerf_ext (x) implicit none real (kind=kind_rq), intent (in) :: x real (kind=kind_rq) :: qerf_ext qerf_ext = erf(x) end function qerf_ext # endif !> FIXME : Add documentation elemental function slgamma_ext (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs):: slgamma_ext slgamma_ext = log_gamma(x) end function slgamma_ext !> FIXME : Add documentation elemental function dlgamma_ext (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd):: dlgamma_ext dlgamma_ext = log_gamma(x) end function dlgamma_ext # ifdef QUAD_PRECISION !> FIXME : Add documentation elemental function qlgamma_ext (x) implicit none real (kind=kind_rq), intent (in) :: x real (kind=kind_rq):: qlgamma_ext qlgamma_ext = log_gamma(x) end function qlgamma_ext # endif # else /* if SPFUNC != _SPLOCAL_ && SPFUNC != _SPNAG_ && ! SPFUNC != _SPF200X_ */ # if FCOMPILER == _PATHSCALE_ !> FIXME : Add documentation elemental function sj0 (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj0 sj0 = besj0(x) end function sj0 !> FIXME : Add documentation elemental function dj0 (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj0 dj0 = dbesj0(x) end function dj0 !> FIXME : Add documentation elemental function sj1 (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj1 if (x == 0.0) then sj1 = 0.5 else sj1 = besj1(x) / x end if end function sj1 !> FIXME : Add documentation elemental function dj1 (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj1 if (x == 0.0) then dj1 = 0.5 else dj1 = dbesj1(x) / x end if end function dj1 !> FIXME : Add documentation elemental function serf_ext (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: serf_ext serf_ext = erf(x) end function serf_ext !> FIXME : Add documentation elemental function derf_ext (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: derf_ext derf_ext = derf(x) end function derf_ext # elif (FCOMPILER == _INTEL_ || FCOMPILER == _IFX_ )/* not _PATHSCALE_ */ !> FIXME : Add documentation function sj0 (x) use ifport, only: besj0 ! single precision version (not generic) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj0 sj0 = besj0(x) end function sj0 !> FIXME : Add documentation function dj0 (x) use ifport, only: dbesj0 implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj0 dj0 = dbesj0(x) end function dj0 !> FIXME : Add documentation function sj1 (x) use ifport, only: besj1 ! single precision version (not generic) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj1 if (x == 0.0) then sj1 = 0.5 else sj1 = besj1(x) / x end if end function sj1 !> FIXME : Add documentation function dj1 (x) use ifport, only: dbesj1 implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj1 if (x == 0.0) then dj1 = 0.5 else dj1 = dbesj1(x) / x end if end function dj1 !> FIXME : Add documentation function serf_ext (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: serf_ext serf_ext = erf(x) ! erf is generic end function serf_ext !> FIXME : Add documentation function derf_ext (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: derf_ext derf_ext = erf(x) ! erf is generic end function derf_ext # elif FCOMPILER == _GFORTRAN_ /* not _INTEL_, _PATHSCALE_ */ !> FIXME : Add documentation elemental function sj0 (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj0 sj0 = besj0(x) end function sj0 !> FIXME : Add documentation elemental function dj0 (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj0 dj0 = besj0(x) end function dj0 # ifdef QUAD_PRECISION !> FIXME : Add documentation elemental function qj0 (x) implicit none real (kind=kind_rq), intent (in) :: x real (kind=kind_rq) :: qj0 qj0 = besj0(x) end function qj0 # endif !> FIXME : Add documentation elemental function sj1 (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj1 if (x == 0.0) then sj1 = 0.5 else sj1 = besj1(x) / x end if end function sj1 !> FIXME : Add documentation elemental function dj1 (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj1 if (x == 0.0) then dj1 = 0.5 else dj1 = besj1(x) / x end if end function dj1 # ifdef QUAD_PRECISION !> FIXME : Add documentation elemental function qj1 (x) implicit none real (kind=kind_rq), intent (in) :: x real (kind=kind_rq) :: qj1 if (x == 0.0_kind_rq) then qj1 = 0.5_kind_rq else qj1 = besj1(x) / x end if end function qj1 # endif !> FIXME : Add documentation elemental function serf_ext (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: serf_ext serf_ext = erf(x) end function serf_ext !> FIXME : Add documentation elemental function derf_ext (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: derf_ext derf_ext = erf(x) end function derf_ext # ifdef QUAD_PRECISION !> FIXME : Add documentation elemental function qerf_ext (x) implicit none real (kind=kind_rq), intent (in) :: x real (kind=kind_rq) :: qerf_ext qerf_ext = erf(x) end function qerf_ext # endif # elif FCOMPILER == _G95_ /* not _GFORTRAN_, _INTEL_, _PATHSCALE_ */ !> FIXME : Add documentation elemental function sj0 (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj0 sj0 = besj0(x) end function sj0 !> FIXME : Add documentation elemental function dj0 (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj0 dj0 = dbesj0(x) end function dj0 !> FIXME : Add documentation elemental function sj1 (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rs) :: sj1 if (x == 0.0) then sj1 = 0.5 else sj1 = besj1(x) / x end if end function sj1 !> FIXME : Add documentation elemental function dj1 (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dj1 if (x == 0.0) then dj1 = 0.5 else dj1 = dbesj1(x) / x end if end function dj1 !> FIXME : Add documentation function serf_ext(x) implicit none real (kind=kind_rs), intent(in) :: x real (kind=kind_rs) :: serf_ext serf_ext=erf(real(x)) end function serf_ext !> FIXME : Add documentation function derf_ext(x) implicit none real (kind=kind_rd), intent(in) :: x real (kind=kind_rd):: derf_ext derf_ext=erf(real(x)) end function derf_ext !> FIXME : Add documentation function slgamma_ext (x) implicit none real (kind=kind_rs), intent (in) :: x real (kind=kind_rd):: slgamma_ext slgamma_ext = algama(real(x)) end function slgamma_ext !> FIXME : Add documentation function dlgamma_ext (x) implicit none real (kind=kind_rd), intent (in) :: x real (kind=kind_rd) :: dlgamma_ext dlgamma_ext = algama(real(x)) end function dlgamma_ext # elif FCOMPILER == _PGI_ /* not _G95_, _GFORTRAN_, _INTEL_, _PATHSCALE_ */ !> FIXME : Add documentation elemental function sj1 (x) implicit none real (kind=4), intent (in) :: x real (kind=4) :: sj1 interface besj1 elemental function besj1(x) real (kind=4), intent (in) :: x real (kind=4) :: besj1 end function besj1 end interface if (x == 0.0) then sj1 = 0.5 else sj1 = besj1(x) / x end if end function sj1 !> FIXME : Add documentation elemental function dj1 (x) implicit none real (kind=8), intent (in) :: x real (kind=8) :: dj1 interface besj1 elemental function dbesj1(x) real (kind=8), intent (in) :: x real (kind=8) :: dbesj1 end function dbesj1 end interface if (x == 0.0) then dj1 = 0.5 else dj1 = besj1(x) / x end if end function dj1 # elif FCOMPILER == _FUJ_ /* not _PGI_, _G95_, _GFORTRAN_, _INTEL_, _PATHSCALE */ ! SSL2 !> FIXME : Add documentation function j0 (x) implicit none real :: j0 real, intent(in) :: x integer :: icon interface bj0 elemental subroutine bj0(x,j0,icon) real (kind=4), intent(in) :: x real (kind=4), intent(out) :: j0 integer, intent(out) :: icon end subroutine bj0 elemental subroutine dbj0(x,j0,icon) real (kind=8), intent(in) :: x real (kind=8), intent(out) :: j0 integer, intent(out) :: icon end subroutine dbj0 end interface call bj0(x,j0,icon) ! if (icon /= 0 ) then ! end if end function j0 ! SSL2 !> FIXME : Add documentation function j1 (x) implicit none real :: j1 real, intent(in) :: x integer :: icon interface bj1 elemental subroutine bj1(x,j1,icon) real (kind=4), intent(in) :: x real (kind=4), intent(out) :: j1 integer, intent(out) :: icon end subroutine bj1 elemental subroutine dbj1(x,j1,icon) real (kind=8), intent(in) :: x real (kind=8), intent(out) :: j1 integer, intent(out) :: icon end subroutine dbj1 end interface call bj1(x,j1,icon) ! if (icon /= 0) then ! end if if (x == 0.) then j1=.5 else j1=j1/x end if end function j1 !> FIXME : Add documentation function erf_ext(x) implicit none real :: erf_ext real, intent(in) :: x erf_ext=erf(x) end function erf_ext !> FIXME : Add documentation function lgamma_ext(x) implicit none real :: lgamma_ext real, intent(in) :: x lgamma_ext=lgamma(x) end function lgamma_ext # elif FCOMPILER == _NEC_ /* not _FUJ_, _PGI_, _G95_, _GFORTRAN_, _INTEL_, _PATHSCALE */ ! ASL !> FIXME : Add documentation function j0s (x) implicit none real :: j0s real, intent(in) :: x real :: xi(1),xo(1) integer :: ierr interface bj0 subroutine vibj0x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=4), intent(in) :: xi(1) real (kind=4), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine vibj0x subroutine wibj0x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=8), intent(in) :: xi(1) real (kind=8), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine wibj0x end interface xi(1)=x call bj0(1,xi,xo,ierr) j0s=xo(1) ! if (ierr /= 0 ) then ! end if end function j0s function j0a (x) implicit none real, intent(in) :: x(:) real :: j0a(size(x,1)) integer :: ierr interface bj0 subroutine vibj0x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=4), intent(in) :: xi(1) real (kind=4), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine vibj0x subroutine wibj0x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=8), intent(in) :: xi(1) real (kind=8), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine wibj0x end interface call bj0(size(x,1),x,j0a,ierr) ! if (ierr /= 0 ) then ! end if end function j0a ! ASL !> FIXME : Add documentation function j1s (x) implicit none real :: j1s real, intent(in) :: x real :: xi(1),xo(1) integer :: ierr interface bj1 subroutine vibj1x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=4), intent(in) :: xi(1) real (kind=4), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine vibj1x subroutine wibj1x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=8), intent(in) :: xi(1) real (kind=8), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine wibj1x end interface xi(1)=x call bj1(1,xi,xo,ierr) ! if (ierr /= 0) then ! end if if (x == 0.) then j1s=.5 else j1s=xo(1)/x end if end function j1s function j1a (x) implicit none real, intent(in) :: x(:) real :: j1a(size(x,1)) integer :: ierr interface bj1 subroutine vibj1x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=4), intent(in) :: xi(1) real (kind=4), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine vibj1x subroutine wibj1x(nv,xi,xo,ierr) integer, intent(in) :: nv real (kind=8), intent(in) :: xi(1) real (kind=8), intent(out) :: xo(1) integer, intent(out) :: ierr end subroutine wibj1x end interface call bj1(size(x,1),x,j1a,ierr) ! if (ierr /= 0) then ! end if where (x == 0.) j1a=.5 elsewhere j1a=j1a/x end where end function j1a !> FIXME : Add documentation function erf_ext(x) implicit none real :: erf_ext real, intent(in) :: x erf_ext=erf(x) end function erf_ext !> FIXME : Add documentation function lgamma_ext(x) implicit none real :: lgamma_ext real, intent(in) :: x lgamma_ext=lgamma(x) end function lgamma_ext # endif /* if FCOMPILER */ # endif /* if SPFUNC */ # if SPFUNC != _SPF200X_ # if (SPFUNC == _SPLOCAL_ \ || ( FCOMPILER != _G95_ && FCOMPILER != _FUJ_ && FCOMPILER != _NEC_ ) ) ! TT: We may be able to exclude other compilers ! TT: but ifort, PGI, pathscale don't seem to have it as I quickly checked ! TT: NAG may have it, which we have to check !-------------------------------------------------------------------------! ! double precision log gamma function (lgamma.f) ! ! http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html ! !-------------------------------------------------------------------------! ! Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp). ! ! You may use, copy, modify this code for any purpose and ! ! without fee. You may distribute this ORIGINAL package. ! !-------------------------------------------------------------------------! ! Modified by Ryusuke NUMATA 2008/06/27 ! ! to fit F90 format ! !-------------------------------------------------------------------------! !> FIXME : Add documentation (or tidy existing) elemental function lgamma_ext(x) use constants, only: pi implicit none real :: lgamma_ext real, intent(in) :: x integer :: k real :: w, v, t, y real, parameter :: a(0:21) = (/ & & 0.00009967270908702825, -0.00019831672170162227, & & -0.00117085315349625822, 0.00722012810948319552, & & -0.00962213009367802970, -0.04219772092994235254, & & 0.16653861065243609743, -0.04200263501129018037, & & -0.65587807152061930091, 0.57721566490153514421, & & 0.99999999999999999764, 0.00004672097259011420, & & -0.00006812300803992063, -0.00132531159076610073, & & 0.00733521178107202770, -0.00968095666383935949, & & -0.04217642811873540280, 0.16653313644244428256, & & -0.04200165481709274859, -0.65587818792782740945, & & 0.57721567315209190522, 0.99999999973565236061 /) real, parameter :: b(0:97) = (/ & & -0.00000000004587497028, 0.00000000019023633960, & & 0.00000000086377323367, 0.00000001155136788610, & & -0.00000002556403058605, -0.00000015236723372486, & & -0.00000316805106385740, 0.00000122903704923381, & & 0.00002334372474572637, 0.00111544038088797696, & & 0.00344717051723468982, 0.03198287045148788384, & & -0.32705333652955399526, 0.40120442440953927615, & & -0.00000000005184290387, -0.00000000083355121068, & & -0.00000000256167239813, 0.00000001455875381397, & & 0.00000013512178394703, 0.00000029898826810905, & & -0.00000358107254612779, -0.00002445260816156224, & & -0.00004417127762011821, 0.00112859455189416567, & & 0.00804694454346728197, 0.04919775747126691372, & & -0.24818372840948854178, 0.11071780856646862561, & & 0.00000000030279161576, 0.00000000160742167357, & & -0.00000000405596009522, -0.00000005089259920266, & & -0.00000002029496209743, 0.00000135130272477793, & & 0.00000391430041115376, -0.00002871505678061895, & & -0.00023052137536922035, 0.00045534656385400747, & & 0.01153444585593040046, 0.07924014651650476036, & & -0.12152192626936502982, -0.07916438300260539592, & & -0.00000000050919149580, -0.00000000115274986907, & & 0.00000001237873512188, 0.00000002937383549209, & & -0.00000030621450667958, -0.00000077409414949954, & & 0.00000816753874325579, 0.00002412433382517375, & & -0.00026061217606063700, -0.00091000087658659231, & & 0.01068093850598380797, 0.11395654404408482305, & & 0.07209569059984075595, -0.10971041451764266684, & & 0.00000000040119897187, -0.00000000013224526679, & & -0.00000001002723190355, 0.00000002569249716518, & & 0.00000020336011868466, -0.00000118097682726060, & & -0.00000300660303810663, 0.00004402212897757763, & & -0.00001462405876235375, -0.00164873795596001280, & & 0.00513927520866443706, 0.13843580753590579416, & & 0.32730190978254056722, 0.08588339725978624973, & & -0.00000000015413428348, 0.00000000064905779353, & & 0.00000000160702811151, -0.00000002655645793815, & & 0.00000007619544277956, 0.00000047604380765353, & & -0.00000490748870866195, 0.00000821513040821212, & & 0.00014804944070262948, -0.00122152255762163238, & & -0.00087425289205498532, 0.14438703699657968310, & & 0.61315889733595543766, 0.55513708159976477557, & & 0.00000000001049740243, -0.00000000025832017855, & & 0.00000000139591845075, -0.00000000021177278325, & & -0.00000005082950464905, 0.00000037801785193343, & & -0.00000073982266659145, -0.00001088918441519888, & & 0.00012491810452478905, -0.00049171790705139895, & & -0.00425707089448266460, 0.13595080378472757216, & & 0.89518356003149514744, 1.31073912535196238583 /) real, parameter :: c(0:64) = (/ & & 0.0000000116333640008, -0.0000000833156123568, & & 0.0000003832869977018, -0.0000015814047847688, & & 0.0000065010672324100, -0.0000274514060128677, & & 0.0001209015360925566, -0.0005666333178228163, & & 0.0029294103665559733, -0.0180340086069185819, & & 0.1651788780501166204, 1.1031566406452431944, & & 1.2009736023470742248, 0.0000000013842760642, & & -0.0000000069417501176, 0.0000000342976459827, & & -0.0000001785317236779, 0.0000009525947257118, & & -0.0000052483007560905, 0.0000302364659535708, & & -0.0001858396115473822, 0.0012634378559425382, & & -0.0102594702201954322, 0.1243625515195050218, & & 1.3888709263595291174, 2.4537365708424422209, & & 0.0000000001298977078, -0.0000000008029574890, & & 0.0000000049454846150, -0.0000000317563534834, & & 0.0000002092136698089, -0.0000014252023958462, & & 0.0000101652510114008, -0.0000774550502862323, & & 0.0006537746948291078, -0.0066014912535521830, & & 0.0996711934948138193, 1.6110931485817511402, & & 3.9578139676187162939, 0.0000000000183995642, & & -0.0000000001353537034, 0.0000000009984676809, & & -0.0000000076346363974, 0.0000000599311464148, & & -0.0000004868554120177, 0.0000041441957716669, & & -0.0000377160856623282, 0.0003805693126824884, & & -0.0045979851178130194, 0.0831422678749791178, & & 1.7929113303999329439, 5.6625620598571415285, & & 0.0000000000034858778, -0.0000000000297587783, & & 0.0000000002557677575, -0.0000000022705728282, & & 0.0000000207024992450, -0.0000001954426390917, & & 0.0000019343161886722, -0.0000204790249102570, & & 0.0002405181940241215, -0.0033842087561074799, & & 0.0713079483483518997, 1.9467574842460867884, & & 7.5343642367587329552 /) real, parameter :: d(0:6) = (/ & & -0.00163312359200500807, 0.00083644533703385956, & & -0.00059518947575728181, 0.00079365057505415415, & & -0.00277777777735463043, 0.08333333333333309869, & & 0.91893853320467274178 /) w = x if (x < 0.) w = 1. - x if (w < 0.5) then k = 0 if (w >= 0.25) k = 11 y = ((((((((((a(k) * w + a(k + 1)) * w + & & a(k + 2)) * w + a(k + 3)) * w + a(k + 4)) * w + & & a(k + 5)) * w + a(k + 6)) * w + a(k + 7)) * w + & & a(k + 8)) * w + a(k + 9)) * w + a(k + 10)) * w y = -log(y) else if (w < 3.5) then t = w - 4.5 / (w + 0.5) k = int(t + 4) t = t - (k - 3.5) k = k * 14 y = ((((((((((((b(k) * t + b(k + 1)) * t + & & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & & b(k + 11)) * t + b(k + 12)) * t + b(k + 13) else if (w < 8.) then k = (int(w)) - 3 t = w - (k + 3.5) k = k * 13 y = (((((((((((c(k) * t + c(k + 1)) * t + & & c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + & & c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + & & c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + & & c(k + 11)) * t + c(k + 12) else v = 1. / w t = v * v y = (((((d(0) * t + d(1)) * t + d(2)) * t + & & d(3)) * t + d(4)) * t + d(5)) * v + d(6) y = y + ((w - 0.5) * log(w) - w) end if if (x < 0.) y = log(pi / sin(pi * x)) - y lgamma_ext = y end function lgamma_ext # endif /* if (SPFUNC == _SPLOCAL_ || ( FCOMPILER != _G95_ && FCOMPILER != _FUJ_ ) */ # endif /* if (SPFUNC != _SPF200X_ */ end module spfunc