Changeset d17228 for qd-gamma.lisp


Ignore:
Timestamp:
03/18/11 02:14:25 (4 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
50fc64
Parents:
46e719
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.

File:
1 edited

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     
Note: See TracChangeset for help on using the changeset viewer.