From e08d5c7655b2f3b9c38e2587b09260139ede9dd7 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 4 Jan 2024 16:06:06 -0500 Subject: [PATCH] incorporate elliptic2 functionality directly into package and remove contour dep --- Manifest.toml | 230 +----------- Project.toml | 14 +- src/LoopFieldCalc.jl | 5 +- src/elliptic2/Elliptic2.jl | 122 +++++++ src/elliptic2/slatec.jl | 343 ++++++++++++++++++ test/autodiff.jl | 53 +++ test/data/readme.txt | 8 + test/data/table_16_1.csv | 343 ++++++++++++++++++ test/data/table_17_1.csv | 102 ++++++ test/data/table_17_2.csv | 92 +++++ test/runtests.jl | 2 + test/runtests_elliptic.jl | 720 +++++++++++++++++++++++++++++++++++++ 12 files changed, 1797 insertions(+), 237 deletions(-) create mode 100644 src/elliptic2/Elliptic2.jl create mode 100644 src/elliptic2/slatec.jl create mode 100644 test/autodiff.jl create mode 100644 test/data/readme.txt create mode 100644 test/data/table_16_1.csv create mode 100644 test/data/table_17_1.csv create mode 100644 test/data/table_17_2.csv create mode 100644 test/runtests_elliptic.jl diff --git a/Manifest.toml b/Manifest.toml index 5c5a03b..05ba538 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,239 +2,11 @@ julia_version = "1.10.0-rc3" manifest_format = "2.0" -project_hash = "dd0d912228d3570dd0b442fd72dbc8a9b466421d" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.1" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" - -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+1" - -[[deps.Contour]] -git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781" -uuid = "d38c429a-6771-53c6-b99e-75d170b6e991" -version = "0.6.2" - -[[deps.Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" - -[[deps.DelimitedFiles]] -deps = ["Mmap"] -git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae" -uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" -version = "1.9.1" - -[[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.3" - -[[deps.Downloads]] -deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" - -[[deps.Elliptic2]] -deps = ["DelimitedFiles", "SpecialFunctions"] -git-tree-sha1 = "fc0b6310c912f838661e19e739667d979683569c" -repo-rev = "master" -repo-url = "https://github.com/archermarx/Elliptic2.jl" -uuid = "0b55928f-c416-41db-a823-c7dfe4e8fda4" -version = "0.2.0" - -[[deps.FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" - -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[deps.IrrationalConstants]] -git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" -uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.2.2" - -[[deps.JLLWrappers]] -deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" -uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.5.0" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.4" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.4.0+0" - -[[deps.LibGit2]] -deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[deps.LibGit2_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] -uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.6.4+0" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.0+1" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[deps.LinearAlgebra]] -deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[deps.LogExpFunctions]] -deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.26" - - [deps.LogExpFunctions.extensions] - LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" - LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" - LogExpFunctionsInverseFunctionsExt = "InverseFunctions" - - [deps.LogExpFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" - -[[deps.Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+1" - -[[deps.Mmap]] -uuid = "a63ad114-7e13-5084-954f-fe012c677804" - -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2023.1.10" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" - -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+2" - -[[deps.OpenLibm_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+2" - -[[deps.OpenSpecFun_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" -uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.5+0" - -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.10.0" - -[[deps.Preferences]] -deps = ["TOML"] -git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e" -uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.4.1" +project_hash = "3ec3cda12b6de964321693f8135f0ce728085bf1" [[deps.Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" -[[deps.REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] -uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[deps.Random]] -deps = ["SHA"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[[deps.SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -version = "0.7.0" - -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[deps.Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" - -[[deps.SpecialFunctions]] -deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d" -uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.3.1" - - [deps.SpecialFunctions.extensions] - SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" - - [deps.SpecialFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - -[[deps.TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.3" - -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.0" - -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+1" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.8.0+1" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.52.0+1" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+2" diff --git a/Project.toml b/Project.toml index c7016d8..f26dea1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,21 +1,21 @@ name = "LoopFieldCalc" uuid = "896acac2-5fe1-47fe-8ead-27e3a9bc5c85" authors = ["Thomas Marks and contributors"] -version = "0.3.1" +version = "0.3.3" [deps] -Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" -Elliptic2 = "0b55928f-c416-41db-a823-c7dfe4e8fda4" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" [compat] -julia = "1" -Elliptic2 = ">=0.1" -Contour = "0.6.2" Printf = "<0.0.1, 1" +julia = "1" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" [targets] -test = ["Test"] +test = ["Test", "SpecialFunctions", "DelimitedFiles", "ForwardDiff", "Zygote"] diff --git a/src/LoopFieldCalc.jl b/src/LoopFieldCalc.jl index c591f6f..4375bec 100644 --- a/src/LoopFieldCalc.jl +++ b/src/LoopFieldCalc.jl @@ -1,6 +1,9 @@ module LoopFieldCalc -using Elliptic2, Printf, Contour +include("elliptic2//Elliptic2.jl") + +using Printf +using .Elliptic2 const μ₀ = π * 4e-7 const μ0 = μ₀ diff --git a/src/elliptic2/Elliptic2.jl b/src/elliptic2/Elliptic2.jl new file mode 100644 index 0000000..6c3eb85 --- /dev/null +++ b/src/elliptic2/Elliptic2.jl @@ -0,0 +1,122 @@ +# NOTE: this is a differentiable version of Elliptic.jl +module Elliptic2 +# elliptic integrals of 1st/2nd/3rd kind +export E, F, K, Pi + +# matlab compatible +export ellipj, ellipke + +include("slatec.jl") + +function E(phi, m) + if isnan(phi) || isnan(m) return NaN end + if m < 0. || m > 1. throw(DomainError(m, "argument m not in [0,1]")) end + if abs(phi) > pi/2 + phi2 = phi + pi/2 + return 2*fld(phi2,pi)*E(m) - _E(cos(mod(phi2,pi)), m) + end + _E(sin(phi), m) +end +function _E(sinphi, m) + sinphi2 = sinphi^2 + cosphi2 = 1. - sinphi2 + y = 1. - m*sinphi2 + drf,ierr1 = SLATEC.DRF(cosphi2, y, 1.) + drd,ierr2 = SLATEC.DRD(cosphi2, y, 1.) + if ierr1 == ierr2 == 0 + return sinphi*(drf - m*sinphi2*drd/3) + elseif ierr1 == ierr2 == 2 + # 2 - (1+m)*sinphi2 < tol + return sinphi + end + NaN +end + +""" +`ellipke(m::Real)` +returns `(K(m), E(m))` for scalar `0 ≤ m ≤ 1` +""" +function ellipke(m) + if m < 1. + y = 1. - m + drf,ierr1 = SLATEC.DRF(0., y, 1.) + drd,ierr2 = SLATEC.DRD(0., y, 1.) + @assert ierr1 == 0 && ierr2 == 0 + return drf, drf - m*drd/3 + elseif m == 1. + return Inf, 1. + elseif isnan(m) + return NaN, NaN + else + throw(DomainError(m, "argument m not <= 1")) + end +end + +E(m) = ellipke(m)[2] + +# assumes 0 ≤ m ≤ 1 +function rawF(sinphi, m) + if abs(sinphi) == 1. && m == 1. return sign(sinphi)*Inf end + sinphi2 = sinphi^2 + drf,ierr = SLATEC.DRF(1. - sinphi2, 1. - m*sinphi2, 1.) + @assert ierr == 0 + sinphi*drf +end + +function F(phi, m) + if isnan(phi) || isnan(m) return NaN end + if m < 0. || m > 1. throw(DomainError(m, "argument m not in [0,1]")) end + if abs(phi) > pi/2 + # Abramowitz & Stegun (17.4.3) + phi2 = phi + pi/2 + return 2*fld(phi2,pi)*K(m) - rawF(cos(mod(phi2,pi)), m) + end + rawF(sin(phi), m) +end + +function K(m) + if m < 1. + drf,ierr = SLATEC.DRF(0., 1. - m, 1.) + @assert ierr == 0 + return drf + elseif m == 1. + return Inf + elseif isnan(m) + return NaN + else + throw(DomainError(m, "argument m not <= 1")) + end +end + +function Pi(n, phi, m) + if isnan(n) || isnan(phi) || isnan(m) return NaN end + if m < 0. || m > 1. throw(DomainError(m, "argument m not in [0,1]")) end + sinphi = sin(phi) + sinphi2 = sinphi^2 + cosphi2 = 1. - sinphi2 + y = 1. - m*sinphi2 + drf,ierr1 = SLATEC.DRF(cosphi2, y, 1.) + drj,ierr2 = SLATEC.DRJ(cosphi2, y, 1., 1. - n*sinphi2) + if ierr1 == 0 && ierr2 == 0 + return sinphi*(drf + n*sinphi2*drj/3) + elseif ierr1 == 2 && ierr2 == 2 + # 2 - (1+m)*sinphi2 < tol + return Inf + elseif ierr1 == 0 && ierr2 == 2 + # 1 - n*sinphi2 < tol + return Inf + end + NaN +end +Π = Pi + +function ellipj(u, m, tol) + phi = Jacobi.am(u, m, tol) + s = sin(phi) + c = cos(phi) + d = sqrt(1. - m*s^2) + s, c, d +end +ellipj(u, m) = ellipj(u, m, eps(Float64)) + +end # module diff --git a/src/elliptic2/slatec.jl b/src/elliptic2/slatec.jl new file mode 100644 index 0000000..21f36b4 --- /dev/null +++ b/src/elliptic2/slatec.jl @@ -0,0 +1,343 @@ +module SLATEC + +export DRC, DRD, DRF, DRJ + +const D1MACH1 = floatmin(Float64) +const D1MACH2 = floatmax(Float64) +const D1MACH3 = eps(Float64)/2 +const D1MACH4 = eps(Float64) +const D1MACH5 = log10(2.) + +#***BEGIN PROLOGUE DRF +#***PURPOSE Compute the incomplete or complete elliptic integral of the +# 1st kind. For X, Y, and Z non-negative and at most one of +# them zero, RF(X,Y,Z) = Integral from zero to infinity of +# -1/2 -1/2 -1/2 +# (1/2)(t+X) (t+Y) (t+Z) dt. +# If X, Y or Z is zero, the integral is complete. +#***LIBRARY SLATEC +#***CATEGORY C14 +#***TYPE DOUBLE PRECISION (RF-S, DRF-D) +#***KEYWORDS COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM, +# INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE FIRST KIND, +# TAYLOR SERIES +#***AUTHOR Carlson, B. C. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Notis, E. M. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Pexton, R. L. +# Lawrence Livermore National Laboratory +# Livermore, CA 94550 + +function DRF(X, Y, Z) + + ERRTOL = (4.0*D1MACH3)^(1.0/6.0) + LOLIM = 5.0 * D1MACH1 + UPLIM = D1MACH2/5.0 + C1 = 1.0/24.0 + C2 = 3.0/44.0 + C3 = 1.0/14.0 + + ans = 0.0 + if min(X,Y,Z) < 0.0 + return ans, 1 + end + + if max(X,Y,Z) > UPLIM + return ans, 3 + end + + if min(X+Y,X+Z,Y+Z) < LOLIM + return ans, 2 + end + + XN = X + YN = Y + ZN = Z + MU = 0. + XNDEV = 0. + YNDEV = 0. + ZNDEV = 0. + + while true + MU = (XN+YN+ZN)/3.0 + XNDEV = 2.0 - (MU+XN)/MU + YNDEV = 2.0 - (MU+YN)/MU + ZNDEV = 2.0 - (MU+ZN)/MU + EPSLON = max(abs(XNDEV),abs(YNDEV),abs(ZNDEV)) + if (EPSLON < ERRTOL) break end + XNROOT = sqrt(XN) + YNROOT = sqrt(YN) + ZNROOT = sqrt(ZN) + LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT + XN = (XN+LAMDA)*0.250 + YN = (YN+LAMDA)*0.250 + ZN = (ZN+LAMDA)*0.250 + end + + E2 = XNDEV*YNDEV - ZNDEV*ZNDEV + E3 = XNDEV*YNDEV*ZNDEV + S = 1.0 + (C1*E2-0.10-C2*E3)*E2 + C3*E3 + ans = S/sqrt(MU) + + return ans, 0 +end + +#***BEGIN PROLOGUE DRD +#***PURPOSE Compute the incomplete or complete elliptic integral of +# the 2nd kind. For X and Y nonnegative, X+Y and Z positive, +# DRD(X,Y,Z) = Integral from zero to infinity of +# -1/2 -1/2 -3/2 +# (3/2)(t+X) (t+Y) (t+Z) dt. +# If X or Y is zero, the integral is complete. +#***LIBRARY SLATEC +#***CATEGORY C14 +#***TYPE DOUBLE PRECISION (RD-S, DRD-D) +#***KEYWORDS COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM, +# INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE SECOND KIND, +# TAYLOR SERIES +#***AUTHOR Carlson, B. C. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Notis, E. M. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Pexton, R. L. +# Lawrence Livermore National Laboratory +# Livermore, CA 94550 + +function DRD(X, Y, Z) + + ERRTOL = (D1MACH3/3.0)^(1.0/6.0) + LOLIM = 2.0/(D1MACH2)^(2.0/3.0) + TUPLIM = D1MACH1^(1.0E0/3.0E0) + TUPLIM = (0.10*ERRTOL)^(1.0E0/3.0E0)/TUPLIM + UPLIM = TUPLIM^2.0 + C1 = 3.0/14.0 + C2 = 1.0/6.0 + C3 = 9.0/22.0 + C4 = 3.0/26.0 + + ans = 0.0 + if min(X,Y) < 0.0 + return ans, 1 + end + + if max(X,Y,Z) > UPLIM + return ans, 3 + end + + if min(X+Y,Z) < LOLIM + return ans, 2 + end + + XN = X + YN = Y + ZN = Z + SIGMA = 0.0 + POWER4 = 1.0 + MU = 0. + XNDEV = 0. + YNDEV = 0. + ZNDEV = 0. + + while true + MU = (XN+YN+3.0*ZN)*0.20 + XNDEV = (MU-XN)/MU + YNDEV = (MU-YN)/MU + ZNDEV = (MU-ZN)/MU + EPSLON = max(abs(XNDEV), abs(YNDEV), abs(ZNDEV)) + if (EPSLON < ERRTOL) break end + XNROOT = sqrt(XN) + YNROOT = sqrt(YN) + ZNROOT = sqrt(ZN) + LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT + SIGMA = SIGMA + POWER4/(ZNROOT*(ZN+LAMDA)) + POWER4 = POWER4*0.250 + XN = (XN+LAMDA)*0.250 + YN = (YN+LAMDA)*0.250 + ZN = (ZN+LAMDA)*0.250 + end + + EA = XNDEV*YNDEV + EB = ZNDEV*ZNDEV + EC = EA - EB + ED = EA - 6.0*EB + EF = ED + EC + EC + S1 = ED*(-C1+0.250*C3*ED-1.50*C4*ZNDEV*EF) + S2 = ZNDEV*(C2*EF+ZNDEV*(-C3*EC+ZNDEV*C4*EA)) + ans = 3.0*SIGMA + POWER4*(1.0+S1+S2)/(MU*sqrt(MU)) + + return ans, 0 +end + +#***BEGIN PROLOGUE DRC +#***PURPOSE Calculate a double precision approximation to +# DRC(X,Y) = Integral from zero to infinity of +# -1/2 -1 +# (1/2)(t+X) (t+Y) dt, +# where X is nonnegative and Y is positive. +#***LIBRARY SLATEC +#***CATEGORY C14 +#***TYPE DOUBLE PRECISION (RC-S, DRC-D) +#***KEYWORDS DUPLICATION THEOREM, ELEMENTARY FUNCTIONS, +# ELLIPTIC INTEGRAL, TAYLOR SERIES +#***AUTHOR Carlson, B. C. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Notis, E. M. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Pexton, R. L. +# Lawrence Livermore National Laboratory +# Livermore, CA 94550 + +function DRC(X, Y) + + ERRTOL = (D1MACH3/16.0)^(1.0/6.0) + LOLIM = 5.0 * D1MACH1 + UPLIM = D1MACH2 / 5.0 + C1 = 1.0/7.0 + C2 = 9.0/22.0 + + ans = 0.0 + if X < 0.0 || Y <= 0.0 + return ans, 1 + end + + if max(X,Y) > UPLIM + return ans, 3 + end + + if X+Y < LOLIM + return ans, 2 + end + + XN = X + YN = Y + MU = 0. + SN = 0. + + while true + MU = (XN+YN+YN)/3.0 + SN = (YN+MU)/MU - 2.0 + if abs(SN) < ERRTOL break end + LAMDA = 2.0*sqrt(XN)*sqrt(YN) + YN + XN = (XN+LAMDA)*0.250 + YN = (YN+LAMDA)*0.250 + end + + S = SN*SN*(0.30+SN*(C1+SN*(0.3750+SN*C2))) + ans = (1.0+S)/sqrt(MU) + + return ans,0 +end + +#***BEGIN PROLOGUE DRJ +#***PURPOSE Compute the incomplete or complete (X or Y or Z is zero) +# elliptic integral of the 3rd kind. For X, Y, and Z non- +# negative, at most one of them zero, and P positive, +# RJ(X,Y,Z,P) = Integral from zero to infinity of +# -1/2 -1/2 -1/2 -1 +# (3/2)(t+X) (t+Y) (t+Z) (t+P) dt. +#***LIBRARY SLATEC +#***CATEGORY C14 +#***TYPE DOUBLE PRECISION (RJ-S, DRJ-D) +#***KEYWORDS COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM, +# INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE THIRD KIND, +# TAYLOR SERIES +#***AUTHOR Carlson, B. C. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Notis, E. M. +# Ames Laboratory-DOE +# Iowa State University +# Ames, IA 50011 +# Pexton, R. L. +# Lawrence Livermore National Laboratory +# Livermore, CA 94550 + +function DRJ(X, Y, Z, P) + + ERRTOL = (D1MACH3/3.0)^(1.0/6.0) + LOLIM = (5.0 * D1MACH1)^(1.0/3.0) + UPLIM = 0.30*( D1MACH2 / 5.0)^(1.0/3.0) + + C1 = 3.0/14.0 + C2 = 1.0/3.0 + C3 = 3.0/22.0 + C4 = 3.0/26.0 + + ans = 0.0 + if min(X,Y,Z) < 0.0 + return ans, 1 + end + + if max(X,Y,Z,P) > UPLIM + return ans, 3 + end + + if min(X+Y,X+Z,Y+Z,P) < LOLIM + return ans, 2 + end + + IER = 0 + XN = X + YN = Y + ZN = Z + PN = P + SIGMA = 0.0 + POWER4 = 1.0 + MU = 0. + XNDEV = 0. + YNDEV = 0. + ZNDEV = 0. + PNDEV = 0. + + while true + MU = (XN+YN+ZN+PN+PN)*0.20 + XNDEV = (MU-XN)/MU + YNDEV = (MU-YN)/MU + ZNDEV = (MU-ZN)/MU + PNDEV = (MU-PN)/MU + EPSLON = max(abs(XNDEV), abs(YNDEV), abs(ZNDEV), abs(PNDEV)) + if EPSLON < ERRTOL break end + XNROOT = sqrt(XN) + YNROOT = sqrt(YN) + ZNROOT = sqrt(ZN) + LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT + ALFA = PN*(XNROOT+YNROOT+ZNROOT) + XNROOT*YNROOT*ZNROOT + ALFA = ALFA*ALFA + BETA = PN*(PN+LAMDA)*(PN+LAMDA) + drc,IER = DRC(ALFA,BETA) + SIGMA = SIGMA + POWER4*drc + POWER4 = POWER4*0.250 + XN = (XN+LAMDA)*0.250 + YN = (YN+LAMDA)*0.250 + ZN = (ZN+LAMDA)*0.250 + PN = (PN+LAMDA)*0.250 + end + + EA = XNDEV*(YNDEV+ZNDEV) + YNDEV*ZNDEV + EB = XNDEV*YNDEV*ZNDEV + EC = PNDEV*PNDEV + E2 = EA - 3.0*EC + E3 = EB + 2.0*PNDEV*(EA-EC) + S1 = 1.0 + E2*(-C1+0.750*C3*E2-1.50*C4*E3) + S2 = EB*(0.50*C2+PNDEV*(-C3-C3+PNDEV*C4)) + S3 = PNDEV*EA*(C2-PNDEV*C3) - C2*PNDEV*EC + ans = 3.0*SIGMA + POWER4*(S1+S2+S3)/(MU* sqrt(MU)) + + return ans, IER +end + +end # module diff --git a/test/autodiff.jl b/test/autodiff.jl new file mode 100644 index 0000000..5c5a186 --- /dev/null +++ b/test/autodiff.jl @@ -0,0 +1,53 @@ +using Zygote +using ForwardDiff +using SpecialFunctions + +@testset "Zygote and ForwardDiff" begin + + num_trials = 1 + ks = rand(num_trials) + ϕs = rand(num_trials) .* 2π + ns = rand(num_trials) + + @show ks, ϕs, ns + + # Test several known derivative identities across a wide range of values of ϕ and m + # in order to verify that derivatives work correctly + for (k, ϕ, n) in zip(ks, ϕs, ns) + m = k^2 + dk_dm = 0.5 / k + + # I. Tests for complete integrals, from https://en.wikipedia.org/wiki/Elliptic_integral + # 1. K'(m) == K'(k) * dk/dm == E(k) / (k * (1 - k^2)) - K(k)/k + @test Zygote.gradient(K, m)[1] ≈ (E(m) / (k * (1 - k^2)) - K(m) / k) * dk_dm + @test ForwardDiff.derivative(K, m) ≈ (E(m) / (k * (1 - k^2)) - K(m) / k) * dk_dm + + # 2. E'(m) == E'(k) * dk/dm == (E(k) - K(k))/k + @test Zygote.gradient(E, m)[1] ≈ (E(m) - K(m))/k * dk_dm + @test ForwardDiff.derivative(E, m) ≈ (E(m) - K(m))/k * dk_dm + + # II. Tests for incomplete integrals, from https://functions.wolfram.com/EllipticIntegrals/EllipticF/introductions/IncompleteEllipticIntegrals/ShowAll.html + # 3. ∂ϕ(F(ϕ, m)) == 1 / √(1 - m*sin(ϕ)^2) + @test Zygote.gradient(ϕ -> F(ϕ, m), ϕ)[1] ≈ 1 / √(1 - m*sin(ϕ)^2) + @test ForwardDiff.derivative(ϕ -> F(ϕ, m), ϕ) ≈ 1 / √(1 - m*sin(ϕ)^2) + + # 4. ∂m(F(ϕ, m)) == E(ϕ, m) / (2 * m * (1 - m)) - F(ϕ, m) / 2m - sin(2ϕ) / (4 * (1-m) * √(1 - m * sin(ϕ)^2)) + @test Zygote.gradient(m -> F(ϕ, m), m)[1] ≈ + E(ϕ, m) / (2 * m * (1 - m)) - + F(ϕ, m) / 2 / m - + sin(2*ϕ) / (4 * (1 - m) * √(1 - m * sin(ϕ)^2)) + @test ForwardDiff.derivative(m -> F(ϕ, m), m) ≈ + E(ϕ, m) / (2 * m * (1 - m)) - + F(ϕ, m) / 2 / m - + sin(2*ϕ) / (4 * (1 - m) * √(1 - m * sin(ϕ)^2)) + + # 5. ∂ϕ(E(ϕ, m)) == √(1 - m * sin(ϕ)^2) + @test Zygote.gradient(ϕ -> E(ϕ, m), ϕ)[1] ≈ √(1 - m * sin(ϕ)^2) + @test ForwardDiff.derivative(ϕ -> E(ϕ, m), ϕ) ≈ √(1 - m * sin(ϕ)^2) + + # 6. ∂m(E(ϕ, m)) == (E(ϕ, m) - F(ϕ, m)) / 2m + @test Zygote.gradient(m -> E(ϕ, m), m)[1] ≈ (E(ϕ, m) - F(ϕ, m)) / 2m + @test ForwardDiff.derivative(m -> E(ϕ, m), m) ≈ (E(ϕ, m) - F(ϕ, m)) / 2m + + end +end diff --git a/test/data/readme.txt b/test/data/readme.txt new file mode 100644 index 0000000..0b34f3f --- /dev/null +++ b/test/data/readme.txt @@ -0,0 +1,8 @@ +data and tables (under fair use of unlcensed work) from +Abramowitz, Milton & Stegun, Irene A. +Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables +https://digital.library.unt.edu/ark:/67531/metadc40301/ + +table 16.1: pages 582-583 +table 17.1: pages 608-609 +table 17.2: pages 610-611 diff --git a/test/data/table_16_1.csv b/test/data/table_16_1.csv new file mode 100644 index 0000000..a95c186 --- /dev/null +++ b/test/data/table_16_1.csv @@ -0,0 +1,343 @@ +α,ϵ,θs,θn + 0, 0,0.000000000,1.0000000000 + 0, 5,0.087155743,1.0000000000 + 0,10,0.173648178,1.0000000000 + 0,15,0.258819045,1.0000000000 + 0,20,0.342020143,1.0000000000 + 0,25,0.422618262,1.0000000000 + 0,30,0.500000000,1.0000000000 + 0,35,0.573576436,1.0000000000 + 0,40,0.642787610,1.0000000000 + 0,45,0.707106781,1.0000000000 + 0,50,0.766044443,1.0000000000 + 0,55,0.819152044,1.0000000000 + 0,60,0.866025404,1.0000000000 + 0,65,0.906307787,1.0000000000 + 0,70,0.939692621,1.0000000000 + 0,75,0.965925826,1.0000000000 + 0,80,0.984807753,1.0000000000 + 0,85,0.996194698,1.0000000000 + 0,90,1.000000000,1.0000000000 + 5, 0, 0.000000000, 1.0000000000 + 5, 5, 0.087321966, 1.0000144942 + 5,10, 0.173979362, 1.0000575362 + 5,15, 0.259312677, 1.0001278184 + 5,20, 0.342672476, 1.0002232051 + 5,25, 0.423424343, 1.0003407982 + 5,30, 0.500953708, 1.0004770246 + 5,35, 0.574670526, 1.0006277451 + 5,40, 0.644013768, 1.0007883803 + 5,45, 0.708455688, 1.0009540492 + 5,50, 0.767505843, 1.0011197181 + 5,55, 0.820714821, 1.0012803532 + 5,60, 0.867677668, 1.0014310738 + 5,65, 0.908036964, 1.0015673002 + 5,70, 0.941485546, 1.0016848932 + 5,75, 0.967768848, 1.0017802800 + 5,80, 0.986686836, 1.0018505621 + 5,85, 0.998095528, 1.0018936042 + 5,90, 1.001908098, 1.0019080984 +10, 0, 0.000000000, 1.0000000000 +10, 5, 0.087824152, 1.0000583670 +10,10, 0.174979967, 1.0002316945 +10,15, 0.260804191, 1.0005147160 +10,20, 0.344643695, 1.0008988322 +10,25, 0.425860446, 1.0013723717 +10,30, 0.503836358, 1.0019209464 +10,35, 0.577977994, 1.0025278880 +10,40, 0.647721085, 1.0031747551 +10,45, 0.712534820, 1.0038418928 +10,50, 0.771925893, 1.0045090305 +10,55, 0.825442256, 1.0051558975 +10,60, 0.872676562, 1.0057628392 +10,65, 0.913269273, 1.0063114139 +10,70, 0.946911395, 1.0067849535 +10,75, 0.973346839, 1.0071690696 +10,80, 0.992374367, 1.0074520912 +10,85, 1.003849133, 1.0076254187 +10,90, 1.007683786, 1.0076837857 +15, 0, 0.000000000, 1.0000000000 +15, 5, 0.088673070, 1.0001328199 +15,10, 0.176671584, 1.0005272438 +15,15, 0.263326099, 1.0011712875 +15,20, 0.347977361, 1.0020453820 +15,25, 0.429981306, 1.0031229684 +15,30, 0.508713952, 1.0043713049 +15,35, 0.583576134, 1.0057524612 +15,40, 0.653998067, 1.0072244718 +15,45, 0.719443681, 1.0087426104 +15,50, 0.779414712, 1.0102607491 +15,55, 0.833454505, 1.0117327599 +15,60, 0.881151505, 1.0131139167 +15,65, 0.922142410, 1.0143622536 +15,70, 0.956114956, 1.0154398405 +15,75, 0.982810311, 1.0163139354 +15,80, 1.002025068, 1.0169579795 +15,85, 1.013612807, 1.0173524037 +15,90, 1.017485224, 1.0174852237 +20, 0, 0.000000000, 1.0000000000 +20, 5, 0.089887414, 1.0002399605 +20,10, 0.179091708, 1.0009525510 +20,15, 0.266934892, 1.0021161200 +20,20, 0.352749211, 1.0036953131 +20,25, 0.435882163, 1.0056421475 +20,30, 0.515701435, 1.0078974700 +20,35, 0.591599683, 1.0103927539 +20,40, 0.662999145, 1.0130521815 +20,45, 0.729356053, 1.0157949474 +20,50, 0.790164790, 1.0185377143 +20,55, 0.844961783, 1.0211971444 +20,60, 0.893329083, 1.0236924323 +20,65, 0.934897610, 1.0259477596 +20,70, 0.969350025, 1.0278945992 +20,75, 0.996423213, 1.0294737972 +20,80, 1.015910350, 1.0306373701 +20,85, 1.027662527, 1.0313499632 +20,90, 1.031589925, 1.0315899246 +25, 0, 0.000000000, 1.0000000000 +25, 5, 0.091495034, 1.0003829783 +25,10, 0.182296223, 1.0015202770 +25,15, 0.271714833, 1.0033773404 +25,20, 0.359072325, 1.0058977438 +25,25, 0.443705382, 1.0090049074 +25,30, 0.524970857, 1.0126044231 +25,35, 0.602250597, 1.0165869227 +25,40, 0.674956130, 1.0208314013 +25,45, 0.742533161, 1.0252088930 +25,50, 0.804465863, 1.0295863905 +25,55, 0.860280899, 1.0338308852 +25,60, 0.909551166, 1.0378134098 +25,65, 0.951899199, 1.0414129561 +25,70, 0.987000216, 1.0445201522 +25,75, 1.014584761, 1.0470405862 +25,80, 1.034440908, 1.0488976746 +25,85, 1.046416011, 1.0500349895 +25,90, 1.050417974, 1.0504179735 +30, 0, 0.000000000, 1.0000000000 +30, 5, 0.093534894, 1.0005664294 +30,10, 0.186363367, 1.0022485079 +30,15, 0.277784006, 1.0049951300 +30,20, 0.367105393, 1.0087228461 +30,25, 0.453651078, 1.0133183978 +30,30, 0.536764494, 1.0186421583 +30,35, 0.615813814, 1.0245323743 +30,40, 0.690196708, 1.0308100797 +30,45, 0.759344980, 1.0372845330 +30,50, 0.822729031, 1.0437590125 +30,55, 0.879862121, 1.0500367930 +30,60, 0.930304365, 1.0559271242 +30,65, 0.973666431, 1.0612510260 +30,70, 1.009612870, 1.0658467280 +30,75, 1.037865044, 1.0695745853 +30,80, 1.058203585, 1.0723213226 +30,85, 1.070470366, 1.0740034764 +30,90, 1.074569932, 1.0745699318 +35, 0, 0.000000000, 1.0000000000 +35, 5, 0.096060073, 1.0007966833 +35,10, 0.191399811, 1.0031625308 +35,15, 0.285303629, 1.0070256701 +35,20, 0.377065455, 1.0122687413 +35,25, 0.465993521, 1.0187324599 +35,30, 0.551415176, 1.0262204548 +35,35, 0.632681725, 1.0345052308 +35,40, 0.709173264, 1.0433350787 +35,45, 0.780303503, 1.0524417208 +35,50, 0.845524503, 1.0615484606 +35,55, 0.904331298, 1.0703785902 +35,60, 0.956266326, 1.0786637978 +35,65, 1.000923589, 1.0861523221 +35,70, 1.037952481, 1.0926166042 +35,75, 1.067061179, 1.0978602047 +35,80, 1.088019556, 1.1017237756 +35,85, 1.100661511, 1.1040899048 +35,90, 1.104886686, 1.1048866859 +40, 0, 0.000000000, 1.0000000000 +40, 5, 0.099142353, 1.0010826253 +40,10, 0.197549961, 1.0042976203 +40,15, 0.294492321, 1.0095473402 +40,20, 0.389247478, 1.0166723379 +40,25, 0.481106437, 1.0254562012 +40,30, 0.569377735, 1.0356321191 +40,35, 0.653392178, 1.0468909786 +40,40, 0.732507761, 1.0588907481 +40,45, 0.806114729, 1.0712668617 +40,50, 0.873640739, 1.0836432917 +40,55, 0.934556042, 1.0956439724 +40,60, 0.988378598, 1.1069042279 +40,65, 1.034678996, 1.1170818582 +40,70, 1.073085074, 1.1258675438 +40,75, 1.103286100, 1.1329942539 +40,80, 1.125036391, 1.1382453698 +40,85, 1.138158265, 1.1414612760 +40,90, 1.142544218, 1.1425442177 +45, 0, 0.000000000, 1.0000000000 +45, 5, 0.102879331, 1.0014367802 +45,10, 0.205010420, 1.0057035065 +45,15, 0.305648349, 1.0126706562 +45,20, 0.404054995, 1.0221267193 +45,25, 0.499502749, 1.0337846028 +45,30, 0.591278602, 1.0472903271 +45,35, 0.678688658, 1.0622337524 +45,40, 0.761063101, 1.0781610137 +45,45, 0.837761607, 1.0945882886 +45,50, 0.908179128, 1.1110164844 +45,55, 0.971751955, 1.1269463970 +45,60, 1.027963895, 1.1418938846 +45,65, 1.076352410, 1.1554045920 +45,70, 1.116514503, 1.1670677783 +45,75, 1.148112152, 1.1765288244 +45,80, 1.170877087, 1.1835000363 +45,85, 1.184614727, 1.1877694140 +45,90, 1.189207115, 1.1892071150 +50, 0, 0.000000000, 1.0000000000 +50, 5, 0.107405819, 1.0018771775 +50,10, 0.214053194, 1.0074517850 +50,15, 0.319185434, 1.0165547635 +50,20, 0.422049614, 1.0289100179 +50,25, 0.521899092, 1.0441427466 +50,30, 0.617996720, 1.0617907561 +50,35, 0.709618904, 1.0813184270 +50,40, 0.796060581, 1.1021329153 +50,45, 0.876641114, 1.1236021058 +50,50, 0.950711025, 1.1450737802 +50,55, 1.017659399, 1.1658954205 +50,60, 1.076921759, 1.1854340490 +50,65, 1.127988100, 1.2030954999 +50,70, 1.170410792, 1.2183425328 +50,75, 1.203812008, 1.2307112287 +50,80, 1.227890346, 1.2398251648 +50,85, 1.242426337, 1.2454069243 +50,90, 1.247286586, 1.2472865857 +55, 0, 0.000000000, 1.0000000000 +55, 5, 0.112912907, 1.0024305914 +55,10, 0.225064618, 1.0096488003 +55,15, 0.335693043, 1.0214361311 +55,20, 0.444034769, 1.0374356974 +55,25, 0.549325515, 1.0571629130 +55,30, 0.650801843, 1.0800200285 +55,35, 0.747704387, 1.1053140947 +55,40, 0.839282749, 1.1322778297 +55,45, 0.924802089, 1.1600927802 +55,50, 1.003551297, 1.1879140899 +55,55, 1.074852509, 1.2148961356 +55,60, 1.138071621, 1.2402182552 +55,65, 1.192629342, 1.2631097835 +55,70, 1.238012299, 1.2828736204 +55,75, 1.273783626, 1.2989075994 +55,80, 1.299592533, 1.3107229838 +55,85, 1.315182322, 1.3179595033 +55,90, 1.320396454, 1.3203964540 +60, 0, 0.000000000, 1.0000000000 +60, 5, 0.119681778, 1.0031385295 +60,10, 0.238614577, 1.0124594672 +60,15, 0.356044091, 1.0276816504 +60,20, 0.471206153, 1.0483457003 +60,25, 0.583323727, 1.0738276019 +60,30, 0.691606043, 1.1033571989 +60,35, 0.795250355, 1.1360411010 +60,40, 0.893446594, 1.1708893642 +60,45, 0.985384972, 1.2068451910 +60,50, 1.070266403, 1.2428167937 +60,55, 1.147315349, 1.2777104815 +60,60, 1.215794546, 1.3104639783 +60,65, 1.275020900, 1.3400789457 +60,70, 1.324381718, 1.3656516965 +60,75, 1.363350417, 1.3864011169 +60,80, 1.391500813, 1.4016928947 +60,85, 1.408519209, 1.4110592570 +60,90, 1.414213562, 1.4142135624 +65, 0, 0.000000000, 1.0000000000 +65, 5, 0.128148474, 1.0040692257 +65,10, 0.255589564, 1.0161550083 +65,15, 0.381603032, 1.0358951569 +65,20, 0.505444270, 1.0626975825 +65,25, 0.626335361, 1.0957573598 +65,30, 0.743459784, 1.1340800433 +65,35, 0.855961570, 1.1765106705 +65,40, 0.962949380, 1.2217677148 +65,45, 1.063505669, 1.2684810938 +65,50, 1.156700687, 1.3152331927 +65,55, 1.241610747, 1.3606017261 +65,60, 1.317339855, 1.4032031647 +65,65, 1.383043549, 1.4417353793 +65,70, 1.437953601, 1.4750181348 +65,75, 1.481402159, 1.5020300916 +65,80, 1.512843876, 1.5219410514 +65,85, 1.531874716, 1.5341383232 +65,90, 1.538246269, 1.5382462687 +70, 0, 0.000000000, 1.0000000000 +70, 5, 0.139041489, 1.0053444028 +70,10, 0.277476571, 1.0212195717 +70,15, 0.414672740, 1.0471556657 +70,20, 0.549947578, 1.0823838086 +70,25, 0.682549331, 1.1258571388 +70,30, 0.811643704, 1.1762797795 +70,35, 0.936308263, 1.2321431946 +70,40, 1.055535305, 1.2917691861 +70,45, 1.168243466, 1.3533585717 +70,50, 1.273297730, 1.4150443413 +70,55, 1.369536895, 1.4749478592 +70,60, 1.455807011, 1.5312364694 +70,65, 1.530998883, 1.5821806891 +70,70, 1.594087380, 1.6262090720 +70,75, 1.644170149, 1.6619587940 +70,80, 1.680503336, 1.6883200831 +70,85, 1.702532036, 1.7044727784 +70,90, 1.709913565, 1.7099135651 +75, 0, 0.000000000, 1.0000000000 +75, 5, 0.153720475, 1.0072088997 +75,10, 0.307065715, 1.0286279374 +75,15, 0.459609511, 1.0636390673 +75,20, 0.610827702, 1.1112286903 +75,25, 0.760058920, 1.1700124008 +75,30, 0.906476281, 1.2382696285 +75,35, 1.049072506, 1.3139880140 +75,40, 1.186660037, 1.3949171251 +75,45, 1.317886740, 1.4786307744 +75,50, 1.441266644, 1.5625967789 +75,55, 1.555224175, 1.6442525175 +75,60, 1.658149352, 1.7210841609 +75,65, 1.748460610, 1.7907070015 +75,70, 1.824671332, 1.8509439670 +75,75, 1.885455864, 1.8998992030 +75,80, 1.929710721, 1.9360235909 +75,85, 1.956606998, 1.9581692561 +75,90, 1.965630511, 1.9656305108 +80, 0, 0.000000000, 1.0000000000 +80, 5, 0.175223596, 1.0102606485 +80,10, 0.350639262, 1.0407643440 +80,15, 0.526335260, 1.0906807598 +80,20, 0.702199693, 1.1586411101 +80,25, 0.877838622, 1.2427619421 +80,30, 1.052514778, 1.3406805139 +80,35, 1.225111680, 1.4496033094 +80,40, 1.394126403, 1.5663690138 +80,45, 1.557692334, 1.6875266770 +80,50, 1.713631283, 1.8094288493 +80,55, 1.859532258, 1.9283382823 +80,60, 1.992852358, 2.0405454606 +80,65, 2.111033523, 2.1424929245 +80,70, 2.211627685, 2.2309012139 +80,75, 2.292422061, 2.3028904563 +80,80, 2.351556149, 2.3560912550 +80,85, 2.387622438, 2.3887386793 +80,90, 2.399743837, 2.3997438370 +85, 0, 0.000000000, 1.0000000000 +85, 5, 0.213217690, 1.0166388247 +85,10, 0.428443440, 1.0661838299 +85,15, 0.647434941, 1.1475159063 +85,20, 0.871464767, 1.2587562174 +85,25, 1.101116239, 1.3972525218 +85,30, 1.336123616, 1.5595726706 +85,35, 1.575268297, 1.7415157980 +85,40, 1.816339939, 1.9381519599 +85,45, 2.056167815, 2.1438995792 +85,50, 2.290723417, 2.3526471220 +85,55, 2.515290558, 2.5579212198 +85,60, 2.724694161, 2.7530984351 +85,65, 2.913574159, 2.9316525995 +85,70, 3.076686743, 3.0874247870 +85,75, 3.209212227, 3.2148991220 +85,80, 3.307047313, 3.3094652989 +85,85, 3.367059918, 3.3676482512 +85,90, 3.387287004, 3.3872870037 diff --git a/test/data/table_17_1.csv b/test/data/table_17_1.csv new file mode 100644 index 0000000..f3017f0 --- /dev/null +++ b/test/data/table_17_1.csv @@ -0,0 +1,102 @@ +m,K,E +0.00,1.570796326794897,1.570796327 +0.01,1.574745561517356,1.566861942 +0.02,1.578739912007773,1.562912645 +0.03,1.582780342406373,1.558948244 +0.04,1.586867847454166,1.554968546 +0.05,1.591003453790792,1.550973352 +0.06,1.595188221321610,1.546962456 +0.07,1.599423244658510,1.542935653 +0.08,1.603709654639253,1.538892730 +0.09,1.608048619930513,1.534833465 +0.10,1.612441348720219,1.530757637 +0.11,1.616889090505203,1.526665017 +0.12,1.621393137980658,1.522555369 +0.13,1.625954829038433,1.518428454 +0.14,1.630575548881754,1.514284027 +0.15,1.635256732264580,1.510121831 +0.16,1.639999865864511,1.505941612 +0.17,1.644806490798881,1.501743101 +0.18,1.649678205294514,1.497526026 +0.19,1.654616667522527,1.493290109 +0.20,1.659623598610528,1.489035058 +0.21,1.664700785845692,1.484760581 +0.22,1.669850086083368,1.480466375 +0.23,1.675073429377219,1.476152126 +0.24,1.680372822848361,1.471817514 +0.25,1.685750354812596,1.467462209 +0.26,1.691208199186631,1.463085873 +0.27,1.696748620196168,1.458688155 +0.28,1.702373977410990,1.454268698 +0.29,1.708086731134606,1.449827128 +0.30,1.713889448178791,1.445363064 +0.31,1.719784808056405,1.440876115 +0.32,1.725775609629320,1.436365871 +0.33,1.731864778252098,1.431831919 +0.34,1.738055373456358,1.427273821 +0.35,1.744350597225613,1.422691133 +0.36,1.750753802915753,1.418083394 +0.37,1.757268504882456,1.413450127 +0.38,1.763898388883731,1.408790839 +0.39,1.770647323333534,1.404105019 +0.40,1.777519371491253,1.399392139 +0.41,1.784518804681873,1.394651652 +0.42,1.791650116652966,1.389882992 +0.43,1.798918039187685,1.385085568 +0.44,1.806327559107699,1.380258774 +0.45,1.813883936816983,1.375401972 +0.46,1.821592726556821,1.370514505 +0.47,1.829459798564730,1.365595691 +0.48,1.837491363355796,1.360644814 +0.49,1.845693998374724,1.355661135 +0.50,1.854074677301372,1.350643881 +0.51,1.862640802332739,1.345592245 +0.52,1.871400239811034,1.340505388 +0.53,1.880361359622178,1.335382430 +0.54,1.889533078853096,1.330222453 +0.55,1.898924910271554,1.325024498 +0.56,1.908547016281211,1.319787557 +0.57,1.918410269109912,1.314510576 +0.58,1.928526318114418,1.309192448 +0.59,1.938907665234220,1.303832008 +0.60,1.949567749806026,1.298428034 +0.61,1.960521044165830,1.292979239 +0.62,1.971783161725656,1.287484262 +0.63,1.983370979527821,1.281941668 +0.64,1.995302777664729,1.276349943 +0.65,2.007598398424376,1.270707480 +0.66,2.020279428603592,1.265012576 +0.67,2.033369409152233,1.259263421 +0.68,2.046894077210577,1.253458093 +0.69,2.060881646730131,1.247594538 +0.70,2.075363135292469,1.241670567 +0.71,2.090372746552360,1.235683836 +0.72,2.105948320052758,1.229631828 +0.73,2.122131863157396,1.223511839 +0.74,2.138970183752114,1.217320955 +0.75,2.156515647499643,1.211056028 +0.76,2.174827090246414,1.204713641 +0.77,2.193970925319189,1.198290087 +0.78,2.214022497846332,1.191781311 +0.79,2.235067755260349,1.185182883 +0.80,2.257205326820854,1.178489924 +0.81,2.280549138422770,1.171697053 +0.82,2.305231736877189,1.164798293 +0.83,2.331408567750251,1.157786979 +0.84,2.359263554745007,1.150655629 +0.85,2.389016486325580,1.143395792 +0.86,2.420932960344303,1.135997843 +0.87,2.455338028321380,1.128450735 +0.88,2.492635323239716,1.120741661 +0.89,2.533334546002200,1.112855607 +0.90,2.578092113348173,1.104774733 +0.91,2.627773332084344,1.096477517 +0.92,2.683551406315229,1.087937503 +0.93,2.747073004024667,1.079121407 +0.94,2.820752496755872,1.069986130 +0.95,2.908337248444552,1.060473728 +0.96,3.016112492477648,1.050502227 +0.97,3.155874947891841,1.039946861 +0.98,3.354141445699160,1.028594520 +0.99,3.695637362989875,1.015993546 +1.00,Inf,1.000000000 diff --git a/test/data/table_17_2.csv b/test/data/table_17_2.csv new file mode 100644 index 0000000..ac719ae --- /dev/null +++ b/test/data/table_17_2.csv @@ -0,0 +1,92 @@ +α,K,E + 0,1.570796326794897,1.570796326794897 + 1,1.570915958127243,1.570676709127960 + 2,1.571274952372225,1.570317919897448 + 3,1.571873610514009,1.569720150423979 + 4,1.572712434995227,1.568883719607763 + 5,1.573792130924768,1.567809073977622 + 6,1.575113607777251,1.566496787760132 + 7,1.576677981592838,1.564947562969419 + 8,1.578486577688648,1.563162229518261 + 9,1.580540933895721,1.561141745351334 +10,1.582842804338351,1.558887196601596 +11,1.585394163775538,1.556399797770947 +12,1.588197212527520,1.553680891936509 +13,1.591254382013687,1.550731950984013 +14,1.594568340931825,1.547554575869993 +15,1.598142002112540,1.544150496914673 +16,1.601978530086952,1.540521574127631 +17,1.606081349410364,1.536669797568556 +18,1.610454153789663,1.532597287745636 +19,1.615100916067722,1.528306296054359 +20,1.620025899124204,1.523799205259774 +21,1.625233667758843,1.519078530025531 +22,1.630729101630788,1.514146917493342 +23,1.636517409335819,1.509007147916775 +24,1.642604143712491,1.503662135353715 +25,1.648995218478530,1.498114928422116 +26,1.655696926310344,1.492368711124151 +27,1.662715958491370,1.486426803744253 +28,1.670059426269580,1.480292663827039 +29,1.677734884080745,1.473969887241625 +30,1.685750354812596,1.467462209339427 +31,1.694114357305914,1.460773506213127 +32,1.702835936312341,1.453907796065210 +33,1.711924695155678,1.446869240695183 +34,1.721390831374249,1.439662147115459 +35,1.731245175657058,1.432290969306756 +36,1.741499234426774,1.424760310124890 +37,1.752165236468845,1.417074923371952 +38,1.763256184059342,1.409239716046096 +39,1.774785909105608,1.401259750785523 +40,1.786769134885021,1.393140248523812 +41,1.799221544049811,1.384886591375413 +42,1.812159853662126,1.376504325772082 +43,1.825601898135889,1.367999165873159 +44,1.839566721093652,1.359376997275008 +45,1.854074677301372,1.350643881047676 +46,1.869147546026462,1.341806058129911 +47,1.884808657380404,1.332869954117179 +48,1.901083033463664,1.323842184481263 +49,1.917997546436423,1.314729560264623 +50,1.935581096004722,1.305539094297794 +51,1.953864809251663,1.296278007994134 +52,1.972882266274650,1.286953738783001 +53,1.992669755734209,1.277573948250391 +54,2.013266565205468,1.268146531065206 +55,2.034715312185791,1.258679624779997 +56,2.057062322797365,1.249181620607472 +57,2.080358066691578,1.239661175288672 +58,2.104657658491159,1.230127224185949 +59,2.130021438399325,1.220588995754247 +60,2.156515647499643,1.211056027568459 +61,2.184213216949248,1.201538184113662 +62,2.213194694979374,1.192045676579886 +63,2.243549341698626,1.182589084945384 +64,2.275376429611676,1.173179382683722 +65,2.308786798167196,1.163827964493139 +66,2.343904724446913,1.154546677524465 +67,2.380870190604429,1.145347856680849 +68,2.419841653739137,1.136244364684239 +69,2.460999458304126,1.127249637757702 +70,2.504550079001634,1.118377737969864 +71,2.550731449627254,1.109643413542761 +72,2.599819730061099,1.101062168757941 +73,2.652138004630204,1.092650345537715 +74,2.708067614590486,1.084425219372543 +75,2.768063145368768,1.076405113076403 +76,2.832672582918100,1.068609532978401 +77,2.902564940670027,1.061059333753857 +78,2.978568951181384,1.053776920407046 +79,3.061728612038789,1.046786499344049 +80,3.153385251887839,1.040114395706010 +81,3.255302942143555,1.033789462390754 +82,3.369868026668445,1.027843619740833 +83,3.500422499171838,1.022312588167584 +84,3.651855969478752,1.017236918341019 +85,3.831741999784146,1.012663506234396 +86,4.052758169549437,1.008647956907096 +87,4.338653975999725,1.005258587209152 +88,4.742717265278886,1.002584085527552 +89,5.434909829625564,1.000751577701834 +90,Inf,1.000000000000000 diff --git a/test/runtests.jl b/test/runtests.jl index 1cb38e3..a9887ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -77,3 +77,5 @@ end @test Bz ≈ μ0 * I * R^2 / (2 * (z^2 + R^2)^1.5) end end + +include("runtests_elliptic.jl") diff --git a/test/runtests_elliptic.jl b/test/runtests_elliptic.jl new file mode 100644 index 0000000..29e7b94 --- /dev/null +++ b/test/runtests_elliptic.jl @@ -0,0 +1,720 @@ +using Test +using SpecialFunctions: gamma +using LoopFieldCalc.Elliptic2 +using DelimitedFiles: readdlm + +@testset "NaNs" begin + @test E(NaN) === NaN + @test K(NaN) === NaN + + @test E(0,NaN) === NaN + @test E(NaN,0) === NaN + @test E(NaN,NaN) === NaN + + @test F(0,NaN) === NaN + @test F(NaN,0) === NaN + @test F(NaN,NaN) === NaN + + @test Pi(0,0,NaN) === NaN + @test Pi(0,NaN,0) === NaN + @test Pi(NaN,0,0) === NaN + @test Pi(NaN,NaN,NaN) === NaN + + @test ellipke(NaN) === (NaN,NaN) +end + +# test negative arguments +@testset "Negative arguments" begin + for m in 1:10 + e = E(m/(m + 1))*sqrt(m + 1) + k = K(m/(m + 1))/sqrt(m + 1) + @test E(-m) ≈ e + @test K(-m) ≈ k + k2,e2 = ellipke(-m) + @test e2 ≈ e + @test k2 ≈ k + end +end + +# +# K(m), E(m), ellipke(m), F(pi/2,m), E(pi/2,m) +# + +@testset "Table 17.1: K(m) and E(m)" begin + # values from Abramowitz & Stegun, Table 17.1 (p608-609) + # m K(m) K(1-m) E(m) E(1-m) + table17p1 = [ + 0.00 1.57079_63267_94897 Inf 1.57079_6327 1.00000_0000; + 0.01 1.57474_55615_17356 3.69563_73629_89875 1.56686_1942 1.01599_3546; + 0.02 1.57873_99120_07773 3.35414_14456_99160 1.56291_2645 1.02859_4520; + 0.03 1.58278_03424_06373 3.15587_49478_91841 1.55894_8244 1.03994_6861; + 0.04 1.58686_78474_54166 3.01611_24924_77648 1.55496_8546 1.05050_2227; + 0.05 1.59100_34537_90792 2.90833_72484_44552 1.55097_3352 1.06047_3728; + 0.06 1.59518_82213_21610 2.82075_24967_55872 1.54696_2456 1.06998_6130; + 0.07 1.59942_32446_58510 2.74707_30040_24667 1.54293_5653 1.07912_1407; + 0.08 1.60370_96546_39253 2.68355_14063_15229 1.53889_2730 1.08793_7503; + 0.09 1.60804_86199_30513 2.62777_33320_84344 1.53483_3465 1.09647_7517; + 0.10 1.61244_13487_20219 2.57809_21133_48173 1.53075_7637 1.10477_4733; + 0.11 1.61688_90905_05203 2.53333_45460_02200 1.52666_5017 1.11285_5607; + 0.12 1.62139_31379_80658 2.49263_53232_39716 1.52255_5369 1.12074_1661; + 0.13 1.62595_48290_38433 2.45533_80283_21380 1.51842_8454 1.12845_0735; + 0.14 1.63057_55488_81754 2.42093_29603_44303 1.51428_4027 1.13599_7843; + 0.15 1.63525_67322_64580 2.38901_64863_25580 1.51012_1831 1.14339_5792; + 0.16 1.63999_98658_64511 2.35926_35547_45007 1.50594_1612 1.15065_5629; + 0.17 1.64480_64907_98881 2.33140_85677_50251 1.50174_3101 1.15778_6979; + 0.18 1.64967_82052_94514 2.30523_17368_77189 1.49752_6026 1.16479_8293; + 0.19 1.65461_66675_22527 2.28054_91384_22770 1.49329_0109 1.17169_7053; + 0.20 1.65962_35986_10528 2.25720_53268_20854 1.48903_5058 1.17848_9924; + 0.21 1.66470_07858_45692 2.23506_77552_60349 1.48476_0581 1.18518_2883; + 0.22 1.66985_00860_83368 2.21402_24978_46332 1.48046_6375 1.19178_1311; + 0.23 1.67507_34293_77219 2.19397_09253_19189 1.47615_2126 1.19829_0087; + 0.24 1.68037_28228_48361 2.17482_70902_46414 1.47181_7514 1.20471_3641; + 0.25 1.68575_03548_12596 2.15651_56474_99643 1.46746_2209 1.21105_6028; + 0.26 1.69120_81991_86631 2.13897_01837_52114 1.46308_5873 1.21732_0955; + 0.27 1.69674_86201_96168 2.12213_18631_57396 1.45868_8155 1.22351_1839; + 0.28 1.70237_39774_10990 2.10594_83200_52758 1.45426_8698 1.22963_1828; + 0.29 1.70808_67311_34606 2.09037_27465_52360 1.44982_7128 1.23568_3836; + 0.30 1.71388_94481_78791 2.07536_31352_92469 1.44536_3064 1.24167_0567; + 0.31 1.71978_48080_56405 2.06088_16467_30131 1.44087_6115 1.24759_4538; + 0.32 1.72577_56096_29320 2.04689_40772_10577 1.43636_5871 1.25345_8093; + 0.33 1.73186_47782_52098 2.03336_94091_52233 1.43183_1919 1.25926_3421; + 0.34 1.73805_53734_56358 2.02027_94286_03592 1.42727_3821 1.26501_2576; + 0.35 1.74435_05972_25613 2.00759_83984_24376 1.42269_1133 1.27070_7480; + 0.36 1.75075_38029_15753 1.99530_27776_64729 1.41808_3394 1.27634_9943; + 0.37 1.75726_85048_82456 1.98337_09795_27821 1.41345_0127 1.28194_1668; + 0.38 1.76389_83888_83731 1.97178_31617_25656 1.40879_0839 1.28748_4262; + 0.39 1.77064_73233_33534 1.96052_10441_65830 1.40410_5019 1.29297_9239; + 0.40 1.77751_93714_91253 1.94956_77498_06026 1.39939_2139 1.29842_8034; + 0.41 1.78451_88046_81873 1.93890_76652_34220 1.39465_1652 1.30383_2008; + 0.42 1.79165_01166_52966 1.92852_63181_14418 1.38988_2992 1.30919_2448; + 0.43 1.79891_80391_87685 1.91841_02691_09912 1.38508_5568 1.31451_0576; + 0.44 1.80632_75591_07699 1.90854_70162_81211 1.38025_8774 1.31978_7557; + 0.45 1.81388_39368_16983 1.89892_49102_71554 1.37540_1972 1.32502_4498; + 0.46 1.82159_27265_56821 1.88953_30788_53096 1.37051_4505 1.33022_2453; + 0.47 1.82945_97985_64730 1.88036_13596_22178 1.36559_5691 1.33538_2430; + 0.48 1.83749_13633_55796 1.87140_02398_11034 1.36064_4814 1.34050_5388; + 0.49 1.84569_39983_74724 1.86264_08023_32739 1.35566_1135 1.34559_2245; + 0.50 1.85407_46773_01372 1.85407_46773_01372 1.35064_3881 1.35064_3881; + ] + + etol = 2e-9 + ktol = 2e-15 + for i = 1:size(table17p1,1) + m = table17p1[i,1] + k = table17p1[i,2] + e = table17p1[i,4] + @test K(m) ≈ k atol=ktol + @test E(m) ≈ e atol=etol + @test F(pi/2,m) ≈ k atol=ktol + @test E(pi/2,m) ≈ e atol=etol + a,b = ellipke(m) + @test a ≈ k atol=ktol + @test b ≈ e atol=etol + end + + for i = 1:size(table17p1,1) + m = 1. - table17p1[i,1] + k = table17p1[i,3] + e = table17p1[i,5] + a,b = ellipke(m) + if k == Inf + @test K(m) == k + @test F(pi/2,m) == k + @test a == k + else + @test K(m) ≈ k atol=ktol + @test F(pi/2,m) ≈ k atol=ktol + @test a ≈ k atol=ktol + end + @test E(m) ≈ e atol=etol + @test E(pi/2,m) ≈ e atol=etol + @test b ≈ e atol=etol + end + + @test K(0) ≈ pi/2 + @test K(0.5) ≈ (0.25/sqrt(pi))*gamma(0.25)^2 + @test K(1) == Inf + + @test E(0) ≈ pi/2 + @test E(1) ≈ 1 +end + +@testset "Table 17.2: K(α) and E(α)" begin + # Abramowitz & Stegun, Table 17.2, p610-11 + # α K(α) K(90-α) E(α) E(90-α) + table17p2 = [ + 0. 1.57079_63267_94897 Inf 1.57079_63267_94897 1.00000_00000_00000; + 1 1.57091_59581_27243 5.43490_98296_25564 1.57067_67091_27960 1.00075_15777_01834; + 2 1.57127_49523_72225 4.74271_72652_78886 1.57031_79198_97448 1.00258_40855_27552; + 3 1.57187_36105_14009 4.33865_39759_99725 1.56972_01504_23979 1.00525_85872_09152; + 4 1.57271_24349_95227 4.05275_81695_49437 1.56888_37196_07763 1.00864_79569_07096; + 5 1.57379_21309_24768 3.83174_19997_84146 1.56780_90739_77622 1.01266_35062_34396; + 6 1.57511_36077_77251 3.65185_59694_78752 1.56649_67877_60132 1.01723_69183_41019; + 7 1.57667_79815_92838 3.50042_24991_71838 1.56494_75629_69419 1.02231_25881_67584; + 8 1.57848_65776_88648 3.36986_80266_68445 1.56316_22295_18261 1.02784_36197_40833; + 9 1.58054_09338_95721 3.25530_29421_43555 1.56114_17453_51334 1.03378_94623_90754; + 10 1.58284_28043_38351 3.15338_52518_87839 1.55888_71966_01596 1.04011_43957_06010; + 11 1.58539_41637_75538 3.06172_86120_38789 1.55639_97977_70947 1.04678_64993_44049; + 12 1.58819_72125_27520 2.97856_89511_81384 1.55368_08919_36509 1.05377_69204_07046; + 13 1.59125_43820_13687 2.90256_49406_70027 1.55073_19509_84013 1.06105_93337_53857; + 14 1.59456_83409_31825 2.83267_25829_18100 1.54755_45758_69993 1.06860_95329_78401; + 15 1.59814_20021_12540 2.76806_31453_68768 1.54415_04969_14673 1.07640_51130_76403; + 16 1.60197_85300_86952 2.70806_76145_90486 1.54052_15741_27631 1.08442_52193_72543; + 17 1.60608_13494_10364 2.65213_80046_30204 1.53666_97975_68556 1.09265_03455_37715; + 18 1.61045_41537_89663 2.59981_97300_61099 1.53259_72877_45636 1.10106_21687_57941; + 19 1.61510_09160_67722 2.55073_14496_27254 1.52830_62960_54359 1.10964_34135_42761; + 20 1.62002_58991_24204 2.50455_00790_01634 1.52379_92052_59774 1.11837_77379_69864; + 21 1.62523_36677_58843 2.46099_94583_04126 1.51907_85300_25531 1.12724_96377_57702; + 22 1.63072_91016_30788 2.41984_16537_39137 1.51414_69174_93342 1.13624_43646_84239; + 23 1.63651_74093_35819 2.38087_01906_04429 1.50900_71479_16775 1.14534_78566_80849; + 24 1.64260_41437_12491 2.34390_47244_46913 1.50366_21353_53715 1.15454_66775_24465; + 25 1.64899_52184_78530 2.30878_67981_67196 1.49811_49284_22116 1.16382_79644_93139; + 26 1.65569_69263_10344 2.27537_64296_11676 1.49236_87111_24151 1.17317_93826_83722; + 27 1.66271_59584_91370 2.24354_93416_98626 1.48642_68037_44253 1.18258_90849_45384; + 28 1.67005_94262_69580 2.21319_46949_79374 1.48029_26638_27039 1.19204_56765_79886; + 29 1.67773_48840_80745 2.18421_32169_49248 1.47396_98872_41625 1.20153_81841_13662; + 30 1.68575_03548_12596 2.15651_56474_99643 1.46746_22093_39427 1.21105_60275_68459; + 31 1.69411_43573_05914 2.13002_14383_99325 1.46077_35062_13127 1.22058_89957_54247; + 32 1.70283_59363_12341 2.10465_76584_91159 1.45390_77960_65210 1.23012_72241_85949; + 33 1.71192_46951_55678 2.08035_80666_91578 1.44686_92406_95183 1.23966_11752_88672; + 34 1.72139_08313_74249 2.05706_23227_97365 1.43966_21471_15459 1.24918_16206_07472; + 35 1.73124_51756_57058 2.03471_53121_85791 1.43229_09693_06756 1.25867_96247_79997; + 36 1.74149_92344_26774 2.01326_65652_05468 1.42476_03101_24890 1.26814_65310_65206; + 37 1.75216_52364_68845 1.99266_97557_34209 1.41707_49233_71952 1.27757_39482_50391; + 38 1.76325_61840_59342 1.97288_22662_74650 1.40923_97160_46096 1.28695_37387_83001; + 39 1.77478_59091_05608 1.95386_48092_51663 1.40125_97507_85523 1.29627_80079_94134; + 40 1.78676_91348_85021 1.93558_10960_04722 1.39314_02485_23812 1.30553_90942_97794; + 41 1.79922_15440_49811 1.91799_75464_36423 1.38488_65913_75413 1.31472_95602_64623; + 42 1.81215_98536_62126 1.90108_30334_63664 1.37650_43257_72082 1.32384_21844_81263; + 43 1.82560_18981_35889 1.88480_86573_80404 1.36799_91658_73159 1.33286_99541_17179; + 44 1.83956_67210_93652 1.86914_75460_26462 1.35937_69972_75008 1.34180_60581_29911; + 45 1.85407_46773_01372 1.85407_46773_01372 1.35064_38810_47676 1.35064_38810_47676; + ] + + tol = 2e-15 + for i = 1:size(table17p2,1) + alpha = table17p2[i,1] + k = table17p2[i,2] + e = table17p2[i,4] + m = sind(alpha)^2 + @test K(m) ≈ k atol=tol + @test E(m) ≈ e atol=tol + @test F(pi/2,m) ≈ k atol=tol + @test E(pi/2,m) ≈ e atol=tol + a,b = ellipke(m) + @test a ≈ k atol=tol + @test b ≈ e atol=tol + end + + # XXX: why isn't tol smaller? + tol = 1e-12 + for i = 1:size(table17p2,1) + alpha = 90. - table17p2[i,1] + k = table17p2[i,3] + e = table17p2[i,5] + m = sind(alpha)^2 + a,b = ellipke(m) + if k == Inf + @test K(m) == k + @test F(pi/2,m) == k + @test a == k + else + @test K(m) ≈ k atol=tol + @test F(pi/2,m) ≈ k atol=tol + @test a ≈ k atol=tol + end + + @test E(m) ≈ e atol=tol + @test E(pi/2,m) ≈ e atol=tol + @test b ≈ e atol=tol + end +end + +@testset "Table 17.5: F(φ, m)" begin + # Abramowitz & Stegun, Table 17.5, p613 + # F(phi\alpha) + # α\φ 0 5 10 15 20 25 30 + table17p5_p613 = [ + 0. 0 0.08726_646 0.17453_293 0.26179_939 0.34906_585 0.43633_231 0.52359_878; + 2 0 0.08726_660 0.17453_400 0.26180_298 0.34907_428 0.43634_855 0.52362_636; + 4 0 0.08726_700 0.17453_721 0.26181_374 0.34909_952 0.43639_719 0.52370_903; + 6 0 0.08726_767 0.17454_255 0.26183_163 0.34914_148 0.43647_806 0.52384_653; + 8 0 0.08726_860 0.17454_999 0.26185_656 0.34919_998 0.43659_086 0.52403_839; + 10 0 0.08726_980 0.17455_949 0.26188_842 0.34927_479 0.43673_518 0.52428_402; + 12 0 0.08727_124 0.17457_102 0.26192_707 0.34936_558 0.43691_046 0.52458_259; + 14 0 0.08727_294 0.17458_451 0.26197_234 0.34947_200 0.43711_606 0.52493_314; + 16 0 0.08727_487 0.17459_991 0.26202_402 0.34959_358 0.43735_119 0.52533_449; + 18 0 0.08727_703 0.17461_714 0.26208_189 0.34972_983 0.43761_496 0.52578_529; + 20 0 0.08727_940 0.17463_611 0.26214_568 0.34988_016 0.43790_635 0.52628_399; + 22 0 0.08728_199 0.17465_675 0.26221_511 0.35004_395 0.43822_422 0.52682_887; + 24 0 0.08728_477 0.17467_895 0.26228_985 0.35022_048 0.43856_733 0.52741_799; + 26 0 0.08728_773 0.17470_261 0.26236_958 0.35040_901 0.43893_430 0.52804_924; + 28 0 0.08729_086 0.17472_762 0.26245_392 0.35060_870 0.43932_365 0.52872_029; + 30 0 0.08729_413 0.17475_386 0.26254_249 0.35081_868 0.43973_377 0.52942_863; + 32 0 0.08729_755 0.17478_119 0.26263_487 0.35103_803 0.44016_296 0.53017_153; + 34 0 0.08730_108 0.17480_950 0.26273_064 0.35126_576 0.44060_939 0.53094_608; + 36 0 0.08730_472 0.17483_864 0.26282_934 0.35150_083 0.44107_115 0.53174_916; + 38 0 0.08730_844 0.17486_848 0.26293_052 0.35174_218 0.44154_622 0.53257_745; + 40 0 0.08731_222 0.17489_887 0.26303_369 0.35198_869 0.44203_247 0.53342_745; + 42 0 0.08731_606 0.17492_967 0.26313_836 0.35223_920 0.44252_769 0.53429_546; + 44 0 0.08731_992 0.17496_073 0.26324_404 0.35249_254 0.44302_960 0.53517_761; + 46 0 0.08732_379 0.17499_189 0.26335_019 0.35274_748 0.44353_584 0.53606_986; + 48 0 0.08732_765 0.17502_300 0.26345_633 0.35300_280 0.44404_397 0.53696_798; + 50 0 0.08733_149 0.17505_392 0.26356_191 0.35325_724 0.44455_151 0.53786_765; + 52 0 0.08733_528 0.17508_448 0.26366_643 0.35350_955 0.44505_593 0.53876_438; + 54 0 0.08733_901 0.17511_455 0.26376_936 0.35375_845 0.44555_469 0.53965_358; + 56 0 0.08734_265 0.17514_397 0.26387_020 0.35400_269 0.44604_519 0.54053_059; + 58 0 0.08734_620 0.17517_260 0.26396_842 0.35424_101 0.44652_487 0.54139_069; + 60 0 0.08734_962 0.17520_029 0.26406_355 0.35447_217 0.44699_117 0.54222_911; + 62 0 0.08735_291 0.17522_690 0.26415_509 0.35469_497 0.44744_153 0.54304_111; + 64 0 0.08735_605 0.17525_232 0.26424_258 0.35490_823 0.44787_348 0.54382_197; + 66 0 0.08735_902 0.17527_640 0.26432_556 0.35511_081 0.44828_459 0.54456_704; + 68 0 0.08736_182 0.17529_903 0.26440_362 0.35530_160 0.44867_252 0.54527_182; + 70 0 0.08736_442 0.17532_010 0.26447_634 0.35547_959 0.44903_502 0.54593_192; + 72 0 0.08736_681 0.17533_949 0.26454_334 0.35564_377 0.44936_997 0.54654_316; + 74 0 0.08736_898 0.17535_712 0.26460_428 0.35579_326 0.44967_538 0.54710_162; + 76 0 0.08737_092 0.17537_289 0.26465_883 0.35592_721 0.44994_944 0.54760_364; + 78 0 0.08737_262 0.17538_672 0.26470_671 0.35604_488 0.45019_046 0.54804_587; + 80 0 0.08737_408 0.17539_854 0.26474_766 0.35614_560 0.45039_699 0.54842_535; + 82 0 0.08737_528 0.17540_830 0.26478_147 0.35622_881 0.45056_775 0.54873_947; + 84 0 0.08737_622 0.17541_594 0.26480_795 0.35629_402 0.45070_168 0.54898_608; + 86 0 0.08737_689 0.17542_143 0.26482_697 0.35634_086 0.45079_795 0.54916_348; + 88 0 0.08737_730 0.17542_473 0.26483_842 0.35636_908 0.45085_596 0.54927_042; + 90 0 0.08737_744 0.17542_583 0.26484_225 0.35637_851 0.45087_533 0.54930_614; + ] + + # Abramowitz & Stegun, Table 17.5, p614 + # F(phi\alpha) + # α\φ 35 40 45 50 55 60 + table17p5_p614 = [ + 0. 0.61086_524 0.69813_170 0.78539_816 0.87266_463 0.95993_109 1.04719_755; + 2 0.61090_819 0.69819_436 0.78548_509 0.87278_045 0.96008_037 1.04738_465; + 4 0.61103_691 0.69838_220 0.78574_574 0.87312_784 0.96052_821 1.04794_603; + 6 0.61125_108 0.69869_484 0.78617_974 0.87370_649 0.96127_450 1.04888_194; + 8 0.61155_010 0.69913_161 0.78678_644 0.87451_593 0.96231_911 1.05019_278; + 10 0.61193_318 0.69969_159 0.78756_494 0.87555_545 0.96366_180 1.05187_911; + 12 0.61239_927 0.70037_358 0.78851_403 0.87682_412 0.96530_224 1.05394_160; + 14 0.61294_707 0.70117_608 0.78963_221 0.87832_076 0.96723_998 1.05638_099; + 16 0.61357_504 0.70209_730 0.79091_768 0.88004_389 0.96947_438 1.05919_813; + 18 0.61428_140 0.70313_511 0.79236_827 0.88199_174 0.97200_462 1.06239_384; + 20 0.61506_406 0.70428_706 0.79398_143 0.88416_214 0.97482_960 1.06596_891; + 22 0.61592_071 0.70555_037 0.79575_422 0.88655_254 0.97794_790 1.06992_405; + 24 0.61684_871 0.70692_183 0.79768_324 0.88915_992 0.98135_773 1.07425_976; + 26 0.61784_515 0.70839_788 0.79976_461 0.89198_071 0.98505_681 1.07897_628; + 28 0.61890_682 0.70997_451 0.80199_389 0.89501_076 0.98904_227 1.08407_347; + 30 0.62003_018 0.71164_728 0.80436_610 0.89824_524 0.99331_059 1.08955_067; + 32 0.62121_138 0.71341_124 0.80687_558 0.90167_852 0.99785_743 1.09540_656; + 34 0.62244_622 0.71526_098 0.80951_599 0.90530_415 1.00267_749 1.10163_899; + 36 0.62373_019 0.71719_052 0.81228_024 0.90911_465 1.00776_438 1.10824_474; + 38 0.62505_840 0.71919_335 0.81516_039 0.91310_148 1.01311_039 1.11521_933; + 40 0.62642_563 0.72126_235 0.81814_765 0.91725_487 1.01870_633 1.12255_667; + 42 0.62782_630 0.72338_982 0.82123_227 0.92156_370 1.02454_127 1.13024_880; + 44 0.62925_446 0.72556_741 0.82440_346 0.92601_535 1.03060_230 1.13828_546; + 46 0.63070_385 0.72778_615 0.82764_941 0.93059_558 1.03687_427 1.14665_369; + 48 0.63216_783 0.73003_640 0.83095_712 0.93528_835 1.04333_948 1.15533_731; + 50 0.63363_947 0.73230_789 0.83431_247 0.94007_568 1.04997_735 1.16431_637; + 52 0.63511_150 0.73458_970 0.83770_010 0.94493_756 1.05676_412 1.17356_652; + 54 0.63657_639 0.73687_028 0.84110_344 0.94985_177 1.06367_248 1.18305_833; + 56 0.63802_636 0.73913_751 0.84450_468 0.95479_381 1.07067_128 1.19275_650; + 58 0.63945_343 0.74137_870 0.84788_483 0.95973_682 1.07772_516 1.20261_907; + 60 0.64084_944 0.74358_071 0.85122_375 0.96465_156 1.08479_434 1.21259_661; + 62 0.64220_613 0.74572_998 0.85450_024 0.96950_647 1.09183_436 1.22263_139; + 64 0.64351_521 0.74781_266 0.85769_220 0.97426_773 1.09879_601 1.23265_660; + 66 0.64476_839 0.74981_471 0.86077_677 0.97889_946 1.10562_535 1.24259_576; + 68 0.64595_751 0.75172_208 0.86373_057 0.98336_406 1.11226_392 1.25236_238; + 70 0.64707_458 0.75352_078 0.86652_996 0.98762_253 1.11864_920 1.26185_988; + 72 0.64811_189 0.75519_716 0.86915_135 0.99163_507 1.12471_530 1.27098_218; + 74 0.64906_209 0.75673_800 0.87157_159 0.99536_166 1.13039_401 1.27961_482; + 76 0.64991_829 0.75813_076 0.87376_830 0.99876_287 1.13561_610 1.28763_696; + 78 0.65067_415 0.75936_376 0.87572_037 1.00180_067 1.14031_304 1.29492_436; + 80 0.65132_394 0.76042_640 0.87740_833 1.00443_942 1.14441_892 1.30135_321; + 82 0.65186_270 0.76130_931 0.87881_481 1.00664_678 1.14787_262 1.30680_495; + 84 0.65228_622 0.76200_457 0.87992_495 1.00839_470 1.15062_010 1.31117_166; + 86 0.65259_116 0.76250_582 0.88072_675 1.00966_028 1.15261_652 1.31436_170; + 88 0.65277_510 0.76280_846 0.88121_143 1.01042_658 1.15382_828 1.31630_510; + 90 0.65283_658 0.76290_965 0.88137_359 1.01068_319 1.15423_455 1.31695_790; + ] + + # Abramowitz & Stegun, Table 17.5, p615 + # F(phi\alpha) + # α\φ 65 70 75 80 85 90 + table17p5_p615 = [ + 0. 1.13446_401 1.22173_048 1.30899_694 1.39626_340 1.48352_986 1.57079_633; + 2 1.13469_294 1.22200_477 1.30931_959 1.39663_672 1.48395_543 1.57127_495; + 4 1.13537_994 1.22282_810 1.31028_822 1.39775_763 1.48523_342 1.57271_244; + 6 1.13652_576 1.22420_180 1.31190_491 1.39962_909 1.48736_769 1.57511_361; + 8 1.13813_158 1.22612_810 1.31417_314 1.40225_598 1.49036_470 1.57848_658; + 10 1.14019_906 1.22861_010 1.31709_778 1.40564_522 1.49423_361 1.58284_280; + 12 1.14273_032 1.23165_180 1.32068_514 1.40980_577 1.49898_627 1.58819_721; + 14 1.14572_789 1.23525_808 1.32494_296 1.41474_871 1.50463_742 1.59456_834; + 16 1.14919_471 1.23943_470 1.32988_047 1.42048_728 1.51120_474 1.60197_853; + 18 1.15313_409 1.24418_827 1.33550_840 1.42703_700 1.51870_904 1.61045_415; + 20 1.15754_967 1.24952_627 1.34183_901 1.43441_578 1.52717_445 1.62002_590; + 22 1.16244_535 1.25545_700 1.34888_616 1.44264_399 1.53662_865 1.63072_910; + 24 1.16782_525 1.26198_957 1.35666_531 1.45174_466 1.54710_309 1.64260_414; + 26 1.17369_362 1.26913_385 1.36519_359 1.46174_360 1.55863_334 1.65569_693; + 28 1.18005_472 1.27690_045 1.37448_981 1.47266_958 1.57125_942 1.67005_943; + 30 1.18691_274 1.28530_059 1.38457_455 1.48455_455 1.58502_624 1.68575_035; + 32 1.19427_162 1.29434_605 1.39547_013 1.49743_384 1.59998_406 1.70283_594; + 34 1.20213_489 1.30404_906 1.40720_064 1.51134_644 1.61618_906 1.72139_083; + 36 1.21050_542 1.31442_210 1.41979_198 1.52633_523 1.63370_398 1.74149_923; + 38 1.21938_520 1.32547_772 1.43327_179 1.54244_734 1.65259_894 1.76325_618; + 40 1.22877_499 1.33722_824 1.44766_938 1.55973_441 1.67295_226 1.78676_913; + 42 1.23867_392 1.34968_545 1.46301_565 1.57825_301 1.69485_156 1.81215_985; + 44 1.24907_904 1.36286_013 1.47934_287 1.59806_493 1.71839_498 1.83956_672; + 46 1.25998_475 1.37676_148 1.49668_437 1.61923_762 1.74369_264 1.86914_755; + 48 1.27138_210 1.39139_640 1.51507_416 1.64184_453 1.77086_836 1.90108_303; + 50 1.28325_798 1.40676_855 1.53454_619 1.66596_542 1.80006_176 1.93558_110; + 52 1.29559_414 1.42287_717 1.55513_354 1.69168_665 1.83143_068 1.97288_227; + 54 1.30836_604 1.43971_560 1.57686_709 1.71910_125 1.86515_414 2.01326_657; + 56 1.32154_149 1.45726_935 1.59977_378 1.74830_880 1.90143_591 2.05706_232; + 58 1.33507_910 1.47551_372 1.62387_409 1.77941_482 1.94050_873 2.10465_766; + 60 1.34892_643 1.49441_087 1.64917_867 1.81252_953 1.98263_957 2.15651_565; + 62 1.36301_803 1.51390_609 1.67568_359 1.84776_547 2.02813_570 2.21319_470; + 64 1.37727_323 1.53392_332 1.70336_398 1.88523_335 2.07735_219 2.27537_643; + 66 1.39159_384 1.55435_972 1.73216_516 1.92503_509 2.13070_052 2.34390_472; + 68 1.40586_195 1.57507_940 1.76199_085 1.96725_237 2.18865_839 2.41984_165; + 70 1.41993_796 1.59590_624 1.79268_736 2.01192_798 2.25177_995 2.50455_008; + 72 1.43365_925 1.61661_644 1.82402_292 2.05903_582 2.32070_416 2.59981_973; + 74 1.44684_001 1.63693_134 1.85566_175 2.10843_282 2.39615_610 2.70806_762; + 76 1.45927_266 1.65651_218 1.88713_308 2.15978_295 2.47892_739 2.83267_258; + 78 1.47073_163 1.67495_873 1.91779_814 2.21243_977 2.56980_281 2.97856_895; + 80 1.48098_006 1.69181_489 1.94682_231 2.26527_326 2.66935_045 3.15338_525; + 82 1.48977_975 1.70658_456 1.97316_666 2.31643_897 2.77736_748 3.36986_803; + 84 1.49690_410 1.71876_033 1.99562_118 2.36313_736 2.89146_664 3.65185_597; + 86 1.50215_336 1.72786_543 2.01290_452 2.40153_358 3.00370_926 4.05275_817; + 88 1.50537_033 1.73350_464 2.02384_126 2.42718_003 3.09448_898 4.74271_727; + 90 1.50645_424 1.73541_516 2.02758_942 2.43624_605 3.13130_133 Inf; + ] + + table17p5 = [table17p5_p613 table17p5_p614[:,2:end] table17p5_p615[:,2:end]] + + for i = 1:size(table17p5,1) + alpha = table17p5[i,1] + m = sind(alpha)^2 + for j = 2:size(table17p5,2) + phid = 5*(j-2) + phi = deg2rad(phid) + f = table17p5[i,j] + if f == Inf + @test F(phi,m) == f + @test F(-phi,m) == -f + if phid == 90 + @test K(m) == f + @test ellipke(m)[1] == f + end + else + @test F(phi,m) ≈ f atol=1e-8 + @test F(-phi,m) ≈ -f atol=1e-8 + if phid == 90 + @test K(m) ≈ f atol=1e-8 + @test ellipke(m)[1] ≈ f atol=1e-8 + end + end + end + end +end + +@testset "Table 17.6: E(φ, m)" begin + # Abramowitz & Stegun, Table 17.6, p616 + # E(phi\alpha) + # α\φ 0 5 10 15 20 25 30 + table17p6_p616 = [ + 0. 0 0.08726_646 0.17453_293 0.26179_939 0.34906_585 0.43633_231 0.52359_878; + 2 0 0.08726_633 0.17453_185 0.26179_579 0.34905_742 0.43631_608 0.52357_119; + 4 0 0.08726_592 0.17452_864 0.26178_503 0.34903_218 0.43626_745 0.52348_856; + 6 0 0.08726_525 0.17452_330 0.26176_715 0.34899_025 0.43618_665 0.52335_123; + 8 0 0.08726_432 0.17451_587 0.26174_224 0.34893_181 0.43607_403 0.52315_981; + 10 0 0.08726_313 0.17450_636 0.26171_041 0.34885_714 0.43593_011 0.52291_511; + 12 0 0.08726_168 0.17449_485 0.26167_182 0.34876_657 0.43575_552 0.52261_821; + 14 0 0.08725_999 0.17448_137 0.26162_664 0.34866_055 0.43555_106 0.52227_039; + 16 0 0.08725_806 0.17446_599 0.26157_510 0.34853_954 0.43531_765 0.52187_317; + 18 0 0.08725_590 0.17444_879 0.26151_743 0.34840_412 0.43505_633 0.52142_828; + 20 0 0.08725_352 0.17442_985 0.26145_391 0.34825_492 0.43476_831 0.52093_770; + 22 0 0.08725_094 0.17440_926 0.26138_485 0.34809_262 0.43445_488 0.52040_357; + 24 0 0.08724_816 0.17438_712 0.26131_056 0.34791_800 0.43411_749 0.51982_827; + 26 0 0.08724_521 0.17436_353 0.26123_141 0.34773_187 0.43375_767 0.51921_436; + 28 0 0.08724_208 0.17433_862 0.26114_778 0.34753_510 0.43337_709 0.51856_461; + 30 0 0.08723_881 0.17431_250 0.26106_005 0.34732_863 0.43297_749 0.51788_193; + 32 0 0.08723_540 0.17428_529 0.26096_867 0.34711_342 0.43256_075 0.51716_944; + 34 0 0.08723_187 0.17425_714 0.26087_405 0.34689_050 0.43212_880 0.51643_040; + 36 0 0.08722_824 0.17422_817 0.26077_666 0.34666_093 0.43168_368 0.51566_820; + 38 0 0.08722_453 0.17419_852 0.26067_697 0.34642_580 0.43122_748 0.51488_638; + 40 0 0.08722_075 0.17416_835 0.26057_545 0.34618_625 0.43076_236 0.51408_862; + 42 0 0.08721_692 0.17413_779 0.26047_261 0.34594_343 0.43029_055 0.51327_866; + 44 0 0.08721_307 0.17410_700 0.26036_893 0.34569_850 0.42981_431 0.51246_037; + 46 0 0.08720_920 0.17407_613 0.26026_492 0.34545_266 0.42933_594 0.51163_767; + 48 0 0.08720_535 0.17404_531 0.26016_110 0.34520_710 0.42885_776 0.51081_454; + 50 0 0.08720_152 0.17401_472 0.26005_795 0.34496_302 0.42838_212 0.50999_501; + 52 0 0.08719_774 0.17398_449 0.25995_600 0.34472_162 0.42791_134 0.50918_310; + 54 0 0.08719_402 0.17395_477 0.25985_574 0.34448_409 0.42744_775 0.50838_287; + 56 0 0.08719_039 0.17392_571 0.25975_765 0.34425_159 0.42699_368 0.50759_831; + 58 0 0.08718_686 0.17389_745 0.25966_224 0.34402_529 0.42655_138 0.50683_341; + 60 0 0.08718_345 0.17387_013 0.25956_996 0.34380_631 0.42612_308 0.50609_207; + 62 0 0.08718_017 0.17384_388 0.25948_126 0.34359_575 0.42571_097 0.50537_811; + 64 0 0.08717_704 0.17381_883 0.25939_660 0.34339_465 0.42531_712 0.50469_523; + 66 0 0.08717_408 0.17379_511 0.25931_640 0.34320_404 0.42494_358 0.50404_700; + 68 0 0.08717_130 0.17377_283 0.25924_104 0.34302_487 0.42459_224 0.50343_686; + 70 0 0.08716_871 0.17375_210 0.25917_090 0.34285_805 0.42426_495 0.50286_804; + 72 0 0.08716_633 0.17373_302 0.25910_634 0.34270_443 0.42396_339 0.50234_359; + 74 0 0.08716_416 0.17371_568 0.25904_767 0.34256_478 0.42368_913 0.50186_633; + 76 0 0.08716_223 0.17370_018 0.25899_519 0.34243_984 0.42344_363 0.50143_886; + 78 0 0.08716_053 0.17368_659 0.25894_917 0.34233_022 0.42322_817 0.50106_351; + 80 0 0.08715_909 0.17367_498 0.25890_983 0.34223_650 0.42304_389 0.50074_232; + 82 0 0.08715_789 0.17366_539 0.25887_737 0.34215_915 0.42289_175 0.50047_707; + 84 0 0.08715_695 0.17365_789 0.25885_195 0.34209_857 0.42277_258 0.50026_923; + 86 0 0.08715_628 0.17365_250 0.25883_370 0.34205_507 0.42268_700 0.50011_993; + 88 0 0.08715_588 0.17364_926 0.25882_271 0.34202_889 0.42263_547 0.50003_003; + 90 0 0.08715_574 0.17364_818 0.25881_905 0.34202_014 0.42261_826 0.50000_000; + ] + + # Abramowitz & Stegun, Table 17.6, p617 + # E(phi\alpha) + # α\φ 35 40 45 50 55 60 + table17p6_p617 = [ + 0. 0.61086_524 0.69813_170 0.78539_816 0.87266_463 0.95993_109 1.04719_755; + 2 0.61082_230 0.69806_905 0.78531_125 0.87254_883 0.95978_184 1.04701_051; + 4 0.61069_365 0.69788_136 0.78505_085 0.87220_183 0.95933_459 1.04644_996; + 6 0.61047_983 0.69756_935 0.78461_792 0.87162_487 0.95859_083 1.04551_764; + 8 0.61018_171 0.69713_427 0.78401_409 0.87081_998 0.95755_301 1.04421_646; + 10 0.60980_055 0.69657_784 0.78324_162 0.86979_001 0.95622_460 1.04255_047; + 12 0.60933_793 0.69590_226 0.78230_343 0.86853_863 0.95461_005 1.04052_491; + 14 0.60879_577 0.69511_023 0.78120_308 0.86707_031 0.95271_478 1.03814_615; + 16 0.60817_636 0.69420_492 0.77994_473 0.86539_034 0.95054_522 1.03542_177; + 18 0.60748_229 0.69318_999 0.77853_323 0.86350_481 0.94810_878 1.03236_049; + 20 0.60671_652 0.69206_954 0.77697_402 0.86142_062 0.94541_386 1.02897_221; + 22 0.60588_229 0.69084_814 0.77527_316 0.85914_545 0.94246_984 1.02526_804; + 24 0.60498_319 0.68953_083 0.77343_735 0.85668_781 0.93928_709 1.02126_023; + 26 0.60402_308 0.68812_308 0.77147_387 0.85405_695 0.93587_699 1.01696_224; + 28 0.60300_616 0.68663_077 0.76939_059 0.85126_295 0.93225_186 1.01238_873; + 30 0.60193_687 0.68506_023 0.76719_599 0.84831_663 0.92842_504 1.00755_556; + 32 0.60081_994 0.68341_817 0.76489_908 0.84522_958 0.92441_083 1.00247_977; + 34 0.59966_035 0.68171_170 0.76250_947 0.84201_414 0.92022_452 0.99717_966; + 36 0.59846_332 0.67994_830 0.76003_726 0.83868_340 0.91588_234 0.99167_469; + 38 0.59723_431 0.67813_578 0.75749_309 0.83525_115 0.91140_150 0.98598_560; + 40 0.59597_897 0.67628_229 0.75488_809 0.83173_189 0.90680_017 0.98013_430; + 42 0.59470_312 0.67439_630 0.75223_383 0.82814_080 0.90209_742 0.97414_397; + 44 0.59341_278 0.67248_651 0.74954_234 0.82449_369 0.89731_325 0.96803_899; + 46 0.59211_406 0.67056_191 0.74682_605 0.82080_700 0.89246_858 0.96184_497; + 48 0.59081_324 0.66863_167 0.74409_773 0.81709_775 0.88758_513 0.95558_873; + 50 0.58951_664 0.66670_515 0.74137_047 0.81338_346 0.88268_551 0.94929_830; + 52 0.58823_065 0.66479_183 0.73865_766 0.80968_217 0.87779_305 0.94300_285; + 54 0.58696_171 0.66290_130 0.73597_286 0.80601_230 0.87293_184 0.93673_272; + 56 0.58571_622 0.66104_317 0.73332_979 0.80239_262 0.86812_660 0.93051_931; + 58 0.58450_056 0.65922_707 0.73074_229 0.79884_217 0.86340_261 0.92439_505; + 60 0.58332_103 0.65746_255 0.72822_416 0.79538_015 0.85878_561 0.91839_329; + 62 0.58218_382 0.65575_905 0.72578_915 0.79202_582 0.85430_169 0.91254_821; + 64 0.58109_497 0.65412_585 0.72345_085 0.78879_839 0.84997_709 0.90689_460; + 66 0.58006_032 0.65257_197 0.72122_260 0.78571_685 0.84583_811 0.90146_778; + 68 0.57908_549 0.65110_612 0.71911_737 0.78279_987 0.84191_082 0.89630_323; + 70 0.57817_584 0.64973_667 0.71714_767 0.78006_562 0.83822_090 0.89143_642; + 72 0.57733_641 0.64847_154 0.71532_545 0.77753_157 0.83479_335 0.88690_237; + 74 0.57657_189 0.64731_812 0.71366_196 0.77521_434 0.83165_223 0.88273_530; + 76 0.57588_663 0.64628_328 0.71216_766 0.77312_952 0.82882_031 0.87896_810; + 78 0.57528_450 0.64537_322 0.71085_210 0.77129_143 0.82631_879 0.87563_185; + 80 0.57476_897 0.64459_347 0.70972_381 0.76971_298 0.82416_694 0.87275_520; + 82 0.57434_302 0.64394_879 0.70879_019 0.76840_544 0.82238_177 0.87036_381; + 84 0.57400_912 0.64344_316 0.70805_745 0.76737_830 0.82097_770 0.86847_970; + 86 0.57376_921 0.64307_973 0.70753_050 0.76663_912 0.81996_631 0.86712_068; + 88 0.57362_470 0.64286_075 0.70721_289 0.76619_339 0.81935_604 0.86629_990; + 90 0.57357_644 0.64278_761 0.70710_678 0.76604_444 0.81915_204 0.86602_540; + ] + + # Abramowitz & Stegun, Table 17.6, p618 + # E(phi\alpha) + # 65 70 75 80 85 90 + table17p6_p618 = [ + 1.13446_401 1.22173_048 1.30899_694 1.39626_340 1.48352_986 1.57079_633; + 1.13423_517 1.22145_628 1.30867_442 1.39589_024 1.48310_448 1.57031_792; + 1.13354_929 1.22063_443 1.30770_767 1.39477_165 1.48182_929 1.56888_372; + 1.13240_837 1.21926_717 1.30609_916 1.39291_030 1.47970_717 1.56649_679; + 1.13081_573 1.21735_820 1.30385_297 1.39031_062 1.47674_288 1.56316_223; + 1.12877_602 1.21491_274 1.30097_484 1.38697_886 1.47294_312 1.55888_720; + 1.12629_522 1.21193_748 1.29747_215 1.38292_302 1.46831_652 1.55368_089; + 1.12338_066 1.20844_065 1.29335_393 1.37815_292 1.46287_363 1.54755_458; + 1.12004_099 1.20443_195 1.28863_089 1.37268_017 1.45662_693 1.54052_157; + 1.11628_624 1.19992_262 1.28331_541 1.36651_823 1.44959_085 1.53259_729; + 1.11212_778 1.19492_542 1.27742_153 1.35968_233 1.44178_179 1.52379_921; + 1.10757_834 1.18945_465 1.27096_502 1.35218_961 1.43321_813 1.51414_692; + 1.10265_204 1.18352_618 1.26396_337 1.34405_903 1.42392_023 1.50366_214; + 1.09736_439 1.17715_743 1.25643_578 1.33531_146 1.41391_049 1.49236_871; + 1.09173_228 1.17036_745 1.24840_326 1.32596_967 1.40321_335 1.48029_266; + 1.08577_404 1.16317_686 1.23988_858 1.31605_841 1.39185_532 1.46746_221; + 1.07950_942 1.15560_796 1.23091_635 1.30560_436 1.37986_503 1.45390_780; + 1.07295_961 1.14768_469 1.22151_305 1.29463_629 1.36727_328 1.43966_215; + 1.06614_728 1.13943_273 1.21170_705 1.28318_499 1.35411_306 1.42476_031; + 1.05909_660 1.13087_946 1.20152_870 1.27128_343 1.34041_965 1.40923_972; + 1.05183_322 1.12205_408 1.19101_036 1.25896_675 1.32623_066 1.39314_025; + 1.04438_435 1.11298_760 1.18018_648 1.24627_240 1.31158_614 1.37650_433; + 1.03677_875 1.10371_291 1.16909_366 1.23324_019 1.29652_865 1.35937_700; + 1.02904_677 1.09426_484 1.15777_077 1.21991_241 1.28110_340 1.34180_606; + 1.02122_034 1.08468_023 1.14625_899 1.20633_398 1.26535_837 1.32384_218; + 1.01333_305 1.07499_796 1.13460_200 1.19255_255 1.24934_449 1.30553_909; + 1.00542_010 1.06525_908 1.12284_604 1.17861_873 1.23311_580 1.28695_374; + 0.99751_835 1.05550_682 1.11104_010 1.16458_621 1.21672_971 1.26814_653; + 0.98966_632 1.04578_671 1.09923_604 1.15051_210 1.20024_724 1.24918_162; + 0.98190_414 1.03614_663 1.08748_883 1.13645_710 1.18373_339 1.23012_722; + 0.97427_354 1.02663_689 1.07585_669 1.12248_590 1.16725_747 1.21105_603; + 0.96681_780 1.01731_023 1.06440_132 1.10866_752 1.15089_364 1.19204_568; + 0.95958_158 1.00822_192 1.05318_814 1.09507_580 1.13472_145 1.17317_938; + 0.95261_084 0.99942_966 1.04228_653 1.08178_986 1.11882_658 1.15454_668; + 0.94595_256 0.99099_354 1.03176_998 1.06889_476 1.10330_172 1.13624_437; + 0.93965_447 0.98297_583 1.02171_634 1.05648_221 1.08824_773 1.11837_774; + 0.93376_462 0.97544_068 1.01220_781 1.04465_133 1.07377_505 1.10106_217; + 0.92833_088 0.96845_360 1.00333_091 1.03350_951 1.06000_556 1.08442_522; + 0.92340_024 0.96208_074 0.99517_606 1.02317_331 1.04707_504 1.06860_953; + 0.91901_802 0.95638_776 0.98783_670 1.01376_904 1.03513_640 1.05377_692; + 0.91522_691 0.95143_847 0.98140_781 1.00543_295 1.02436_393 1.04011_440; + 0.91206_588 0.94729_297 0.97598_331 0.99831_000 1.01495_896 1.02784_362; + 0.90956_905 0.94400_544 0.97165_228 0.99255_019 1.00715_650 1.01723_692; + 0.90776_445 0.94162_171 0.96849_392 0.98830_025 1.00123_026 1.00864_796; + 0.90667_305 0.94017_677 0.96657_142 0.98568_915 0.99748_392 1.00258_409; + 0.90630_779 0.93969_262 0.96592_583 0.98480_775 0.99619_470 1.00000_000; + ] + + table17p6 = [table17p6_p616 table17p6_p617[:,2:end] table17p6_p618] + + for i = 1:size(table17p6,1) + alpha = table17p6[i,1] + m = sind(alpha)^2 + for j = 2:size(table17p6,2) + phid = 5*(j-2) + phi = deg2rad(phid) + e = table17p6[i,j] + @test E(phi,m) ≈ e atol=1e-8 + @test E(-phi,m) ≈ -e atol=1e-8 + if phid == 90 + @test E(m) ≈ e atol=1e-8 + @test ellipke(m)[2] ≈ e atol=1e-8 + end + end + end + + @test E(0,0.5) ≈ 0 + @test E(0.5,0) ≈ 0.5 + @test E(0.5,1) ≈ sin(0.5) + + @test F(0,0.5) ≈ 0 + @test F(0.5,0) ≈ 0.5 + @test F(0.5,1) ≈ log(sec(0.5) + tan(0.5)) + + # values from Abramowitz & Stegun, Table 17.2 (p610-611) + m20 = sind(20)^2 + K20 = 1.62002_58991_24204 + E20 = 1.52379_92052_59774 + @test K(m20) ≈ K20 atol=1e-15 + @test E(m20) ≈ E20 atol=1e-15 + + # values from Abramowitz & Stegun, Table 17.5 (p613-615) + # values from Abramowitz & Stegun, Table 17.6 (p616-618) + E2020 = 0.34825_492 + F2020 = 0.34988_016 + for i = -2:2 + @test E(deg2rad(20+180i), m20) ≈ E2020 + 2i*E20 atol=1e-8 + @test F(deg2rad(20+180i), m20) ≈ F2020 + 2i*K20 atol=1e-8 + end +end + + +@testset "Table 17.9: Π(n; φ, m)" begin + # Abramowitz & Stegun, Table 17.9, p625-6 + # Π(n;phi\alpha) + # n α\φ 0 15 30 45 60 75 90 + table17p9 = [ + 0.0 0 0 0.26180 0.52360 0.78540 1.04720 1.30900 1.57080; + 0.0 15 0 0.26200 0.52513 0.79025 1.05774 1.32733 1.59814; + 0.0 30 0 0.26254 0.52943 0.80437 1.08955 1.38457 1.68575; + 0.0 45 0 0.26330 0.53562 0.82602 1.14243 1.48788 1.85407; + 0.0 60 0 0.26406 0.54223 0.85122 1.21260 1.64918 2.15651; + 0.0 75 0 0.26463 0.54736 0.87270 1.28371 1.87145 2.76806; + 0.0 90 0 0.26484 0.54931 0.88137 1.31696 2.02759 Inf; + 0.1 0 0 0.26239 0.52820 0.80013 1.07949 1.36560 1.65576; + 0.1 15 0 0.26259 0.52975 0.80514 1.09058 1.38520 1.68536; + 0.1 30 0 0.26314 0.53412 0.81972 1.12405 1.44649 1.78030; + 0.1 45 0 0.26390 0.54041 0.84210 1.17980 1.55739 1.96326; + 0.1 60 0 0.26467 0.54712 0.86817 1.25393 1.73121 2.29355; + 0.1 75 0 0.26524 0.55234 0.89040 1.32926 1.97204 2.96601; + 0.1 90 0 0.26545 0.55431 0.89939 1.36454 2.14201 Inf; + 0.2 0 0 0.26299 0.53294 0.81586 1.11534 1.43078 1.75620; + 0.2 15 0 0.26319 0.53452 0.82104 1.12705 1.45187 1.78850; + 0.2 30 0 0.26374 0.53896 0.83612 1.16241 1.51792 1.89229; + 0.2 45 0 0.26450 0.54535 0.85928 1.22139 1.63775 2.09296; + 0.2 60 0 0.26527 0.55217 0.88629 1.30003 1.82643 2.45715; + 0.2 75 0 0.26585 0.55747 0.90934 1.38016 2.08942 3.20448; + 0.2 90 0 0.26606 0.55948 0.91867 1.41777 2.27604 Inf; + 0.3 0 0 0.26359 0.53784 0.83271 1.15551 1.50701 1.87746; + 0.3 15 0 0.26379 0.53945 0.83808 1.16791 1.52988 1.91309; + 0.3 30 0 0.26434 0.54396 0.85370 1.20543 1.60161 2.02779; + 0.3 45 0 0.26511 0.55046 0.87771 1.26812 1.73217 2.25038; + 0.3 60 0 0.26588 0.55739 0.90574 1.35193 1.93879 2.65684; + 0.3 75 0 0.26646 0.56278 0.92969 1.43759 2.22876 3.49853; + 0.3 90 0 0.26667 0.56483 0.93938 1.47789 2.43581 Inf; + 0.4 0 0 0.26420 0.54291 0.85084 1.20098 1.59794 2.02789; + 0.4 15 0 0.26440 0.54454 0.85641 1.21419 1.62298 2.06774; + 0.4 30 0 0.26495 0.54912 0.87262 1.25419 1.70165 2.19629; + 0.4 45 0 0.26572 0.55573 0.89756 1.32117 1.84537 2.44683; + 0.4 60 0 0.26650 0.56278 0.92670 1.41098 2.07413 2.90761; + 0.4 75 0 0.26708 0.56827 0.95162 1.50309 2.39775 3.87214; + 0.4 90 0 0.26729 0.57035 0.96171 1.54653 2.63052 Inf; + 0.5 0 0 0.26481 0.54814 0.87042 1.25310 1.70919 2.22144; + 0.5 15 0 0.26501 0.54980 0.87621 1.26726 1.73695 2.26685; + 0.5 30 0 0.26557 0.55447 0.89307 1.31017 1.82433 2.41367; + 0.5 45 0 0.26634 0.56119 0.91902 1.38218 1.98464 2.70129; + 0.5 60 0 0.26712 0.56837 0.94939 1.47906 2.24155 3.23477; + 0.5 75 0 0.26770 0.57394 0.97538 1.57881 2.60846 4.36620; + 0.5 90 0 0.26792 0.57606 0.98591 1.62599 2.87468 Inf; + 0.6 0 0 0.26543 0.55357 0.89167 1.31379 1.85002 2.48365; + 0.6 15 0 0.26563 0.55525 0.89770 1.32907 1.88131 2.53677; + 0.6 30 0 0.26619 0.56000 0.91527 1.37544 1.98005 2.70905; + 0.6 45 0 0.26696 0.56684 0.94235 1.45347 2.16210 3.04862; + 0.6 60 0 0.26775 0.57414 0.97406 1.55884 2.45623 3.68509; + 0.6 75 0 0.26833 0.57982 1.00123 1.66780 2.88113 5.05734; + 0.6 90 0 0.26855 0.58198 1.01225 1.71951 3.19278 Inf; + 0.7 0 0 0.26605 0.55918 0.91487 1.38587 2.03720 2.86787; + 0.7 15 0 0.26625 0.56090 0.92116 1.40251 2.07333 2.93263; + 0.7 30 0 0.26681 0.56573 0.93952 1.45309 2.18765 3.14339; + 0.7 45 0 0.26759 0.57270 0.96784 1.53846 2.39973 3.56210; + 0.7 60 0 0.26838 0.58014 1.00104 1.65425 2.74586 4.35751; + 0.7 75 0 0.26897 0.58592 1.02954 1.77459 3.25315 6.11030; + 0.7 90 0 0.26918 0.58812 1.04110 1.83192 3.63042 Inf; + 0.8 0 0 0.26668 0.56501 0.94034 1.47370 2.30538 3.51240; + 0.8 15 0 0.26688 0.56676 0.94694 1.49205 2.34868 3.59733; + 0.8 30 0 0.26745 0.57168 0.96618 1.54790 2.48618 3.87507; + 0.8 45 0 0.26823 0.57877 0.99588 1.64250 2.74328 4.43274; + 0.8 60 0 0.26902 0.58635 1.03076 1.77145 3.16844 5.51206; + 0.8 75 0 0.26961 0.59225 1.06073 1.90629 3.80370 7.96669; + 0.8 90 0 0.26982 0.59449 1.07290 1.97080 4.28518 Inf; + 0.9 0 0 0.26731 0.57106 0.96853 1.58459 2.74439 4.96729; + 0.9 15 0 0.26752 0.57284 0.97547 1.60515 2.79990 5.09958; + 0.9 30 0 0.26808 0.57785 0.99569 1.66788 2.97710 5.53551; + 0.9 45 0 0.26887 0.58508 1.02695 1.77453 3.31210 6.42557; + 0.9 60 0 0.26966 0.59281 1.06372 1.92081 3.87661 8.20086; + 0.9 75 0 0.27025 0.59882 1.09535 2.07487 4.74432 12.46407; + 0.9 90 0 0.27047 0.60110 1.10821 2.14899 5.42125 Inf; + 1.0 0 0 0.26795 0.57735 1.00000 1.73205 3.73205 Inf; + 1.0 15 0 0.26816 0.57916 1.00731 1.75565 3.81655 Inf; + 1.0 30 0 0.26872 0.58428 1.02866 1.82781 4.08864 Inf; + 1.0 45 0 0.26951 0.59165 1.06170 1.95114 4.61280 Inf; + 1.0 60 0 0.27031 0.59953 1.10060 2.12160 5.52554 Inf; + 1.0 75 0 0.27090 0.60566 1.13414 2.30276 7.00372 Inf; + 1.0 90 0 0.27112 0.60799 1.14779 2.39053 8.22356 Inf; + ] + + for i = 1:size(table17p9,1) + n = table17p9[i,1] + alpha = table17p9[i,2] + m = sind(alpha)^2 + for j = 3:size(table17p9,2) + phi = 15.0*(j-3) + x = table17p9[i,j] + y = Pi(n, deg2rad(phi), m) + if x == Inf + @test x == y + else + @test x ≈ y atol=1e-5*max(x,1.) + end + end + end +end + +include("autodiff.jl")