Show
Ignore:
Timestamp:
03/17/11 19:14:25 (3 years ago)
Author:
Raymond Toy <toy.raymond@…>
Children:
50fc6412a3effd9bea8f3e481f7c35c3c911f856
Parents:
46e7193fb5b2fad35e67366ff375d9ebbb524c5c
git-committer:
Raymond Toy <toy.raymond@…> (03/17/11 19:14:25)
Message:

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.

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • qd-gamma.lisp

    r46e719 rd17228  
    293293            (cf-incomplete-gamma a z) 
    294294            (- (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)))))) 
    296298 
    297299(defun erf (z) 
     
    342344  (* (expt z (- v 1)) 
    343345     (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