share/contrib/fresnel/fresnelstudy.mac

   1 matchdeclare(a,NON_ZERO_AND_FREEOF(t));
   2 matchdeclare([b,c],freeof(t));
   3 defrule(c1,'integrate(cos(t *b +a* t^2 +c),t),
   4    sqrt(%pi/(2*a))*(cos((b^2-a*c)/a)*cfresnel(sqrt(2/(a*%pi))*(a*t+b))
   5                    +sin((b^2-a*c)/a)*sfresnel(sqrt(2/(a*%pi))*(a*t+b))));
   6 defrule(s1,'integrate(sin(t *b +a* t^2 +c),t),
   7    sqrt(%pi/(2*a))*(cos((b^2-a*c)/a)*sfresnel(sqrt(2/(a*%pi))*(a*t+b))
   8                    -sin((b^2-a*c)/a)*cfresnel(sqrt(2/(a*%pi))*(a*t+b))));
   9 matchdeclare(q,freeof(%i));
  10 defrule(cfi,cfresnel(%i*q),%i*cfresnel(q));
  11 defrule(sfi,sfresnel(%i*q),-%i*sfresnel(q));
  12 trigreduce(demoivre(subst([ia=%i*w,ib= - 15840000 * %pi, ic=16000000 * %pi,
  13  id=16000000 * %pi ], (exp(-ia*t)*(cos(ib*t+ic*t^2)+cos(id*t))))));
  14 map(lambda([x],'integrate(x,t)),%);
  15
  16 fresnelF(u):=block([x:abs(u)],signum(u)*(1.0+0.926*x)/(2.0+1.792*x+3.104*x^2));
  17 fresnelG(u):=block([x:abs(u)],signum(u)/(2.0+4.142*x+3.492*x^2+6.670*X^3));
  18 cfresnel(x):=block([fg:dfloat(((1+%i)/2)*werf(abs(x)*(1+%i)*sqrt(%pi)/2)),f,g],
  19    f:realpart(fg),g:imagpart(fg),
  20    signum(x)*(0.5+f*sin((%pi/2)*x^2)-g*cos((%pi/2)*x^2))
  21 )$
  22 sfresnel(x):=block([fg:dfloat(((1+%i)/2)*werf(abs(x)*(1+%i)*sqrt(%pi)/2)),f,g],
  23    f:realpart(fg),g:imagpart(fg),
  24    signum(x)*(0.5d0-f*cos((%pi/2)*x^2)-g*sin((%pi/2)*x^2))
  25 )$
  26 addem(l):=block([p:0,n:0],
  27    l:sort(l,lambda([x,y],if abs(x) < abs(y) then true else false)),
  28    for i in l do if i > 0 then p:p+i else n:n+i,
  29    p+n);
  30
  31 matchdeclare(mv,any);
  32 defrule(fg1,cfresnel(mv),(1/2+f(mv)*sin((%pi/2)*mv^2)-g(mv)*cos((%pi/2)*mv^2)));
  33 defrule(fg2,sfresnel(mv),(1/2-f(mv)*cos((%pi/2)*mv^2)-g(mv)*sin((%pi/2)*mv^2)));
  34 g:realpart(dfloat(expand((1+%i)/2*werf(dfloat((sqrt(%pi)/2)*(1+%i)*40*sqrt(2))))));
  35 1/((%pi*x)*(%pi*(x)^2));
  36 defrule(ga,fresnelG(mv),1/((%pi*mv)*(%pi*(mv)^2)));
  37 defrule(fa,fresnelF(mv),1/(%pi*mv)+3/(%pi*(mv^2))^2);
  38 maplist(dfloat,expand(taylor(d17,w,16000000 * %pi,0)));
  39 funcargs(d78,[fresnelF,fresnelG]);
  40 atvalue(f(x),[x=0],1/2);
  41 atvalue(g(x),[x=0],1/2);
  42 atvalue(fresnelf(x),[x=0],1/2);
  43 atvalue(fresnelg(x),[x=0],1/2);
  44 atvalue(fresnelf(x),x=((1596000001 * sqrt(2))/(5000 * sqrt(5))),0);
  45
  46 array(fresnel_array,65);
  47 fresnelG(z):=block(
  48    if abs(z) > 4.0d0 then (
  49       [ORDER:65,df:1,pz:dfloat(%pi*z),pz2:dfloat(%pi*z^2),pz2m,pzpz2,
  50           m1:1,term,n:65],
  51       fresnel_array[0]:1/(pz*pz2),
  52       pz2m:1,
  53       pzpz2:pz*pz2,
  54       for m:1 thru ORDER do (
  55          m1:m1*(-1),
  56          df:df*(4*m-1)*(4*m+1),
  57          pz2m:pz2m*(pz2^2),
  58          term:m1*dfloat((df)/(pzpz2*pz2m)),
  59          if abs(term) > abs(fresnel_array[m-1]) then (n:(m-1),return(n)),
  60          if abs(term/fresnel_array[0]) < DFLOAT_EPSILON then (n:m,return(n))
  61          else fresnel_array[m]:term
  62       ),
  63       sum(fresnel_array[n-i],i,0,n))
  64    else
  65       realpart(dfloat(expand((1+%i)/2*werf(dfloat((sqrt(%pi)/2)*(1+%i)*z)))))
  66 )$
  67 fresnelF(z):=block(
  68    if abs(z) > 4.0d0 then (
  69       [ORDER:65,df:1,pz:dfloat(%pi*z),pz2:dfloat(%pi*z^2),pz2m,pzpz2,
  70           m1:1,term,n:65],
  71       fresnel_array[0]:1/pz,
  72       pz2m:1,
  73       pzpz2:pz*pz2,
  74       for m:1 thru ORDER do (
  75          m1:m1*(-1),
  76          df:df*(4*m-3)*(4*m-1),
  77          pz2m:pz2m*(pz2^2),
  78          term:m1*dfloat((df)/(pz*pz2m)),
  79          if abs(term) > abs(fresnel_array[m-1]) then (n:(m-1),return(n)),
  80          if abs(term/fresnel_array[0]) < DFLOAT_EPSILON then (n:m,return(n))
  81          else fresnel_array[m]:term
  82       ),
  83       sum(fresnel_array[n-i],i,0,n))
  84    else
  85       realpart(dfloat(expand((1+%i)/2*werf(dfloat((sqrt(%pi)/2)*(1+%i)*z)))))
  86 )$
  87
  88 'diff(fresnelF,x)=-%pi*x*fresnelG;
  89 'diff(fresnelG,x)=%pi*x*fresnelF-1;
  90 taylor_ode(['diff(fresnelF,x)=-%pi*x*fresnelG,'diff(fresnelG,x)=%pi*x*fresnelF-1],
  91 [fresnelF,fresnelG],x,3,[0,[1/2],[1/2]]);
  92 reverse(maplist(lambda([u],dfloat(ev(u,x=1/%pi))),(1/2) - x + ((%pi * x^2)/4) - (((%pi)^2 * x^4)/16) + (((%pi)^2 * x^5)/15) - (((%pi)^3 * x^6)/96) + (((%pi)^4 * x^8)/768) - (((%pi)^4 * x^9)/945) + (((%pi)^5 * x^10)/7680) ));
  93
  94 'diff(fresnelg,x,2)*x = -%pi^2*fresnelg*x^3+'diff(fresnelg,x,1)+1;
  95 ev(((%pi)^(2*n)/((4*n+1)!!))*(z^(4*n+1)),n=2);
  96 map(lambda([n],dfloat(ev(((%pi)^(2*n)/((4*n+1)!!))*(z^(4*n+1)),z=1.0d0))),[0,1,2,4]);
  97 makelist(dfloat(((-1)^m)*ev((4*m+1)!!/((%pi*z)*((%pi*z^2)^(2*m+1))),z=1/0.3d0)),m,0,8);
  98
  99 fresnelA(z):=block([nn,BFLOAT_PRECISION:24],
 100    for n:0 step 1 do (
 101       if abs(dfloat(((%pi)^(2*n)/((4*n+1)!!))*(z^(4*n+1)))) < DFLOAT_EPSILON
 102          then (nn:n,return(nn))
 103    ),
 104    sum(bfloat(((-1)^n)*((%pi)^(2*n)/((4*n+1)!!))*(z^(4*n+1))),n,0,nn)
 105 )$
 106
 107 fresnelB(z):=block([nn,BFLOAT_PRECISION:24],
 108    for n:0 step 1 do (
 109       if abs(dfloat(((%pi)^(2*n+1)/((4*n+3)!!))*(z^(4*n+3)))) < DFLOAT_EPSILON
 110          then (nn:n,return(nn))
 111    ),
 112    sum(bfloat(((-1)^n)*((%pi)^(2*n+1)/((4*n+3)!!))*(z^(4*n+3))),n,0,nn)
 113 )$
 114 fresnelC(z):=block([nn,BFLOAT_PRECISION,term,az:dfloat(abs(z)),
 115    BL:ceiling(abs(log10(DFLOAT_PRECISION))),FLOAT2BF:true],
 116    BFLOAT_PRECISION:BL,
 117    for n:0 step 1 do (
 118       term:dfloat(((%pi/2)^(2*n)/(((2*n)!)*(4*n+1)))*(az^(4*n+1))),
 119       lterm:ceiling(log10(term)),
 120       if (lterm+BL) > BFLOAT_PRECISION then BFLOAT_PRECISION:lterm+BL,
 121       if term < DFLOAT_EPSILON then (nn:n,return(nn))
 122    ),
 123    BFLOAT_PRECISION:BFLOAT_PRECISION+2,
 124    z:bfloat(z),
 125    sum(bfloat((-1)^n*((%pi/2)^(2*n)/(((2*n)!)*(4*n+1)))*(z^(4*n+1))),n,0,nn)
 126 )$
 127 fresnelS(z):=block([nn,BFLOAT_PRECISION,term,az:dfloat(abs(z)),
 128    BL:ceiling(abs(log10(DFLOAT_PRECISION))),FLOAT2BF:true],
 129    BFLOAT_PRECISION:BL,
 130    for n:0 step 1 do (
 131       term:dfloat(((%pi/2)^(2*n+1)/(((2*n+1)!)*(4*n+3)))*(az^(4*n+3))),
 132       lterm:ceiling(log10(term)),
 133       if (lterm+BL) > BFLOAT_PRECISION then BFLOAT_PRECISION:lterm+BL,
 134       if term < DFLOAT_EPSILON then (nn:n,return(nn))
 135    ),
 136    BFLOAT_PRECISION:BFLOAT_PRECISION+2,
 137    z:bfloat(z),
 138    sum(bfloat((-1)^n*((%pi/2)^(2*n+1)/(((2*n+1)!)*(4*n+3)))*(z^(4*n+3))),n,0,nn)
 139 )$
 140 rncombine(expand(linsolve([fc = f * sz2 - cz2 * g + (1/2), fs =  - g * sz2 - cz2 * f + (1/2)], [f, g])));
 141 fresnelFG(z):=block([fc,fs,cz2,sz2,zp2,FLOAT2BF:true],
 142    z:bfloat(z),
 143    zp2:bfloat(%pi)*z*z/2.0b0,
 144    cz2:cos(zp2),
 145    sz2:sin(zp2),
 146    fc:fresnelC(z),
 147    fs:fresnelS(z),
 148    [z,
 149     ((2 * fc * sz2 - sz2 - 2 * cz2 * fs + cz2)/(2 * sz2^2 + 2 * cz2^2)),
 150     (( - 2 * fs * sz2 + sz2 - 2 * cz2 * fc + cz2)/(2 * sz2^2 + 2 * cz2^2))]
 151 )$
 152 ODNS_EPS:1.0d-14;
 153 ode_numsol(['diff(fresnelF,x)=-%pi*x*fresnelG,'diff(fresnelG,x)=%pi*x*fresnelF-1],[fresnelF,fresnelG],x,['at(fresnelg,x = 1) = 0.06174085260964d0, 'at(fresnelf,x = 1) = 0.27989340037682d0],1.0d0,1.005d0,1.005d0,'rk5a);
 154
 155 subst([sin(%pi*u/25)=%pi*u/25,sin(%pi/1250)=%pi/1250,cos(%pi*u/25)=1,cos(%pi/1250)=1],d93);