share/contrib/fresnel/fresnel1.mac

   1 /* This function computes the fresnel F and G functions with dfloat accuracy
   2 over -inf,inf.  It divides the real line into 3 regions.  The first (0,1) uses
   3 a division of 1/100, and a taylor series.  The second region (1,5)
   4 uses reciprocal divisions, and a taylor series also.  The third region (5,inf)
   5 uses the asymptotic series.  This routine could be improved by using
   6 lazy evaluation to increase speed in the initialization.
   7 */
   8
   9 load("fresnel_linear_values.mac")$
  10 define_variable(FRESNEL_RECIPROCAL,1.0d0,any_check,
  11    "Point to switch to reciprocal values")$
  12 define_variable(FRESNEL_ASYMPTOTIC,5.0d0,any_check,
  13    "Point to switch to asymptotic series")$
  14 define_variable(FRESNEL_DIVISION,100,any_check,
  15    "Number of divisions")$
  16 define_variable(FRESNEL_X,0,any_check, "Array constant")$
  17 define_variable(FRESNEL_F,1,any_check, "Array constant")$
  18 define_variable(FRESNEL_G,2,any_check, "Array constant")$
  19 array(fresnel_a,2,2*FRESNEL_DIVISION+1)$
  20 define_variable(fresnel_Ftaylor, 8.33333333333334d-3 *
  21    (u *
  22      (u *
  23        (u *
  24          (u *
  25            (u * x0 *
  26              (x0 *
  27               (x0 * (9.74090910340024d+2 * f0 - 3.0601968478528d+2 * g0 * x0^2)
  28                - 2.79056490122698d+2)
  29               + 4.65094150204496d+2 * g0)
  30              + x0^2 *
  31                  (x0 * (4.87045455170012d+2 * f0 * x0 - 1.55031383401498d+2)
  32                   + 9.30188300408994d+2 * g0)
  33             - 1.4804406601634d+2 * f0)
  34            + x0 * (6.20125533605996d+2 * g0 * x0^2 - 5.9217626406536d+2 * f0)
  35           + 1.25663706143592d+2)
  36          + x0 * (1.88495559215388d+2 - 5.9217626406536d+2 * f0 * x0)
  37         - 1.88495559215388d+2 * g0)
  38       - 3.76991118430774d+2 * g0 * x0)
  39     + 1.2d+2 * f0),any_check,"Taylor series for F about x0")$
  40 define_variable(fresnel_Gtaylor, 8.33333333333334d-3 *
  41    (u *
  42     (u *
  43      (u *
  44       (u *
  45        (x0 *
  46         (x0 *
  47          (4.87045455170012d+2 * g0 * x0^2 - 9.30188300408994d+2 * f0)
  48          + 2.46740110027234d+2)
  49         + u *
  50           (x0 *
  51            (x0^2 *
  52             (x0 *
  53              (3.0601968478528d+2 * f0 * x0 - 9.74090910340024d+1)
  54              + 9.74090910340024d+2 * g0)
  55             - 4.65094150204496d+2 * f0)
  56            + 7.89568352087148d+1)
  57         - 1.4804406601634d+2 * g0)
  58        + x0 * (x0 * (1.97392088021786d+2 - 6.20125533605996d+2 * f0 * x0)
  59        - 5.9217626406536d+2 * g0))
  60       - 5.9217626406536d+2 * g0 * x0^2 + 1.88495559215388d+2 * f0)
  61      + 3.76991118430774d+2 * f0 * x0 - 1.2d+2)
  62     + 1.2d+2 * g0),any_check,"Taylor series for G about x0")$
  63
  64 fresnel_InitLinear():=block([i:0],
  65    for l in fresnel_linear_values do (
  66       fresnel_a[FRESNEL_X,i]:first(l),
  67       block([x0:first(l),f0:second(l),g0:third(l)],
  68          fresnel_a[FRESNEL_F,i]:
  69             funmake(lambda,[[x],block([u:x-x0],ev(fresnel_Ftaylor))]),
  70          fresnel_a[FRESNEL_G,i]:
  71             funmake(lambda,[[x],block([u:x-x0],ev(fresnel_Gtaylor))])),
  72       i:i+1
  73    )
  74 )$
  75 fresnel_init():=block(
  76    fresnel_InitLinear()
  77 )$
  78 fresnelF(x):=block()$
  79 fresnelG(x):=block()$
  80 fresnel_init()$