|
! T. Fukushima, 2015, |
|
! Precise and fast computation of complete elliptic integrals by piecewise minimax rational function approximation, |
|
! J. Comp. Appl. Math., 282, 71-76 |
|
! |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
program xkebd |
|
! |
|
! Test driver of minimax rational function approximations |
|
! of complete elliptic integrals, K(m), E(m), B(m), and D(m) |
|
! |
|
implicit integer (i-n) |
|
implicit real*8 (a-h,o-z) |
|
|
|
iend=10 |
|
dmc=1.d0/dble(iend) |
|
|
|
write (*,"(5x,a10,a15,a25)") "m_c","single","double" |
|
|
|
do i=1,iend |
|
fmc=dmc*dble(i) |
|
sf=sceik(fmc) |
|
f=ceik(fmc) |
|
write (*,"(a5,0pf10.5,1pe15.8,1pe25.16)") "K(m)",fmc,sf,f |
|
enddo |
|
|
|
write (*,*) " " |
|
|
|
do i=1,iend |
|
fmc=dmc*dble(i) |
|
sf=sceie(fmc) |
|
f=ceie(fmc) |
|
write (*,"(a5,0pf10.5,1pe15.8,1pe25.16)") "E(m)",fmc,sf,f |
|
enddo |
|
|
|
write (*,*) " " |
|
|
|
do i=1,iend |
|
fmc=dmc*dble(i) |
|
sf=sceib(fmc) |
|
f=ceib(fmc) |
|
write (*,"(a5,0pf10.5,1pe15.8,1pe25.16)") "B(m)",fmc,sf,f |
|
enddo |
|
|
|
write (*,*) " " |
|
|
|
do i=1,iend |
|
fmc=dmc*dble(i) |
|
sf=sceid(fmc) |
|
f=ceid(fmc) |
|
write (*,"(a5,0pf10.5,1pe15.8,1pe25.16)") "D(m)",fmc,sf,f |
|
enddo |
|
|
|
end program xkebd |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function ceik(mc) |
|
! |
|
! Double precision minimax rational approximation of K(m): |
|
! the complete elliptic integral of the first kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(ceik) negative m: mc=",mc |
|
elseif(mc.gt.0.592990d0) then |
|
t=2.45694208987494165d0*mc-1.45694208987494165d0 |
|
ceik=(3703.75266375099019d0 & |
|
+t*(5462.47093231923466d0 & |
|
+t*(2744.82029097576810d0 & |
|
+t*(543.839017382099411d0 & |
|
+t*(36.2381612593459565d0 & |
|
+t*0.393188651542789784d0)))))/ & |
|
(2077.94377067058435d0 & |
|
+t*(3398.00069767755460d0 & |
|
+t*(1959.05960044399275d0 & |
|
+t*(472.794455487539279d0 & |
|
+t*(43.5464368440078942d0 & |
|
+t))))) |
|
elseif(mc.gt.0.350756d0) then |
|
t=4.12823963605439369d0*mc-1.44800482178389491d0 |
|
ceik=(4264.28203103974630d0 & |
|
+t*(6341.90978213264024d0 & |
|
+t*(3214.59187442783167d0 & |
|
+t*(642.790566685354573d0 & |
|
+t*(43.2589626155454993d0 & |
|
+t*0.475223892294445943d0)))))/ & |
|
(2125.06914237062279d0 & |
|
+t*(3479.95663350926514d0 & |
|
+t*(2006.03187933518870d0 & |
|
+t*(482.900172581418890d0 & |
|
+t*(44.1848041560412224d0 & |
|
+t))))) |
|
elseif(mc.gt.0.206924d0) then |
|
t=6.95255575949719117d0*mc-1.43865064797819679d0 |
|
ceik=(4870.25402224986382d0 & |
|
+t*(7307.18826377416591d0 & |
|
+t*(3738.29369283392307d0 & |
|
+t*(754.928587580583704d0 & |
|
+t*(51.3609902253065926d0 & |
|
+t*0.571948962277566451d0)))))/ & |
|
(2172.51745704102287d0 & |
|
+t*(3565.04737778032566d0 & |
|
+t*(2056.13612019430497d0 & |
|
+t*(493.962405117599400d0 & |
|
+t*(44.9026847057686146d0 & |
|
+t))))) |
|
elseif(mc.gt.0.121734d0) then |
|
t=11.7384669562155183d0*mc-1.42897053644793990d0 |
|
ceik=(5514.8512729127464d0 & |
|
+t*(8350.4595896779631d0 & |
|
+t*(4313.60788246750934d0 & |
|
+t*(880.27903031894216d0 & |
|
+t*(60.598720224393536d0 & |
|
+t*0.68504458747933773d0)))))/ & |
|
(2218.41682813309737d0 & |
|
+t*(3650.41829123846319d0 & |
|
+t*(2107.97379949034285d0 & |
|
+t*(505.74295207655096d0 & |
|
+t*(45.6911096775045314d0 & |
|
+t))))) |
|
elseif(mc.gt.0.071412d0) then |
|
t=19.8720241643813839d0*mc-1.41910098962680339d0 |
|
ceik=(6188.8743957372448d0 & |
|
+t*(9459.3331440432847d0 & |
|
+t*(4935.41351498551527d0 & |
|
+t*(1018.21910476032105d0 & |
|
+t*(70.981049144472361d0 & |
|
+t*0.81599895108245948d0)))))/ & |
|
(2260.73112539748448d0 & |
|
+t*(3732.66955095581621d0 & |
|
+t*(2159.68721749761492d0 & |
|
+t*(517.86964191812384d0 & |
|
+t*(46.5298955058476510d0 & |
|
+t))))) |
|
elseif(mc.gt.0.041770d0) then |
|
t=33.7359152553808785d0*mc-1.40914918021725929d0 |
|
ceik=(6879.5170681289562d0 & |
|
+t*(10615.0836403687221d0 & |
|
+t*(5594.8381504799829d0 & |
|
+t*(1167.26108955935542d0 & |
|
+t*(82.452856129147838d0 & |
|
+t*0.96592719058503951d0)))))/ & |
|
(2296.88303450660439d0 & |
|
+t*(3807.37745652028212d0 & |
|
+t*(2208.74949754945558d0 & |
|
+t*(529.79651353072921d0 & |
|
+t*(47.3844470709989137d0 & |
|
+t))))) |
|
elseif(mc.gt.0.024360d0) then |
|
t=57.4382538770821367d0*mc-1.39919586444572085d0 |
|
ceik=(7570.6827538712100d0 & |
|
+t*(11792.9392624454532d0 & |
|
+t*(6279.2661370014890d0 & |
|
+t*(1325.01058966228180d0 & |
|
+t*(94.886883830605940d0 & |
|
+t*1.13537029594409690d0)))))/ & |
|
(2324.04824540459984d0 & |
|
+t*(3869.56755306385732d0 & |
|
+t*(2252.22250562615338d0 & |
|
+t*(540.85752251676412d0 & |
|
+t*(48.2089280211559345d0 & |
|
+t))))) |
|
elseif(mc.gt.0.014165d0) then |
|
t=98.0872976949485042d0*mc-1.38940657184894556d0 |
|
ceik=(8247.2601660137746d0 & |
|
+t*(12967.7060124572914d0 & |
|
+t*(6974.7495213178613d0 & |
|
+t*(1488.54008220335966d0 & |
|
+t*(108.098282908839979d0 & |
|
+t*1.32411616748380686d0)))))/ & |
|
(2340.47337508405427d0 & |
|
+t*(3915.63324533769906d0 & |
|
+t*(2287.70677154700516d0 & |
|
+t*(550.45072377717361d0 & |
|
+t*(48.9575432570382154d0 & |
|
+t))))) |
|
elseif(mc.gt.0.008213d0) then |
|
t=168.010752688172043d0*mc-1.37987231182795699d0 |
|
ceik=(8894.2961573611293d0 & |
|
+t*(14113.7038749808951d0 & |
|
+t*(7666.5611739483371d0 & |
|
+t*(1654.60731579994159d0 & |
|
+t*(121.863474964652041d0 & |
|
+t*1.53112170837206117d0)))))/ & |
|
(2344.88618943372377d0 & |
|
+t*(3942.81065054556536d0 & |
|
+t*(2313.28396270968662d0 & |
|
+t*(558.07615380622169d0 & |
|
+t*(49.5906602613891184d0 & |
|
+t))))) |
|
elseif(mc.gt.0.d0) then |
|
t=1.d0-121.758188238159016d0*mc |
|
ceik=-log(mc*0.0625d0) & |
|
*(34813.4518336350547d0 & |
|
+t*(235.767716637974271d0 & |
|
+t*0.199792723884069485d0))/ & |
|
(69483.5736412906324d0 & |
|
+t*(614.265044703187382d0 & |
|
+t)) & |
|
-mc*(9382.53386835986099d0 & |
|
+t*(51.6478985993381223d0 & |
|
+t*0.00410754154682816898d0))/ & |
|
(37327.7262507318317d0 & |
|
+t*(408.017247271148538d0 & |
|
+t)) |
|
else |
|
write(*,"(a20,1pe15.7)") "(ceik) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(ceik) mc,K=",mc,ceik |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function ceie(mc) |
|
! |
|
! Double precision minimax rational approximation of E(m): |
|
! the complete elliptic integral of the second kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(ceie) negative m: mc=",mc |
|
elseif(mc.gt.0.566638d0) then |
|
t=2.30753965506897236d0*mc-1.30753965506897236d0 |
|
ceie=(19702.2363352671642d0 & |
|
+t*(31904.1559574281609d0 & |
|
+t*(18177.1879313824040d0 & |
|
+t*(4362.94760768571862d0 & |
|
+t*(409.975559128654710d0 & |
|
+t*10.3244775335024885d0)))))/ & |
|
(14241.2135819448616d0 & |
|
+t*(20909.9899599927367d0 & |
|
+t*(10266.4884503526076d0 & |
|
+t*(1934.86289070792954d0 & |
|
+t*(117.162100771599098d0 & |
|
+t))))) |
|
elseif(mc.gt.0.315153d0) then |
|
t=3.97638030101198879d0*mc-1.25316818100483130d0 |
|
ceie=(16317.0721393008221d0 & |
|
+t*(26627.8852140835023d0 & |
|
+t*(15129.4009798463159d0 & |
|
+t*(3574.15857605556033d0 & |
|
+t*(326.113727011739428d0 & |
|
+t*7.93163724081373477d0)))))/ & |
|
(13047.1505096551210d0 & |
|
+t*(19753.5762165922376d0 & |
|
+t*(9964.25173735060361d0 & |
|
+t*(1918.72232033637537d0 & |
|
+t*(117.670514069579649d0 & |
|
+t))))) |
|
elseif(mc.gt.0.171355d0) then |
|
t=6.95419964116329852d0*mc-1.19163687951153702d0 |
|
ceie=(13577.3850240991520d0 & |
|
+t*(22545.4744699553993d0 & |
|
+t*(12871.9137872656293d0 & |
|
+t*(3000.74575264868572d0 & |
|
+t*(263.964361648520708d0 & |
|
+t*6.08522443139677663d0)))))/ & |
|
(11717.3306408059832d0 & |
|
+t*(18431.1264424290258d0 & |
|
+t*(9619.40382323874064d0 & |
|
+t*(1904.06010727307491d0 & |
|
+t*(118.690522739531267d0 & |
|
+t))))) |
|
elseif(mc.gt.0.090670d0) then |
|
t=12.3938774245522712d0*mc-1.12375286608415443d0 |
|
ceie=(11307.9485341543712d0 & |
|
+t*(19328.6173704569489d0 & |
|
+t*(11208.6068472959372d0 & |
|
+t*(2596.54874477084334d0 & |
|
+t*(219.253495956962613d0 & |
|
+t*4.66931143174036616d0)))))/ & |
|
(10307.6837501971393d0 & |
|
+t*(16982.2450249024383d0 & |
|
+t*(9241.7604666150102d0 & |
|
+t*(1893.41905403040679d0 & |
|
+t*(120.498555754227847d0 & |
|
+t))))) |
|
elseif(mc.gt.0.046453d0) then |
|
t=22.6157360291290680d0*mc-1.05056878576113260d0 |
|
ceie=(9383.1490856819874d0 & |
|
+t*(16718.9730458676860d0 & |
|
+t*(9977.2498973537718d0 & |
|
+t*(2323.49987246555537d0 & |
|
+t*(188.618148076418837d0 & |
|
+t*3.59313532204509922d0)))))/ & |
|
(8877.1964704758383d0 & |
|
+t*(15450.0537230364062d0 & |
|
+t*(8840.2771293410661d0 & |
|
+t*(1889.13672102820913d0 & |
|
+t*(123.422125687316355d0 & |
|
+t))))) |
|
elseif(mc.gt.0.022912d0) then |
|
t=42.4790790535661187d0*mc-0.973280659275306911d0 |
|
ceie=(7719.1171817802054d0 & |
|
+t*(14521.7363804934985d0 & |
|
+t*(9045.3996063894006d0 & |
|
+t*(2149.92068078627829d0 & |
|
+t*(169.386557799782496d0 & |
|
+t*2.78515570453129137d0)))))/ & |
|
(7479.7539074698012d0 & |
|
+t*(13874.4978011497847d0 & |
|
+t*(8420.3848818926324d0 & |
|
+t*(1892.69753150329759d0 & |
|
+t*(127.802109608726363d0 & |
|
+t))))) |
|
elseif(mc.gt.0.010809d0) then |
|
t=82.6241427745187144d0*mc-0.893084359249772784d0 |
|
ceie=(6261.6095608987273d0 & |
|
+t*(12593.0874916293982d0 & |
|
+t*(8304.3265605809870d0 & |
|
+t*(2048.68391263416822d0 & |
|
+t*(159.371262600702237d0 & |
|
+t*2.18867046462858104d0)))))/ & |
|
(6156.4532048239501d0 & |
|
+t*(12283.8373999680518d0 & |
|
+t*(7979.7435857665227d0 & |
|
+t*(1903.60556312663537d0 & |
|
+t*(133.911640385965187d0 & |
|
+t))))) |
|
elseif(mc.gt.0.004841d0) then |
|
t=167.560321715817694d0*mc-0.811159517426273458d0 |
|
ceie=(4978.06146583586728d0 & |
|
+t*(10831.7178150656694d0 & |
|
+t*(7664.6703673290453d0 & |
|
+t*(1995.66437151562090d0 & |
|
+t*(156.689647694892782d0 & |
|
+t*1.75859085945198570d0)))))/ & |
|
(4935.56743322938333d0 & |
|
+t*(10694.5510113880077d0 & |
|
+t*(7506.8028283118051d0 & |
|
+t*(1918.38517009740321d0 & |
|
+t*(141.854303920116856d0 & |
|
+t))))) |
|
elseif(mc.gt.0.d0) then |
|
t=1.d0-206.568890725056806d0*mc |
|
ceie=-mc*log(mc*0.0625d0) & |
|
*(41566.6612602868736d0 & |
|
+t*(154.034981522913482d0 & |
|
+t*0.0618072471798575991d0))/ & |
|
(165964.442527585615d0 & |
|
+t*(917.589668642251803d0 & |
|
+t)) & |
|
+(132232.803956682877d0 & |
|
+t*(353.375480007017643d0 & |
|
-t*1.40105837312528026d0))/ & |
|
(132393.665743088043d0 & |
|
+t*(192.112635228732532d0 & |
|
-t)) |
|
else |
|
write(*,"(a20,1pe15.7)") "(ceie) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(ceie) mc,E=",mc,ceie |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function ceib(mc) |
|
! |
|
! Double precision minimax rational approximation of B(m): |
|
! the first associate complete elliptic integral of the second kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(ceib) negative m: mc=",mc |
|
elseif(mc.gt.0.555073d0) then |
|
t=2.24755971204264969d0*mc-1.24755971204264969d0 |
|
ceib=(2030.25011505956379d0 & |
|
+t*(3223.16236100954529d0 & |
|
+t*(1727.60635612511943d0 & |
|
+t*(361.164121995173076d0 & |
|
+t*(25.0715510300422010d0 & |
|
+t*0.280355207707726826d0)))))/ & |
|
(2420.64907902774675d0 & |
|
+t*(4034.28168313496638d0 & |
|
+t*(2327.48464880306840d0 & |
|
+t*(549.234220839203960d0 & |
|
+t*(47.9870997057202318d0 & |
|
+t))))) |
|
elseif(mc.gt.0.302367d0) then |
|
t=3.95716761770595079d0*mc-1.19651690106289522d0 |
|
ceib=(2209.26925068374373d0 & |
|
+t*(3606.58475322372526d0 & |
|
+t*(1981.37862223307242d0 & |
|
+t*(422.693774742063054d0 & |
|
+t*(29.7612810087709299d0 & |
|
+t*0.334623999861181980d0)))))/ & |
|
(2499.57898767250755d0 & |
|
+t*(4236.30953048456334d0 & |
|
+t*(2467.63998386656941d0 & |
|
+t*(581.879599221457589d0 & |
|
+t*(50.0198090806651216d0 & |
|
+t))))) |
|
elseif(mc.gt.0.161052d0) then |
|
t=7.07638962601280827d0*mc-1.13966670204861480d0 |
|
ceib=(2359.14823394150129d0 & |
|
+t*(3983.28520266051676d0 & |
|
+t*(2254.30785457761760d0 & |
|
+t*(492.601686517364701d0 & |
|
+t*(35.2259786264917876d0 & |
|
+t*0.396605124984359783d0)))))/ & |
|
(2563.95563932625156d0 & |
|
+t*(4450.19076667898892d0 & |
|
+t*(2633.23323959119935d0 & |
|
+t*(622.983787815718489d0 & |
|
+t*(52.6711647124832948d0 & |
|
+t))))) |
|
elseif(mc.gt.0.083522d0) then |
|
t=12.8982329420869341d0*mc-1.07728621178898491d0 |
|
ceib=(2464.65334987833736d0 & |
|
+t*(4333.38639187691528d0 & |
|
+t*(2541.68516994216007d0 & |
|
+t*(571.53606797524881d0 & |
|
+t*(41.5832527504007778d0 & |
|
+t*0.465975784547025267d0)))))/ & |
|
(2600.66956117247726d0 & |
|
+t*(4661.64381841490914d0 & |
|
+t*(2823.69445052534842d0 & |
|
+t*(674.25435972414302d0 & |
|
+t*(56.136001230010910d0 & |
|
+t))))) |
|
elseif(mc.gt.0.041966d0) then |
|
t=24.0639137549331023d0*mc-1.00986620463952257d0 |
|
ceib=(2509.86724450741259d0 & |
|
+t*(4631.12336462339975d0 & |
|
+t*(2835.27071287535469d0 & |
|
+t*(659.86172161727281d0 & |
|
+t*(48.9701196718008345d0 & |
|
+t*0.54158304771955794d0)))))/ & |
|
(2594.15983397593723d0 & |
|
+t*(4848.17491604384532d0 & |
|
+t*(3034.20118545214106d0 & |
|
+t*(737.15143838356850d0 & |
|
+t*(60.652838995496991d0 & |
|
+t))))) |
|
elseif(mc.gt.0.020313d0) then |
|
t=46.1829769546944996d0*mc-0.938114810880709371d0 |
|
ceib=(2480.58307884128017d0 & |
|
+t*(4845.57861173250699d0 & |
|
+t*(3122.00900554841322d0 & |
|
+t*(757.31633816400643d0 & |
|
+t*(57.541132641218839d0 & |
|
+t*0.62119950515996627d0)))))/ & |
|
(2528.85218300581396d0 & |
|
+t*(4979.31783250484768d0 & |
|
+t*(3253.86151324157460d0 & |
|
+t*(812.40556572486862d0 & |
|
+t*(66.496093157522450d0 & |
|
+t))))) |
|
elseif(mc.gt.0.009408d0) then |
|
t=91.7010545621274645d0*mc-0.862723521320495186d0 |
|
ceib=(2365.25385348859592d0 & |
|
+t*(4939.53925884558687d0 & |
|
+t*(3381.09304915246175d0 & |
|
+t*(862.16657576129841d0 & |
|
+t*(67.442026950538221d0 & |
|
+t*0.70143698925710129d0)))))/ & |
|
(2390.48737882063755d0 & |
|
+t*(5015.4675579215077d0 & |
|
+t*(3462.34808443022907d0 & |
|
+t*(898.99542983710459d0 & |
|
+t*(73.934680452209164d0 & |
|
+t))))) |
|
elseif(mc.gt.0.004136d0) then |
|
t=189.681335356600910d0*mc-0.784522003034901366d0 |
|
ceib=(2160.82916040868119d0 & |
|
+t*(4877.14832623847052d0 & |
|
+t*(3584.53058926175721d0 & |
|
+t*(970.53716686804832d0 & |
|
+t*(78.769178005879162d0 & |
|
+t*0.77797110431753920d0)))))/ & |
|
(2172.70451405048305d0 & |
|
+t*(4916.35263668839769d0 & |
|
+t*(3630.52345460629336d0 & |
|
+t*(993.36676027886685d0 & |
|
+t*(83.173163222639080d0 & |
|
+t))))) |
|
elseif(mc.gt.0.d0) then |
|
t=1.d0-106.292517006802721d0*mc |
|
ceib=mc*log(mc*0.0625d0) & |
|
*(6607.46457640413908d0 & |
|
+t*(19.0287633783211078d0 & |
|
-t*0.00625368946932704460d0))/ & |
|
(26150.3443630974309d0 & |
|
+t*(354.603981274536040d0 & |
|
+t)) & |
|
+(26251.5678902584870d0 & |
|
+t*(168.788023807915689d0 & |
|
+t*0.352150236262724288d0))/ & |
|
(26065.7912239203873d0 & |
|
+t*(353.916840382280456d0 & |
|
+t)) |
|
else |
|
write(*,"(a20,1pe15.7)") "(ceib) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(ceib) mc,B=",mc,ceib |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function ceid(mc) |
|
! |
|
! Double precision minimax rational approximation of D(m): |
|
! the second associate complete elliptic integral of the second kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(ceid) negative m: mc=",mc |
|
elseif(mc.gt.0.599909d0) then |
|
t=2.49943137936119533d0*mc-1.49943137936119533d0 |
|
ceid=(1593.39813781813498d0 & |
|
+t*(2233.25576544961714d0 & |
|
+t*(1058.56241259843217d0 & |
|
+t*(195.247394601357872d0 & |
|
+t*(11.7584241242587571d0 & |
|
+t*0.101486443490307517d0)))))/ & |
|
(1685.47865546030468d0 & |
|
+t*(2756.20968383181114d0 & |
|
+t*(1604.88100543517015d0 & |
|
+t*(397.504162950935944d0 & |
|
+t*(38.6743012128666717d0 & |
|
+t))))) |
|
elseif(mc.gt.0.359180d0) then |
|
t=4.15404874360795750d0*mc-1.49205122772910617d0 |
|
ceid=(1967.01442513777287d0 & |
|
+t*(2779.87604145516343d0 & |
|
+t*(1329.30058268219177d0 & |
|
+t*(247.475085945854673d0 & |
|
+t*(15.0447805948342760d0 & |
|
+t*0.130547566005491628d0)))))/ & |
|
(1749.70634057327467d0 & |
|
+t*(2853.92630369567765d0 & |
|
+t*(1654.40804288486242d0 & |
|
+t*(406.925098588378587d0 & |
|
+t*(39.1895256017535337d0 & |
|
+t))))) |
|
elseif(mc.gt.0.214574d0) then |
|
t=6.91534237860116454d0*mc-1.48385267554596628d0 |
|
ceid=(2409.64196912091452d0 & |
|
+t*(3436.40744503228691d0 & |
|
+t*(1659.30176823041376d0 & |
|
+t*(312.186468430688790d0 & |
|
+t*(19.1942111405094383d0 & |
|
+t*0.167847673021897479d0)))))/ & |
|
(1824.89205701262525d0 & |
|
+t*(2971.02216287936566d0 & |
|
+t*(1715.38574780156913d0 & |
|
+t*(418.929791715319490d0 & |
|
+t*(39.8798253173462218d0 & |
|
+t))))) |
|
elseif(mc.gt.0.127875d0) then |
|
t=11.5341584101316047d0*mc-1.47493050669557896d0 |
|
ceid=(2926.81143179637839d0 & |
|
+t*(4214.52119721241319d0 & |
|
+t*(2056.45624281065334d0 & |
|
+t*(391.420514384925370d0 & |
|
+t*(24.3811986813439843d0 & |
|
+t*0.215574280659075512d0)))))/ & |
|
(1910.33091918583314d0 & |
|
+t*(3107.04531802441481d0 & |
|
+t*(1787.99942542734799d0 & |
|
+t*(433.673494280825971d0 & |
|
+t*(40.7663012893484449d0 & |
|
+t))))) |
|
elseif(mc.gt.0.076007d0) then |
|
t=19.2797100331611013d0*mc-1.46539292049047582d0 |
|
ceid=(3520.63614251102960d0 & |
|
+t*(5121.2842239226937d0 & |
|
+t*(2526.67111759550923d0 & |
|
+t*(486.926821696342529d0 & |
|
+t*(30.7739877519417978d0 & |
|
+t*0.276315678908126399d0)))))/ & |
|
(2003.81997889501324d0 & |
|
+t*(3259.09205279874214d0 & |
|
+t*(1871.05914195570669d0 & |
|
+t*(451.007555352632053d0 & |
|
+t*(41.8489850490387023d0 & |
|
+t))))) |
|
elseif(mc.gt.0.045052d0) then |
|
t=32.3049588111775157d0*mc-1.45540300436116944d0 |
|
ceid=(4188.00087087025347d0 & |
|
+t*(6156.0080960857764d0 & |
|
+t*(3072.05695847158556d0 & |
|
+t*(599.76666155374012d0 & |
|
+t*(38.5070211470790031d0 & |
|
+t*0.352955925261363680d0)))))/ & |
|
(2101.60113938424690d0 & |
|
+t*(3421.55151253792527d0 & |
|
+t*(1961.76794074710108d0 & |
|
+t*(470.407158843118117d0 & |
|
+t*(43.0997999502743622d0 & |
|
+t))))) |
|
elseif(mc.gt.0.026626d0) then |
|
t=54.2711386084880061d0*mc-1.44502333658960165d0 |
|
ceid=(4916.74442376570733d0 & |
|
+t*(7304.6632479558695d0 & |
|
+t*(3688.12811638360551d0 & |
|
+t*(729.75841970840314d0 & |
|
+t*(47.6447145147811350d0 & |
|
+t*0.448422756936257635d0)))))/ & |
|
(2197.49982676612397d0 & |
|
+t*(3584.94502590860852d0 & |
|
+t*(2055.19657857622715d0 & |
|
+t*(490.880160668822953d0 & |
|
+t*(44.4576261146308645d0 & |
|
+t))))) |
|
elseif(mc.gt.0.015689d0) then |
|
t=91.4327512114839536d0*mc-1.43448843375697175d0 |
|
ceid=(5688.7542903989517d0 & |
|
+t*(8542.6096475195826d0 & |
|
+t*(4364.21513060078954d0 & |
|
+t*(875.35992968472914d0 & |
|
+t*(58.159468141567195d0 & |
|
+t*0.56528145509695951d0)))))/ & |
|
(2285.44062680812883d0 & |
|
+t*(3739.30422133833258d0 & |
|
+t*(2145.80779422696555d0 & |
|
+t*(511.23253971875808d0 & |
|
+t*(45.8427480379028781d0 & |
|
+t))))) |
|
elseif(mc.gt.0.009216d0) then |
|
t=154.487872701992894d0*mc-1.42376023482156651d0 |
|
ceid=(6475.3392225234969d0 & |
|
+t*(9829.1138694605662d0 & |
|
+t*(5081.2997108708577d0 & |
|
+t*(1033.32687775311981d0 & |
|
+t*(69.910123337464043d0 & |
|
+t*0.70526087421186325d0)))))/ & |
|
(2357.74885505777295d0 & |
|
+t*(3872.32565152553360d0 & |
|
+t*(2226.89527217032394d0 & |
|
+t*(530.03943432061149d0 & |
|
+t*(47.1609071069631012d0 & |
|
+t))))) |
|
elseif(mc.gt.0.d0) then |
|
t=1.d0-108.506944444444444d0*mc |
|
ceid=-log(mc*0.0625d0) & |
|
*(6.2904323649908115d6 & |
|
+t*(58565.284164780476d0 & |
|
+t*(131.176674599188545d0 & |
|
+t*0.0426826410911220304d0)))/ & |
|
(1.24937550257219890d7 & |
|
+t*(203580.534005225410d0 & |
|
+t*(921.17729845011868d0 & |
|
+t))) & |
|
-(27356.1090344387530d0 & |
|
+t*(107.767403612304371d0 & |
|
-t*0.0827769227048233593d0))/ & |
|
(27104.0854889805978d0 & |
|
+t*(358.708172147752755d0 & |
|
+t)) |
|
else |
|
write(*,"(a20,1pe15.7)") "(ceid) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(ceid) mc,D=",mc,ceid |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function sceik(mc) |
|
! |
|
! Single precision minimax rational approximation of K(m): |
|
! the complete elliptic integral of the first kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(sceik) negative m: mc=",mc |
|
elseif(mc.gt.0.265401d0) then |
|
t=1.36128691d0*mc-0.361286906d0 |
|
sceik=(4.52449784d0 & |
|
+t*(13.0783131d0 & |
|
+t*(8.08551381d0 & |
|
+t*0.613554928d0)))/ & |
|
(2.12438551d0 & |
|
+t*(7.37227783d0 & |
|
+t*(6.24763233d0 & |
|
+t))) |
|
elseif(mc.gt.0.068497d0) then |
|
t=5.07861699d0*mc-0.347870028d0 |
|
sceik=(6.02468129d0 & |
|
+t*(18.2518070d0 & |
|
+t*(11.8544683d0 & |
|
+t*0.949905745d0)))/ & |
|
(2.18489368d0 & |
|
+t*(7.70591189d0 & |
|
+t*(6.51975355d0 & |
|
+t))) |
|
elseif(mc.gt.0.017146d0) then |
|
t=19.4738175d0*mc-0.333898074d0 |
|
sceik=(7.50380458d0 & |
|
+t*(23.8982304d0 & |
|
+t*(16.3931558d0 & |
|
+t*1.41140436d0)))/ & |
|
(2.18785514d0 & |
|
+t*(7.90708625d0 & |
|
+t*(6.75017953d0 & |
|
+t))) |
|
elseif(mc.gt.0.d0) then |
|
t=1.d0-58.3226408d0*mc |
|
sceik=-log(mc*0.0625d0) & |
|
*(0.502164162d0 & |
|
+t*(-0.00218520030d0 & |
|
+t*0.0000210458447d0)) & |
|
+mc*(-0.252848894d0 & |
|
+t*(0.00288501596d0 & |
|
-t*0.0000361375649d0)) |
|
else |
|
write(*,"(a20,1pe15.7)") "(sceik) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(sceik) mc,K=",mc,sceik |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function sceie(mc) |
|
! |
|
! Single precision minimax rational approximation of E(m): |
|
! the complete elliptic integral of the second kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(sceie) negative m: mc=",mc |
|
elseif(mc.gt.0.208783d0) then |
|
t=1.26387578d0*mc-0.263875776d0 |
|
sceie=(12.7011016d0 & |
|
+t*(42.7873202d0 & |
|
+t*(34.3113605d0 & |
|
+t*5.70390793d0)))/ & |
|
(10.7238982d0 & |
|
+t*(31.3579945d0 & |
|
+t*(17.7176485d0 & |
|
+t))) |
|
elseif(mc.gt.0.031393d0) then |
|
t=5.63729635d0*mc-0.176971644d0 |
|
sceie=(7.48427789d0 & |
|
+t*(29.6203092d0 & |
|
+t*(23.0865151d0 & |
|
+t*2.92339720d0)))/ & |
|
(7.18634090d0 & |
|
+t*(27.1192312d0 & |
|
+t*(17.9837855d0 & |
|
+t))) |
|
elseif(mc.gt.0.002434d0) then |
|
t=34.5315791d0*mc-0.0840498636d0 |
|
sceie=(3.40850957d0 & |
|
+t*(21.1170648d0 & |
|
+t*(20.0060895d0 & |
|
+t*1.57849869d0)))/ & |
|
(3.39241442d0 & |
|
+t*(20.8511381d0 & |
|
+t*(19.0310444d0 & |
|
+t))) |
|
elseif(mc.gt.0.d0) then |
|
t=1.d0-410.846343d0*mc |
|
sceie=-mc*log(mc*0.0625d0) & |
|
*(0.250228492d0 & |
|
-t*0.000228535208d0) & |
|
+(0.999390295d0 & |
|
+t*(0.000610911723d0 & |
|
-t*1.20643348d-6)) |
|
else |
|
write(*,"(a20,1pe15.7)") "(sceie) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(sceie) mc,K=",mc,sceie |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function sceib(mc) |
|
! |
|
! Single precision minimax rational approximation of E(m): |
|
! the first associate complete elliptic integral of the second kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(sceib) negative m: mc=",mc |
|
elseif(mc.gt.0.189887d0) then |
|
t=1.23439570d0*mc-0.234395695d0 |
|
sceib=(2.17259677d0 & |
|
+t*(7.49197895d0 & |
|
+t*(5.17252751d0 & |
|
+t*0.410720891d0)))/ & |
|
(2.38320592d0 & |
|
+t*(8.81630300d0 & |
|
+t*(7.21462256d0 & |
|
+t))) |
|
elseif(mc.gt.0.025580d0) then |
|
t=6.08616797d0*mc-0.155684177d0 |
|
sceib=(2.16824853d0 & |
|
+t*(9.48517355d0 & |
|
+t*(7.56016795d0 & |
|
+t*0.608774578d0)))/ & |
|
(2.21854578d0 & |
|
+t*(9.94995130d0 & |
|
+t*(8.57542327d0 & |
|
+t))) |
|
elseif(mc.gt.0.001697d0) then |
|
t=41.8707868d0*mc-0.0710547251d0 |
|
sceib=(1.43250252d0 & |
|
+t*(10.1250100d0 & |
|
+t*(10.7354585d0 & |
|
+t*0.824405328d0)))/ & |
|
(1.43625827d0 & |
|
+t*(10.1960125d0 & |
|
+t*(11.0213608d0 & |
|
+t))) |
|
elseif(mc.gt.0.d0) then |
|
t=589.275192d0*mc |
|
sceib=mc*log(mc*0.0625d0) & |
|
*(0.250000000d0 & |
|
+t*(0.000477280426d0 & |
|
+t*8.45890619d-7)) & |
|
+(1.00000000d0 & |
|
+t*(0.00127274774d0 & |
|
+t*2.30087306d-6)) |
|
else |
|
write(*,"(a20,1pe15.7)") "(sceib) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(sceib) mc,K=",mc,sceib |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
real*8 function sceid(mc) |
|
! |
|
! Single precision minimax rational approximation of D(m): |
|
! the second associate complete elliptic integral of the second kind |
|
! |
|
real*8 mc |
|
real*8 t |
|
if(mc.gt.1.d0) then |
|
write(*,"(a20,1pe15.7)") "(sceid) negative m: mc=",mc |
|
elseif(mc.gt.0.278803d0) then |
|
t=1.38658369d0*mc-0.386583693d0 |
|
sceid=(2.29325466d0 & |
|
+t*(5.97114960d0 & |
|
+t*(3.22818778d0 & |
|
+t*0.186239786d0)))/ & |
|
(1.88166457d0 & |
|
+t*(6.42045898d0 & |
|
+t*(5.56782608d0 & |
|
+t))) |
|
elseif(mc.gt.0.075964d0) then |
|
t=4.93001839d0*mc-0.374503917d0 |
|
sceid=(3.65159230d0 & |
|
+t*(9.96877208d0 & |
|
+t*(5.67763554d0 & |
|
+t*0.343806034d0)))/ & |
|
(2.07805961d0 & |
|
+t*(7.06978062d0 & |
|
+t*(5.96868037d0 & |
|
+t))) |
|
elseif(mc.gt.0.020126d0) then |
|
t=17.9089509d0*mc-0.360435546d0 |
|
sceid=(5.34561892d0 & |
|
+t*(15.3896942d0 & |
|
+t*(9.31994615d0 & |
|
+t*0.616190094d0)))/ & |
|
(2.25552058d0 & |
|
+t*(7.75097992d0 & |
|
+t*(6.44810211d0 & |
|
+t))) |
|
elseif(mc.gt.0.d0) then |
|
t=49.6869721d0*mc |
|
sceid=-log(mc*0.0625d0) & |
|
*(0.5d0 & |
|
+t*(0.00754726438d0 & |
|
+t*(0.000142330859d0 & |
|
+t*2.89991893d-6))) & |
|
-(1.d0 & |
|
+t*(0.0201260396d0 & |
|
+t*(0.000389035891d0 & |
|
+t*7.98249490d-6))) |
|
else |
|
write(*,"(a20,1pe15.7)") "(sceid) out of range: mc=",mc |
|
endif |
|
! |
|
!write(*,"(1x,a20,0p5f20.15)") "(sceid) mc,K=",mc,sceid |
|
! |
|
return |
|
end |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
m_c single double |
|
K(m) 0.10000 2.57809209E+00 2.5780921133481733E+00 |
|
K(m) 0.20000 2.25720545E+00 2.2572053268208543E+00 |
|
K(m) 0.30000 2.07536319E+00 2.0753631352924695E+00 |
|
K(m) 0.40000 1.94956778E+00 1.9495677498060255E+00 |
|
K(m) 0.50000 1.85407472E+00 1.8540746773013719E+00 |
|
K(m) 0.60000 1.77751926E+00 1.7775193714912532E+00 |
|
K(m) 0.70000 1.71388946E+00 1.7138894481787912E+00 |
|
K(m) 0.80000 1.65962369E+00 1.6596235986105279E+00 |
|
K(m) 0.90000 1.61244128E+00 1.6124413487202194E+00 |
|
K(m) 1.00000 1.57079642E+00 1.5707963267948966E+00 |
|
|
|
E(m) 0.10000 1.10477478E+00 1.1047747327040733E+00 |
|
E(m) 0.20000 1.17848997E+00 1.1784899243278386E+00 |
|
E(m) 0.30000 1.24167064E+00 1.2416705679458226E+00 |
|
E(m) 0.40000 1.29842796E+00 1.2984280350469131E+00 |
|
E(m) 0.50000 1.35064391E+00 1.3506438810476755E+00 |
|
E(m) 0.60000 1.39939221E+00 1.3993921388974324E+00 |
|
E(m) 0.70000 1.44536302E+00 1.4453630644126652E+00 |
|
E(m) 0.80000 1.48903498E+00 1.4890350580958529E+00 |
|
E(m) 0.90000 1.53075771E+00 1.5307576368977633E+00 |
|
E(m) 1.00000 1.57079623E+00 1.5707963267948963E+00 |
|
|
|
B(m) 0.10000 9.41072750E-01 9.4107280152139550E-01 |
|
B(m) 0.20000 9.08811110E-01 9.0881107370458458E-01 |
|
B(m) 0.30000 8.84373724E-01 8.8437375336868840E-01 |
|
B(m) 0.40000 8.64334940E-01 8.6433489187417145E-01 |
|
B(m) 0.50000 8.47213051E-01 8.4721308479397905E-01 |
|
B(m) 0.60000 8.32201256E-01 8.3220129000670062E-01 |
|
B(m) 0.70000 8.18801537E-01 8.1880150229170512E-01 |
|
B(m) 0.80000 8.06680930E-01 8.0668089603715254E-01 |
|
B(m) 0.90000 7.95604191E-01 7.9560423049565732E-01 |
|
B(m) 1.00000 7.85398210E-01 7.8539816339744839E-01 |
|
|
|
D(m) 0.10000 1.63701922E+00 1.6370193118267773E+00 |
|
D(m) 0.20000 1.34839428E+00 1.3483942531162689E+00 |
|
D(m) 0.30000 1.19098945E+00 1.1909893819237805E+00 |
|
D(m) 0.40000 1.08523285E+00 1.0852328579318544E+00 |
|
D(m) 0.50000 1.00686163E+00 1.0068615925073927E+00 |
|
D(m) 0.60000 9.45318027E-01 9.4531808148455243E-01 |
|
D(m) 0.70000 8.95087948E-01 8.9508794588708585E-01 |
|
D(m) 0.80000 8.52942753E-01 8.5294270257337523E-01 |
|
D(m) 0.90000 8.16837090E-01 8.1683711822456173E-01 |
|
D(m) 1.00000 7.85398211E-01 7.8539816339744850E-01 |
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
