123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408 |
- module sfellipi; % Procedures and Rules for Elliptic Integrals.
- % Author: Lisa Temme, ZIB, October 1994
- algebraic <<
- %######################################################################
- %DESCENDING LANDEN TRANSFORMATION
- procedure landentrans(phi,alpha);
- begin scalar alpha_n!+1, alpha_n, phi_n!+1, phi_n, aNtoa0, pNtop0,
- a0toaN, p0topN;
- alpha_n := alpha;
- phi_n := phi;
- aNtoa0 := {alpha_n};
- pNtop0 := {phi_n};
- while alpha_n > 10^(-(Symbolic !:prec!:)) do
- <<
- alpha_n!+1:= asin(2/(1+cos(alpha_n)) -1);
- phi_n!+1 := phi_n + (atan(cos(alpha_n)*tan(phi_n)))
- + floor((floor(phi_n/(pi/2))+1)/2)*pi;
- aNtoa0 := alpha_n!+1.aNtoa0;
- pNtop0 := phi_n!+1.pNtop0;
- alpha_n := alpha_n!+1;
- phi_n := phi_n!+1
- >>;
- a0toaN := reverse(aNtoa0);
- p0topN := reverse(pNtop0);
- return list(p0topN, a0toaN)
- end;
- %######################################################################
- %VALUE OF EllipticF(phi,m)
- procedure F_function(phi,m);
-
- begin scalar alpha, bothlists, a0toaN, a1toaN, p0topN, phi_n, y,
- elptF;
- alpha := asin(sqrt(m));
- bothlists := landentrans(phi,alpha);
- a0toaN := PART(bothlists,2);
- a1toaN := REST(a0toaN);
- p0topN := PART(bothlists,1);
- phi_n := PART(reverse(p0topN),1);
- if phi = (pi/2)
- then
- elptF := K_function(m)
- else
- elptF :=
- phi_n *for each y in a1toaN PRODUCT(1/2)*(1+sin(y));
- return elptF
- end;
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %EllipticF definition
- %====================
- operator EllipticF;
- EllipticFrules :=
- {
- EllipticF(~phi,0) => phi,
- EllipticF(i*~phi,0) => i*phi,
- EllipticF(~phi,1) => ln(sec(phi)+tan(phi)),
- EllipticF(i*~phi,1) => i*atan(sinh(phi)),
- EllipticF(~phi,~m) => Num_Elliptic(F_function,phi,m)
- when lisp !*rounded and numberp phi
- and numberp m
- };
- let EllipticFrules;
- %######################################################################
- %VALUE OF K(m)
- procedure K_function(m);
- begin scalar AGM, aN;
- AGM := AGM_function(1,sqrt(1-m),sqrt(m));
- aN := PART(AGM,2);
- return (pi / (2*aN));
- end;
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %EllipticK definition
- %====================
- EllipticKrules :=
- {
- EllipticK(~m) => K_function(m) when lisp !*rounded
- and numberp m,
- EllipticK!'(~m) => K_function(1-m) when lisp !*rounded
- and numberp m
- };
- let EllipticKrules;
- %######################################################################
- %VALUE OF EllipticE(phi,m)
- procedure E_function(phi,m);
- begin scalar F, N, alpha, bothlists, a0toaN, p0topN, a1toaN, p1topN,
- sinalist, sinplist, b, s, blist, c, allz, w, z, allx,
- h, x, elptE;
- F := F_function(phi,m);
- alpha := asin(sqrt(m));
- bothlists := landentrans(phi,alpha);
- a0toaN := PART(bothlists, 2);
- p0topN := PART(bothlists, 1);
- a1toaN := REST(a0toaN);
- p1topN := REST(p0topN);
- N := LENGTH(a1toaN);
- sinalist := sin(a1toaN);
- sinplist := sin(p1topN);
- b := PART(sinalist,1);
- s := b;
- blist := for each c in rest sinalist collect << b := b*c >>;
- blist := s.blist;
- allz := 0;
- for w := 1:N do
- <<
- z := (1/(2^w))*PART(blist,w);
- allz := allz + z
- >>;
- allx := 0;
- for h := 1:N do
- <<
- x := (1/(2^h))*((PART(blist,h))^(1/2))
- * PART(sinplist,h);
- allx := allx + x
- >>;
- elptE := F * (1 - (1/2)*((sin(PART(a0toaN,1)))^2)*(1 + allz))
- + sin(PART(a0toaN,1))*allx ;
- return elptE;
- end;
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %EllipticE(phi,m) definition
- %====================
- operator EllipticE;
- JacobiErules :=
- {
- EllipticE(0,~m) => 0,
- EllipticE(~phi,0) => phi,
- EllipticE(i*~phi,0) => i*phi,
- EllipticE(~phi,1) => sin(phi),
- EllipticE(i*~phi,1) => i*sinh phi,
- EllipticE(-~phi,~m) => -EllipticE(phi,m),
- EllipticE(~phi,-~m) => EllpiticE(phi,m),
- df(EllipticE(~phi,~m),~phi) => Jacobidn(phi,m)^2,
- df(EllipticE(~phi,~m),~m) =>
- m * (Jacobisn(phi,m) * Jacobicn(phi,m) * Jacobidn(phi,m)
- - EllipticE(phi,m) * Jacobicn(phi,m)^2) / (1-m^2)
- - m * phi * Jacobisn(phi,m)^2,
- EllipticE(~phi,~m) => Num_Elliptic(E_function,phi,m)
- when lisp !*rounded and numberp phi
- and numberp m,
- EllipticE(~m) => Num_Elliptic(E_function,pi/2,m)
- when lisp !*rounded and numberp m
- };
- let JacobiErules;
- %######################################################################
- %CALCULATING THE FOUR THETA FUNCTIONS
- %Theta 1 (often written H(u) - and has period 4K)
- %Theta 2 (often written H1(u) -and has period 4K)
- %Theta 3 (often written Theta1(u) - and has period 2K)
- %Theta 4 (often written Theta(u) - and has period 2K)
- procedure num_theta(a,u,m);
- begin scalar n, new, all, z, q, total;
- n := if a>2 then 1 else 0;
- new := 100; % To initiate loop
- all := 0;
- z := (pi*u)/(2*EllipticK(m));
- q := EXP(-pi*EllipticK(1-m)/EllipticK(m));
- while new > 10^(-(Symbolic !:prec!:)) do
- << new := if a =1 then
- ((-1)^n)*(q^(n*(n+1)))*sin((2*n+1)*z)
- else if a=2 then (q^(n*(n+1)))*cos((2*n+1)*z)
- else if a=3 then (q^(n*n))*cos(2*n*z)
- else if a=4 then ((-1)^n)*(q^(n*n))*cos(2*n*z);
- all := new + all;
- n := n+1
- >>;
- return if a > 2 then (1 + 2*all)
- else (2*(q^(1/4))*all);
- end;
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %Theta Functions
- operator EllipticTheta;
- EllipticTHETArules :=
- {
- %Theta1rules
- %-----------
- EllipticTheta(1,~u,~m) =>
- Num_Elliptic(num_theta,1,u,m) when lisp !*rounded
- and numberp u and numberp m,
- EllipticTheta(1,-~u,~m) => -EllipticTheta(1,u,m),
- EllipticTheta(1,~u+EllipticK(~m),~m) => EllipticTheta(2,u,m),
- EllipticTheta(1,~u+(2*EllipticK(~m)),~m) =>
- -EllipticTheta(1,u,m),
- EllipticTheta(1,~u+i*EllipticK!'(~m),~m) =>
- i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(4,u,m),
- EllipticTheta(1,~u+2*i*EllipticK!'(~m),~m) =>
- -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(1,u,m),
- EllipticTheta(1,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(3,u,m),
- EllipticTheta(1,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(1,u,m),
- %Theta2rules
- %-----------
- EllipticTheta(2,~u,~m) =>
- Num_Elliptic(num_theta,2,u,m) when lisp !*rounded
- and numberp u and numberp m,
- EllipticTheta(2,-~u,~m) => EllipticTheta(2,u,m),
- EllipticTheta(2,~u+EllipticK(~m),~m) => -EllipticTheta(1,u,m),
- EllipticTheta(2,~u+(2*EllipticK(~m)),~m) =>
- -EllipticTheta(2,u,m),
- EllipticTheta(2,~u+i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(3,u,m),
- EllipticTheta(2,~u+2*i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(2,u,m),
- EllipticTheta(2,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
- -i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(4,u,m),
- EllipticTheta(2,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
- -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(2,u,m),
- %Theta3rules
- %-----------
- EllipticTheta(3,~u,~m) =>
- Num_Elliptic(num_theta,3,u,m) when lisp !*rounded
- and numberp u and numberp m,
- EllipticTheta(3,-~u,~m) => EllipticTheta(3,u,m),
- EllipticTheta(3,~u+EllipticK(~m),~m) => EllipticTheta(4,u,m),
- EllipticTheta(3,~u+(2*EllipticK(~m)),~m) =>
- EllipticTheta(3,u,m),
- EllipticTheta(3,~u+i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(2,u,m),
- EllipticTheta(3,~u+2*i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(3,u,m),
- EllipticTheta(3,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
- i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(1,u,m),
- EllipticTheta(3,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(3,u,m),
- %Theta4rules
- %-----------
- EllipticTheta(4,~u,~m) =>
- Num_Elliptic(num_theta,4,u,m) when lisp !*rounded
- and numberp u and numberp m,
- EllipticTheta(4,-~u,~m) => EllipticTheta(4,u,m),
- EllipticTheta(4,~u+EllipticK(~m),~m) => EllipticTheta(3,u,m),
- EllipticTheta(4,~u+(2*EllipticK(~m)),~m)=>EllipticTheta(4,u,m),
- EllipticTheta(4,~u+i*EllipticK!'(~m),~m) =>
- i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(1,u,m),
- EllipticTheta(4,~u+2*i*EllipticK!'(~m),~m) =>
- -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(4,u,m),
- EllipticTheta(4,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
- (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
- *EllipticTheta(2,u,m),
- EllipticTheta(4,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
- -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
- *EllipticTheta(4,u,m),
- %Error
- %-----
- EllipticTheta(~a,~u,~m) =>
- printerr ("In EllipticTheta(a,u,m); a = 1,2,3 or 4.")
- when numberp a
- and not(fixp a and a<5 and a>0)
- };
- let EllipticTHETArules;
- %######################################################################
- %CALCULATING ZETA
- procedure ZETA_function(u,m);
- begin scalar phi_list, clist, L, j, z, cn, phi_n;
- phi_list := PHI_function(1,sqrt(1-m),sqrt(m),u);
- clist := PART(AGM_function(1,sqrt(1-m),sqrt(m)),5);
- L := LENGTH(phi_list);
- j := 1;
- z := 0;
- while j < L do
- <<
- cn := PART(clist,L-j);
- phi_n := PART(phi_list,1+j);
- z := cn*sin(phi_n) + z;
- j := j+1
- >>;
- return z
- end;
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %JacobiZETA definition
- %=====================
- operator JacobiZeta;
- JacobiZETArules :=
- {
- JacobiZeta(~u,0) => 0,
- JacobiZeta(~u,1) => tanh(u),
- JacobiZeta(-~u,~m) => -JacobiZeta(u,m),
- JacobiZeta(~u+~v,~m) => JacobiZeta(u,m) + JacobiZeta(v,m) -
- (m*Jacobisn(u,m)*Jacobisn(v,m)
- *Jacobisn(u+v,m)),
- JacobiZeta(~u+2*EllipticK(~m),m) => JacobiZeta(u,m),
- JacobiZeta(EllipticK(~m) - ~u,m) =>
- -JacobiZeta(EllipticK(m)+u,m),
- % JacobiZeta(~u,~m) => JacobiZeta(u - EllipticK(m),m) -
- % m * Jacobisn(u - EllipticK(m),m)
- % * Jacobicd(u - EllipticK(m),m),
- JacobiZeta(~u,~m) => Num_Elliptic(ZETA_function,u,m)
- when lisp !*rounded and numberp u
- and numberp m
- };
- let JacobiZETArules;
- %######################################################################
- >>;
- endmodule;
- end;
|