Code revision

Project.toml

@@ -4,14 +4,15 @@ authors = ["Kiran Shila\n <[email protected]>"]
 version = "0.1.0"
 
 [deps]
+IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 julia = "1"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","QuadGK"]
+test = ["Test", "QuadGK"]

src/FresnelIntegrals.jl

@@ -1,27 +1,97 @@
 module FresnelIntegrals
 using SpecialFunctions
-export fresnelc
-export fresnels
-export fresnel
+using IrrationalConstants: sqrtπ
+export fresnelc, fresnels, fresnel
 
 """
-        fresnelc(z)
-Calculates the Fresnel cosine integral for the number z for
-    ``C(z) = \\int_{0}^{z} \\cos{\\left(\\frac{\\pi t^2}{2}\\right)}dt``
+    fresnelc(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``.
 """
-fresnelc(z::Number) = 0.25*(1-1im)*(1im*erf(0.5*(1-1im)*z*√(π)) + erf(0.5*(1+1im)*z*√(π)))
+function fresnelc(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 fresnelc(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
 
 """
-        fresnels(z)
-Calculates the Fresnel sine integral for the number z for
-    ``S(z) = \\int_{0}^{z} \\sin{\\left(\\frac{\\pi t^2}{2}\\right)}dt``
+    fresnels(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``.
 """
-fresnels(z::Number) = 0.25*(1+1im)*(-1im*erf(0.5*(1-1im)*z*√(π)) + erf(0.5*(1+1im)*z*√(π)))
+function fresnels(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 fresnels(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
 
 """
-        fresnel(z)
-Calculates the cosine and sine fresnel integrals
+    fresnel(z::Number)
+
+Calculate the normalized cosine and sine fresnel integrals.
+
+See also [`fresnels`](@ref), [`fresnelc`](@ref).
 """
-fresnel(z::Number) = (fresnelc(z),fresnels(z))
+function fresnel(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 fresnel(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
+
+
 
 end # module

test/runtests.jl

@@ -10,4 +10,20 @@ using QuadGK
     @test fresnels(z) ≈ quadgk(t->sin(π*t^2/2),0,z)[1]
     # Test just for code coverage 😄
     @test (fresnelc(z),fresnels(z)) == fresnel(z)
+
+    # Generate random real number
+    z = randn(Float64)
+    # Test by comparing to numeric solution
+    @test fresnelc(z) ≈ quadgk(t->cos(π*t^2/2),0,z)[1]
+    @test fresnels(z) ≈ quadgk(t->sin(π*t^2/2),0,z)[1]
+    # Test just for code coverage 😄
+    @test (fresnelc(z),fresnels(z)) == fresnel(z)
+
+    # Precise values come from WolframAlpha calculator
+    # One could add more decimals and more tests if needed
+
+    @test fresnels(1.)      ≈ 0.4382591473903
+    @test fresnelc(1.)      ≈ 0.7798934003768
+    @test fresnels(sqrt(2)*im) ≈ -0.7139722140219*im
+    @test fresnelc(sqrt(2)*im) ≈ 0.5288915951112*im
 end