Artifact cbddf7c46e00bc8fe088394eaa1ea6832aefb5fa1ca7468f3ed66796a6395380:
- Executable file
r37/packages/specfn/sfellipi.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 10906) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/sfellipi.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 10906) [annotate] [blame] [check-ins using]
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;