diff --git a/Project.toml b/Project.toml index fad59a8..516b602 100644 --- a/Project.toml +++ b/Project.toml @@ -4,14 +4,15 @@ authors = ["Kiran Shila\n "] 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"] diff --git a/src/FresnelIntegrals.jl b/src/FresnelIntegrals.jl index c37f2e8..e6636b2 100644 --- a/src/FresnelIntegrals.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index d70e55c..a48b2b9 100644 --- a/test/runtests.jl +++ b/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