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;



