Commit d17228b86105b8a6ac652459b1100266ad0311b3

Authored by Raymond Toy
1 parent 46e7193f

Add Fresnel integrals; fix issue in incomplete-gamma for large args.

o INCOMPLETE-GAMMA was returning bad values for large (complex)
  arguments.  Fix this by using incomplete gamma tail function since
  the incomplete gamma function approaches gamma for large arguments.

o Implement Fresnel S and C functions.
Showing 1 changed file with 48 additions and 1 deletions
qd-gamma.lisp
... ... @@ -292,7 +292,9 @@
292 292 (if (< z (- a 1))
293 293 (cf-incomplete-gamma a z)
294 294 (- (gamma a) (cf-incomplete-gamma-tail a z)))
295   - (cf-incomplete-gamma a z))))
  295 + (if (< (abs z) (abs a))
  296 + (cf-incomplete-gamma a z)
  297 + (- (gamma a) (cf-incomplete-gamma-tail a z))))))
296 298  
297 299 (defun erf (z)
298 300 "Error function:
... ... @@ -341,3 +343,48 @@
341 343 ;; Wolfram gives E(v,z) = z^(v-1)*gamma_incomplete_tail(1-v,z)
342 344 (* (expt z (- v 1))
343 345 (incomplete-gamma-tail (- 1 v) z)))
  346 +
  347 +(defun fresnel-s (z)
  348 + "Fresnel S:
  349 +
  350 + S(z) = integrate(sin(%pi*t^2/2), t, 0, z) "
  351 + (let ((sqrt-pi (sqrt (float-pi z))))
  352 + (flet ((fs (z)
  353 + ;; Wolfram gives
  354 + ;;
  355 + ;; S(z) = (1+%i)/4*(erf(c*z) - %i*erf(conjugate(c)*z))
  356 + ;;
  357 + ;; where c = sqrt(%pi)/2*(1+%i).
  358 + (* #c(1/4 1/4)
  359 + (- (erf (* #c(1/2 1/2) sqrt-pi z))
  360 + (* #c(0 1)
  361 + (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
  362 + (if (realp z)
  363 + ;; FresnelS is real for a real argument. And it is odd.
  364 + (if (minusp z)
  365 + (- (realpart (fs (- z))))
  366 + (realpart (fs z)))
  367 + (fs z)))))
  368 +
  369 +(defun fresnel-c (z)
  370 + "Fresnel C:
  371 +
  372 + C(z) = integrate(cos(%pi*t^2/2), t, 0, z) "
  373 + (let ((sqrt-pi (sqrt (float-pi z))))
  374 + (flet ((fs (z)
  375 + ;; Wolfram gives
  376 + ;;
  377 + ;; C(z) = (1-%i)/4*(erf(c*z) + %i*erf(conjugate(c)*z))
  378 + ;;
  379 + ;; where c = sqrt(%pi)/2*(1+%i).
  380 + (* #c(1/4 -1/4)
  381 + (+ (erf (* #c(1/2 1/2) sqrt-pi z))
  382 + (* #c(0 1)
  383 + (erf (* #c(1/2 -1/2) sqrt-pi z)))))))
  384 + (if (realp z)
  385 + ;; FresnelS is real for a real argument. And it is odd.
  386 + (if (minusp z)
  387 + (- (realpart (fs (- z))))
  388 + (realpart (fs z)))
  389 + (fs z)))))
  390 +
344 391 \ No newline at end of file
... ...