diff --git a/Project.toml b/Project.toml index 35dd1b93..bc8cabd4 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" OpenLibm_jll = "05823500-19ac-5b8b-9628-191a04bc5112" OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" [compat] ChainRulesCore = "0.9.44, 0.10, 1" diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index c4fe6120..1058e509 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -78,8 +78,12 @@ export expintx, sinint, cosint, - lbinomial - + lbinomial, + fresnelcos, + fresnelsin, + fresnelsincos + +include("fresnel.jl") include("bessel.jl") include("erf.jl") include("ellip.jl") diff --git a/src/fresnel.jl b/src/fresnel.jl new file mode 100644 index 00000000..d265930b --- /dev/null +++ b/src/fresnel.jl @@ -0,0 +1,92 @@ +# Code ported from https://github.com/kiranshila/FresnelIntegrals.jl + + +""" + fresnelcos(z::Number) + +Calculate the normalized Fresnel cosine integral +```math +C(z) = \\int_{0}^{z} \\cos{\\left(\\frac{\\pi t^2}{2}\\right)} \\, \\mathrm{d}t +``` +for the number ``z``. +""" +function fresnelcos(z::Number) + x = (z * sqrtπ) / 2 + re_x, im_x = reim(x) + a = (re_x + im_x) + (im_x - re_x) * im + b = (re_x - im_x) + (im_x + re_x) * im + re_erf_a, im_erf_a = reim(erf(a)) + re_erf_b, im_erf_b = reim(erf(b)) + re_y = (re_erf_a - im_erf_a + re_erf_b + im_erf_b) / 4 + im_y = (im_erf_a + re_erf_a - re_erf_b + im_erf_b) / 4 + y = re_y + im_y * im + return y +end +function fresnelcos(z::Real) + x = (z * sqrtπ) / 2 + a = x + x * im + re_erf_a, im_erf_a = reim(erf(a)) + y = (re_erf_a + im_erf_a) / 2 + return y +end + +""" + fresnelsin(z::Number) + +Calculate the normalized Fresnel sine integral +```math +S(z) = \\int_{0}^{z} \\sin{\\left(\\frac{\\pi t^2}{2}\\right)} \\, \\mathrm{d}t +``` +for the number ``z``. +""" +function fresnelsin(z::Number) + x = (z * sqrtπ) / 2 + re_x, im_x = reim(x) + a = (re_x + im_x) + (im_x - re_x) * im + b = (re_x - im_x) + (im_x + re_x) * im + re_erf_a, im_erf_a = reim(erf(a)) + re_erf_b, im_erf_b = reim(erf(b)) + re_y = (re_erf_a + im_erf_a + re_erf_b - im_erf_b) / 4 + im_y = (im_erf_a - re_erf_a + re_erf_b + im_erf_b) / 4 + y = re_y + im_y * im + return y +end +function fresnelsin(z::Real) + x = (z * sqrtπ) / 2 + a = x + x * im + re_erf_a, im_erf_a = reim(erf(a)) + y = (re_erf_a - im_erf_a) / 2 + return y +end + +""" + fresnelsincos(z::Number) + +Calculate the normalized cosine and sine fresnel integrals. + +See also [`fresnelsin`](@ref), [`fresnelcos`](@ref). +""" +function fresnelsincos(z::Number) + x = (z * sqrtπ) / 2 + re_x, im_x = reim(x) + a = (re_x + im_x) + (im_x - re_x) * im + b = (re_x - im_x) + (im_x + re_x) * im + re_erf_a, im_erf_a = reim(erf(a)) + re_erf_b, im_erf_b = reim(erf(b)) + re_y_sin = (re_erf_a + im_erf_a + re_erf_b - im_erf_b) / 4 + im_y_sin = (im_erf_a - re_erf_a + re_erf_b + im_erf_b) / 4 + re_y_cos = (re_erf_a - im_erf_a + re_erf_b + im_erf_b) / 4 + im_y_cos = (im_erf_a + re_erf_a - re_erf_b + im_erf_b) / 4 + y_sin = re_y_sin + im_y_sin * im + y_cos = re_y_cos + im_y_cos * im + return (y_cos, y_sin) +end +function fresnelsincos(z::Real) + x = (z * sqrtπ) / 2 + a = x + x * im + re_erf_a, im_erf_a = reim(erf(a)) + y_sin = (re_erf_a - im_erf_a) / 2 + y_cos = (re_erf_a + im_erf_a) / 2 + return (y_cos, y_sin) +end + diff --git a/test/fresnel.jl b/test/fresnel.jl new file mode 100644 index 00000000..87acb167 --- /dev/null +++ b/test/fresnel.jl @@ -0,0 +1,26 @@ +using QuadGK +@testset "fresnel" begin + # Generate random complex number + z = randn(ComplexF64) + # Test by comparing to numeric solution + @test fresnelcos(z) ≈ quadgk(t->cos(π*t^2/2),0,z)[1] + @test fresnelsin(z) ≈ quadgk(t->sin(π*t^2/2),0,z)[1] + # Test just for code coverage 😄 + @test (fresnelcos(z),fresnelsin(z)) == fresnelsincos(z) + + # Generate random real number + z = randn(Float64) + # Test by comparing to numeric solution + @test fresnelcos(z) ≈ quadgk(t->cos(π*t^2/2),0,z)[1] + @test fresnelsin(z) ≈ quadgk(t->sin(π*t^2/2),0,z)[1] + # Test just for code coverage 😄 + @test (fresnelcos(z),fresnelsin(z)) == fresnelsincos(z) + + # Precise values come from WolframAlpha calculator + # One could add more decimals and more tests if needed + + @test fresnelsin(1.) ≈ 0.4382591473903 + @test fresnelcos(1.) ≈ 0.7798934003768 + @test fresnelsin(sqrt(2)*im) ≈ -0.7139722140219*im + @test fresnelcos(sqrt(2)*im) ≈ 0.5288915951112*im +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index f4a30360..f81b1a2c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,7 @@ checktol(err::Float64) = err ≤ 1e-13 tests = [ + "fresnel", "bessel", "beta_inc", "betanc",