share/contrib/fresnel/fresnel.mac

   1 /* This function computes the fresnel F and G functions with dfloat accuracy
   2 over -inf,inf.  It divides the real line into 2 regions.  The first (0,4) uses
   3 a division of 1/100, and a taylor series.  The second region (4,inf)
   4 uses the asymptotic series.  This routine could be improved by using
   5 lazy evaluation to increase speed in the initialization.
   6 */
   7
   8 load("fresnel_linear_values.mac")$
   9 define_variable(FRESNEL_ASYMPTOTIC,4.0d0,any_check,
  10    "Point to switch to asymptotic series")$
  11 define_variable(FRESNEL_DIVISION,100,any_check,
  12    "Number of divisions")$
  13 define_variable(FRESNEL_X,0,any_check, "Array constant")$
  14 define_variable(FRESNEL_F,1,any_check, "Array constant")$
  15 define_variable(FRESNEL_G,2,any_check, "Array constant")$
  16 define_variable(FRESNEL_A,3,any_check, "Array constant")$
  17 array(fresnel_a,3,4*FRESNEL_DIVISION+1)$
  18 define_variable(fresnel_Ftaylor, 8.33333333333334d-3 *
  19    (u *
  20      (u *
  21        (u *
  22          (u *
  23            (u * x0 *
  24              (x0 *
  25               (x0 * (9.74090910340024d+2 * f0 - 3.0601968478528d+2 * g0 * x0^2)
  26                - 2.79056490122698d+2)
  27               + 4.65094150204496d+2 * g0)
  28              + x0^2 *
  29                  (x0 * (4.87045455170012d+2 * f0 * x0 - 1.55031383401498d+2)
  30                   + 9.30188300408994d+2 * g0)
  31             - 1.4804406601634d+2 * f0)
  32            + x0 * (6.20125533605996d+2 * g0 * x0^2 - 5.9217626406536d+2 * f0)
  33           + 1.25663706143592d+2)
  34          + x0 * (1.88495559215388d+2 - 5.9217626406536d+2 * f0 * x0)
  35         - 1.88495559215388d+2 * g0)
  36       - 3.76991118430774d+2 * g0 * x0)
  37     + 1.2d+2 * f0),any_check,"Taylor series for F about x0")$
  38 define_variable(fresnel_Gtaylor, 8.33333333333334d-3 *
  39    (u *
  40     (u *
  41      (u *
  42       (u *
  43        (x0 *
  44         (x0 *
  45          (4.87045455170012d+2 * g0 * x0^2 - 9.30188300408994d+2 * f0)
  46          + 2.46740110027234d+2)
  47         + u *
  48           (x0 *
  49            (x0^2 *
  50             (x0 *
  51              (3.0601968478528d+2 * f0 * x0 - 9.74090910340024d+1)
  52              + 9.74090910340024d+2 * g0)
  53             - 4.65094150204496d+2 * f0)
  54            + 7.89568352087148d+1)
  55         - 1.4804406601634d+2 * g0)
  56        + x0 * (x0 * (1.97392088021786d+2 - 6.20125533605996d+2 * f0 * x0)
  57        - 5.9217626406536d+2 * g0))
  58       - 5.9217626406536d+2 * g0 * x0^2 + 1.88495559215388d+2 * f0)
  59      + 3.76991118430774d+2 * f0 * x0 - 1.2d+2)
  60     + 1.2d+2 * g0),any_check,"Taylor series for G about x0")$
  61
  62 fresnel_InitLinear():=block([i:0],
  63    for l in fresnel_linear_values do (
  64       fresnel_a[FRESNEL_X,i]:first(l),
  65       block([x0:first(l),f0:second(l),g0:third(l)],
  66          fresnel_a[FRESNEL_F,i]:
  67             funmake(lambda,[[x],block([u:x-x0],ev(fresnel_Ftaylor))]),
  68          fresnel_a[FRESNEL_G,i]:
  69             funmake(lambda,[[x],block([u:x-x0],ev(fresnel_Gtaylor))])),
  70       i:i+1
  71    )
  72 )$
  73 fresnel_init():=block(
  74    fresnel_InitLinear()
  75 )$
  76 fresnelF(z):=block([z:abs(z)],
  77    if z < FRESNEL_ASYMPTOTIC then block([d:round(FRESNEL_DIVISION*z)],
  78       fresnel_a[FRESNEL_F,d](z)
  79    )
  80    else block(
  81       [ORDER:65,df:1,pz:dfloat(%pi*z),pz2:dfloat(%pi*z^2),pz2m,pzpz2,
  82           m1:1,term,n:65],
  83       fresnel_a[FRESNEL_A,0]:1/pz,
  84       pz2m:1,
  85       pzpz2:pz*pz2,
  86       for m:1 thru ORDER do (
  87          m1:m1*(-1),
  88          df:df*(4*m-3)*(4*m-1),
  89          pz2m:pz2m*(pz2^2),
  90          term:m1*dfloat((df)/(pz*pz2m)),
  91          if abs(term) > abs(fresnel_a[FRESNEL_A,m-1]) then (n:(m-1),return(n)),
  92          if abs(term/fresnel_a[FRESNEL_A,0]) < DFLOAT_EPSILON then
  93              (fresnel_a[FRESNEL_A,m]:0.0d0,n:m,return(n))
  94          else fresnel_a[FRESNEL_A,m]:term
  95       ),
  96       sum(fresnel_a[FRESNEL_A,n-i],i,0,n)
  97    )
  98 )$
  99 fresnelG(z):=block([z:abs(z)],
 100    if z < FRESNEL_ASYMPTOTIC then block([d:round(FRESNEL_DIVISION*z)],
 101       fresnel_a[FRESNEL_G,d](z)
 102    ) else block(
 103       [ORDER:65,df:1,pz:dfloat(%pi*z),pz2:dfloat(%pi*z^2),pz2m,pzpz2,
 104           m1:1,term,n:65],
 105       fresnel_a[FRESNEL_A,0]:1/(pz*pz2),
 106       pz2m:1,
 107       pzpz2:pz*pz2,
 108       for m:1 thru ORDER do (
 109          m1:m1*(-1),
 110          df:df*(4*m-1)*(4*m+1),
 111          pz2m:pz2m*(pz2^2),
 112          term:m1*dfloat((df)/(pzpz2*pz2m)),
 113          if abs(term) > abs(fresnel_a[FRESNEL_A,m-1]) then (n:(m-1),return(n)),
 114          if abs(term/fresnel_a[FRESNEL_A,0]) < DFLOAT_EPSILON then
 115              (fresnel_a[FRESNEL_A,m]:0.0d0,n:m,return(n))
 116          else fresnel_a[FRESNEL_A,m]:term
 117       ),
 118       sum(fresnel_a[FRESNEL_A,n-i],i,0,n)
 119    )
 120 )$
 121 fresnel_init()$