It turns out that SciPy's Fresnel values are wrong for complex arguments and large enough absolute value. I'm trying to fix that.
The implementation is based on Zhang/Jin, Computation of special functions, which in turn is based on Abramowitz/Stegun, Handbook of Mathematical Functions. There we find for the Fresnel S integral (7.3.10)
$$
S(z) = \frac{1}{2} – f(z) \cos\left(\frac{\pi}{2} z^2\right) – g(z) \sin\left(\frac{\pi}{2} z^2\right)
$$
for all $z$ with the auxiliary functions (7.3.5), (7.3.6)
$$
\begin{split}
f(z) &= \left[\frac{1}{2} – S(z) \right] \cos\left(\frac{\pi}{2} z^2\right) – \left[\frac{1}{2} – C(z) \right] \sin\left(\frac{\pi}{2} z^2\right),\\
g(z) &= \left[\frac{1}{2} – C(z) \right] \cos\left(\frac{\pi}{2} z^2\right) + \left[\frac{1}{2} – S(z) \right] \sin\left(\frac{\pi}{2} z^2\right).
\end{split}
$$
Computation of $S$ for large values is done via the asymptotic expansion of $f$ (7.3.27)
$$
\DeclareMathOperator\arg{arg}
\pi z f(z)\sim 1 + \sum_{m=1}^\infty (-1)^m \frac{1\cdot 3\cdot\cdots \cdot (4m-1)}{(\pi z^2)^{2m}}, \quad z\to\infty, |\arg(z)|<\frac{\pi}{2}.
$$
This is where I have problems understanding the approximation.
Consider $f(iz)$; for the asymptotic $f_a$, we have $f_a(iz) = -if_a(z)$. However no such thing is true for $f$ itself. Numerical computation via the representations
$$
\DeclareMathOperator\erf{erf}
\begin{split}
S(z) &= \frac{1 + i}{4} \left[\erf\left(\frac{1 + i}{2} \sqrt{\pi} z\right) – i \erf\left(\frac{1 – i}{2} \sqrt{\pi} z\right)\right]\\
C(z) &= \frac{1 – i}{4} \left[\erf\left(\frac{1 + i}{2} \sqrt{\pi} z\right) + i \erf\left(\frac{1 – i}{2} \sqrt{\pi} z\right)\right]
\end{split}
$$
shows that the approximation incorrect for everything off of the real axis, but the issue here could be numerical instability in the $\erf$ represenation too.
My current guess is that the above infinite sum is valid only in the area $|\arg{z}|<\pi/4$, which would already change how we compute the Fresnel integral values significantly.
To finish things off, here's a cplot of $f$ (for smaller $|z|$):
Best Answer
I will consider the function $\operatorname{f}(z)$, the treatment of $\operatorname{g}(z)$ is analogous. By http://dlmf.nist.gov/7.12.ii, we have $$ \operatorname{f}(z) = \frac{1}{{\pi z}}\sum\limits_{m = 0}^{N - 1} {( - 1)^m \left( {\frac{1}{2}} \right)_{2m} \frac{1}{{(\pi z^2 /2)^{2m} }}} + R_N^{(\operatorname{f})} (z) $$ where \begin{align*} R_N^{(\operatorname{f})} (z) & = \frac{{( - 1)^N }}{{\pi \sqrt 2 }}\int_0^{ + \infty } {\frac{{e^{ - \pi z^2 t/2} t^{2N - 1/2} }}{{1 + t^2 }}dt} \\ & = \frac{1}{{\pi z}}( - 1)^N \left( {\frac{1}{2}} \right)_{2N} \frac{1}{{(\pi z^2 /2)^{2N} }}\frac{1}{{\Gamma \left( {2N + \frac{1}{2}} \right)}}\int_0^{ + \infty } {\frac{{e^{ - s} s^{2N - 1/2} }}{{1 + s^2 /(\pi z^2 /2)^2 }}ds} \\ & = \frac{1}{{\pi z}}( - 1)^N \left( {\frac{1}{2}} \right)_{2N} \frac{1}{{(\pi z^2 /2)^{2N} }}\Pi _{2N + 1/2} (\pi z^2 /2), \end{align*} provided that $|\arg z|<\frac{\pi}{4}$ and $N\geq 0$. Here $\Pi_p(w)$ denotes one of Dingle's basic terminants: $$ \Pi _p (w) = \frac{1}{{\Gamma (p)}}\int_0^{ + \infty } {\frac{{e^{ - s} s^{p - 1} }}{{1 + (s/w)^2 }}ds} $$ for $|\arg w|<\frac{\pi}{2}$ and by analytic continuation elswhere. Using the expression for $R_N^{(\operatorname{f})} (z)$ in terms of this terminant, we can extend $R_N^{(\operatorname{f})} (z)$ to the universal covering of $\mathbb C \setminus \left\{ 0\right\}$. Now employing the estimates for the basic terminant established in https://doi.org/10.1007/s10440-017-0099-0, we obtain the bound \begin{align*} \left| {R_N^{(\operatorname{f})} (z)} \right| \le &\; \left| {\frac{1}{{\pi z}}( - 1)^N \left( {\frac{1}{2}} \right)_{2N} \frac{1}{{(\pi z^2 /2)^{2N} }}} \right| \\ & \times \begin{cases} 1 & \text{ if } \; \left|\arg z\right| \leq \frac{\pi}{8}, \\ \min\!\Big(\left|\csc ( 4\arg z)\right|,1 + \cfrac{1}{2}\chi(2N+1/2)\Big) & \text{ if } \; \frac{\pi}{8} < \left|\arg z\right| \leq \frac{\pi}{4}, \\ \cfrac{\sqrt {2\pi (2N + 1/2)} }{2\left| {\sin (2\arg z)} \right|^{2N+1/2} } + 1 + \cfrac{1}{2}\chi (2N +1/2) & \text{ if } \; \frac{\pi}{4} < \left|\arg z\right| < \frac{\pi}{2}. \end{cases} \end{align*} Here $\chi(p) =\sqrt{\pi}\Gamma(p/2+1)/\Gamma(p/2+1/2)$ for $p>0$. It is seen that the asymptotic expansion of $\operatorname{f}(z)$ is valid in every closed sub-sector of $|\arg z|<\frac{\pi}{2}$ in the sense of Poincaré. The Stokes lines are $\arg z =\pm \frac{\pi}{4}$ where terms are swiched on that remain exponentially small compared to any negative power of $z$ as long as we stay away from the rays $\arg z =\pm \frac{\pi}{2}$ (the anti-Stokes lines).
To obtain a better result and reveal the exponentially small terms, we can use the functional equation $$ \Pi _p (w) = \pm \pi i\frac{{e^{ \mp \frac{\pi }{2}ip} }}{{\Gamma (p)}}w^p e^{ \pm iw} + \Pi _p (we^{ \mp \pi i} ) $$ where $p>0$ and $w$ is any element of the universal covering of $\mathbb C \setminus \left\{ 0\right\}$ (the Riemann surface of the logarithm). With this functional equation, we find, after some algebra, $$ R_N^{(\operatorname{f})} (ze^{ \mp \frac{\pi }{2}i} ) = \pm iR_N^{(\operatorname{f})} (z) + \frac{{1 \mp i}}{2}e^{ \pm \frac{\pi }{2}iz^2 } . $$ This result is valid for all $N\geq 0$ and $z$ on the universal covering of $\mathbb C \setminus \left\{ 0\right\}$. You can see that if we omit the second term (which is exponentially small when $ 0 < \pm \arg z < \frac{\pi}{2}$) then, with $N=0$, we get the false result $$ \operatorname{f}(ze^{ \mp \frac{\pi }{2}i} ) = \pm i\operatorname{f}(z). $$ See also http://dlmf.nist.gov/7.4.
In summary, use $$ \operatorname{f}(z) \sim \frac{1}{{\pi z}}\sum\limits_{m = 0}^\infty {( - 1)^m \left( {\frac{1}{2}} \right)_{2m} \frac{1}{{(\pi z^2 /2)^{2m} }}} $$ when $\left| {\arg z} \right| \le \frac{\pi }{4}$, and use $$ \operatorname{f}(z) \sim \frac{{1 \pm i}}{2}e^{ \pm \frac{\pi }{2}iz^2 } + \frac{1}{{\pi z}}\sum\limits_{m = 0}^\infty {( - 1)^m \left( {\frac{1}{2}} \right)_{2m} \frac{1}{{(\pi z^2 /2)^{2m} }}} $$ when $\frac{\pi }{4} < \pm \arg z \le \frac{{3\pi }}{4}$. Of couse by relaying on the symmetry relation http://dlmf.nist.gov/7.4.E8, we can assume $|\arg z|\leq \frac{\pi}{2}$.