verifying the asymptotic expansion of Fresnel auxiliary functions

fresnel-aux.py

from numpy import pi, cos, sin


def f(z, s_z, c_z):
    return (0.5 - s_z) * cos(pi / 2 * z ** 2) - (0.5 - c_z) * sin(pi / 2 * z ** 2)


def f_asymp(z):
    zp2 = (pi * z ** 2) ** 2
    r = 1.0
    f = 1.0
    for k in range(1, 21):
        r *= -(4 * k - 1) * (4 * k - 3) / zp2
        f += r
    f /= pi * z
    return f


z0 = 5.0
# values from wolframalpha
c_z0 = 0.5636311887040122311021074044130139641207537623099921078616593412
s_z0 = 0.4991913819171168867519283804659916554084319970723881534101411151
print(f(z0, s_z0, c_z0))
print(f_asymp(z0))
# 0.06363118870401219
# 0.06363118870401226


z0 = 5.0 + 0.3j
# values from wolframalpha
c_z0 = (
    4.02556471559977327707822210694337432408377169490431616402495978
    + 0.332688495001379877552515347849647714872236656530633626755045188j
)
s_z0 = (
    0.167190623257434661649247894355254546820257708244723444212293547
    + 3.52500790048219018446948210309003170353997178015674671343643381j
)
print()
print(f(z0, s_z0, c_z0))
print(f_asymp(z0))
# (0.06340444527984346-0.003797046505276569j)
# (0.06340444527985087-0.0037970465052788557j)


z0 = 4.0 + 1.0j
# values from wolframalpha
c_z0 = (
    -10871.7039066376346912813725171464799373397451050087823687267604
    + 2519.85253846801037430128602658852786575926239432986235040331346j
)
s_z0 = (
    -2519.35253839923345767133409785512608007673244754274649087284911
    - 10872.2039063798973572036565705336849997516940364531952381900864j
)
print()
print(f(z0, s_z0, c_z0))
print(f_asymp(z0))
# (0.07486844062805176-0.018648505210876465j)
# (0.07486835668512495-0.018648524070775562j)


z0 = 4.0 - 1.0j
# values from wolframalpha
c_z0 = (
    -10871.7039066376346912813725171464799373397451050087823687267604
    - 2519.85253846801037430128602658852786575926239432986235040331346j
)
s_z0 = (
    -2519.35253839923345767133409785512608007673244754274649087284911
    + 10872.2039063798973572036565705336849997516940364531952381900864j
)
print("x")
print(f(z0, s_z0, c_z0))
print(f_asymp(z0))
# whoops! only 4 digits agree
# (0.07486844062805176+0.018648505210876465j)
# (0.07486835668512495+0.018648524070775562j)


z0 = 1.0 + 4.0j
# values from wolframalpha
c_z0 = (
    2519.852538468010374301286026588527865759262394329862350403313465
    - 10871.70390663763469128137251714647993733974510500878236872676043j
)
s_z0 = (
    10872.20390637989735720365657053368499975169403645319523819008643
    + 2519.352538399233457671334097855126080076732447542746490872849118j
)
print()
print(f(z0, s_z0, c_z0))
print(f_asymp(z0))
# whoops! only 4 digits agree
# (0.01864677667617798-0.07486653327941895j)
# (0.018648524070775562-0.07486835668512495j)


z0 = 0.1 + 5.0j
# values from wolframalpha
c_z0 = (
    0.000971856588047086154178458248331847140322418344278944596995681
    + 0.659704145460982834686734591751993165227356600071220819810856j
)
s_z0 = (
    -0.14650022466696344385063124801975259656544735691449836752404069
    - 0.49838903151412768700839683673807747651830440950230724938834477j
)
print()
print(f(z0, s_z0, c_z0))
print(f_asymp(z0))
# yikes!
# (0.1068292542731959-0.16590026661344703j)
# (0.0012696736573447598-0.06360591646598417j)