Added Fresnel integrals #398

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"

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")

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
+

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

test/runtests.jl

@@ -22,6 +22,7 @@ checktol(err::Float64) = err ≤ 1e-13
 
 
 tests = [
+    "fresnel",
     "bessel",
     "beta_inc",
     "betanc",