#	Program DIRAC1
#
#	written by T.C. Scott, G.J. Fee and M.B. Monagan	(1988)
#
#	This program solves both the wavefunction and eigenergies of the 
#	Dirac equation according to Moore's decoupling procedure.
#
#	This approach obtains a series solution by perturbation theory 
#	and solves the wavefunction to order O(j) and the energies
#	to order O(2*j+1).


fix := proc(f) local a,g,i,terms,tt,tc,coffs,s;
	g := subs(exp(t/2) = exp(t)^(1/2), f);
	g := collect(g,t,distributed,normal);
g;
end:

read `Integration`;

Inner:=proc(Ru1,Rd1,Ru2,Rd2,t)
local arg; option remember,system;
	arg:=fix((t^2)*(Ru1*Ru2 + Rd1*Rd2)/lam3):
	fix(A2*IIntd(arg,t)):
end:
M:=proc(i,j)
option remember;
	M(i,j):=Inner(R1.i,R2.i,R1.j,R2.j,t):
	M(j,i):=M(i,j):
end:

#   Generate Zero-order Solution

lag:=proc(n,l)
local s;
	Lnl:=0;
	Anl:=sqrt( lam3*( (n-l-1)!/(2*n*((n+l)!^3)) ) );
	for s from 0 to (n-l-1) do
		la.s:=((n+l)!^2)/((n-l-1-s)!*(2*l+1+s)!*(s!));
		Lnl:=Lnl + ((-1)^s)*( la.s )*t^s;
	od;
	fix(Anl*1/exp(t)^(1/2)*(t^l)*Lnl);
end;
read `Ghost`;
Inhomo:=proc(f)
local f1,f2,u,v;
#
#   solve Inhomogeneous equation by method of Variation of parameters
#
	f1:=fix(t^2*f*y1):
	f2:=fix(t^2*f*y2):
	u:=-fix(IInt(f2,t)):
	v:= fix(IInt(f1,t)):
	fix(u*y1 + v*y2):
end;

#   Main Program
	veri:=1;
	if veri=1 then read Dirac1veri fi;

#   Inititialise Zero-order solution

	n:=1; l:=0; choice:=1;
	finalloop:=4;
#
	if (choice<>1) and (l<>0) then 
		k:= -l; label:=m.(abs(k)); 
	else 
		k:=l+1; label:=k;
	fi;
#
#   Define constants
#
#    ... Without any loss of generality we can set:
#
	Z:=1;			# the "Z" just becomes absorbed into "alp"
	a0:=1; 			# it's value will not affect the final energies
#
	A:=1/sqrt(1+(Z*alp/(2*n))^2);
	A2:=fix(A^2);
	lam:=fix(2*Z*A2/(n*a0));	lam3:=fix(lam^3);
	Ze2:=2*Z*a0*RH;			mc2:=2*RH/(alp^2);
	ch:=2*a0*RH/alp;		c02:=fix(lam*ch/(2*mc2*A2));
#
#   Define operators
#
	Pm1:=<diff(R,t) - ((k-1)*R/t)|R>:
	Pp1:=<diff(R,t) + ((k+1)*R/t)|R>:
#
#   Define Zero-order Eigenfunction (upper and lower components)
#
	R10:=lag(n,l); 
	R20:=fix(c02*Pm1(R10)):
	y1:=fix(R10/Anl):
#
#   Define Potential and zero-order energies
#
	V:=-fix(Ze2*lam/t);	
	eps00:=-RH*((Z/n)^2);	eps0:=fix(eps00*A2);
#
#   Get "Ghost" state 
#
	y2:=fix(Ghost2(y1));
	if veri=1 then Veri2(y1,1); Veri2(y2,2) fi:

#
#   Calculation of first-order Energies
#
	arg:=fix(t^2*R20*V*R20/lam3):
	E1:=fix((A2/alp^2)*IIntd(arg,t));

# reset all integration counters to zero
counter:=0: counter2:=0: counter3:=0: counter4:=0:
#
#   Main Loop
#
for loop from 1 to finalloop do

#  Solve loop'th inhomogeneous equation

	f:=0;
	for i from 1 to loop do
		j:=loop-i;
		temp:=fix(c02*Pp1(R2.j)):
		temp:=fix(E.i*(R1.j - temp));
		f:=f+temp;
	od:
	arg:=fix(V*R2.(loop-1)):
	arg:=fix(((1/alp)^2)*c02*Pp1(arg));
	f:=fix((f+arg)/(4*eps0)):
#	
	yp.i:=Inhomo(f):
	R1.loop:=fix(yp.i + C.loop*R10):
#
#   Define Lower component
#
	R2.loop:=fix(c02*Pm1(R1.loop)):
	catch:=fix(( (V/(alp^2)) - E1)*R2.(loop-1)):
	for i from 2 to loop do
		j:=loop-i:
		catch:= catch - fix(E.i*(R2.j)):
	od:
	R2.loop:=R2.loop + fix((1/(2*mc2*A2))*catch):
	R2.loop:=fix(R2.loop);

	if (loop mod 2)=0 then
		eqC:=M(loop/2,loop/2):
		lasti:=(loop/2) -1:
	else
		eqC:=0:
		lasti:=(loop-1)/2;
	fi;
	for i from 0 to lasti do 
		eqC:=eqC + 2*M(i,loop-i):
	od:
	C.loop:=solve(eqC,C.loop):
	C.loop:=fix(C.loop):
	R1.loop:=fix(R1.loop):
	R2.loop:=fix(R2.loop):
	yp.loop:=R1.loop:
	if veri=1 then Veri3(loop) fi;

#   Calculate Energies
#
#   E(2*n)
#
# 	First part

	arg:=fix(t^2*R2.(loop-1)*V*R2.loop/lam3):
	E.(2*loop).1:=fix((A2/alp^2)*IIntd(arg,t)):
#
#       Second part
#
	E.(2*loop).2 :=0:
	for i from 1 to (loop-1) do
		sm:=0:
		for j from 1 to i do
			sm:=sm+M(loop-j,loop+j-i):
		od:
		E.(2*loop).2:=E.(2*loop).2 + fix(E.i*sm):
	od:
	E.(2*loop):=E.(2*loop).1 - E.(2*loop).2 + (E.loop)*M(0,loop):
	E.(2*loop):=factor(fix(E.(2*loop))):
	print(E.(2*loop));
#
#   E(2*n+1)
# 
# 	First part
#
	arg:=fix(t^2*R2.loop*V*R2.loop/lam3):
	E.(2*loop+1).1:=fix((A2/alp^2)*IIntd(arg,t)):
#
#       Second part
#
	E.(2*loop+1).2:=fix(E1*M(loop,loop)):
	for i from 2 to loop do
		if (i mod 2)=0 then
			lastj:=i/2;
			sm:=0;
		else
			lastj:=(i-1)/2;
			sm:=M(loop-lastj,loop-lastj):
		fi;
		for j from 1 to lastj do
			sm:=sm+2*M(loop-j+1,loop+j-i):
		od:
		E.(2*loop+1).2:=E.(2*loop+1).2 + fix(E.i*sm):
	od:
	E.(2*loop+1):=E.(2*loop+1).1 - E.(2*loop+1).2:
	E.(2*loop+1):=factor(fix(E.(2*loop+1))):
	print(E.(2*loop+1));
	if veri=1 then Veri4(loop) fi;
	# reset all integration counters to zero
	counter:=0: counter2:=0: counter3:=0: counter4:=0:
od:
#
#  Semi-final result
#
	mc2:='mc2':
	Edirac:=mc2/sqrt( (1 + (Z*alp/(n-abs(k) + sqrt(k^2-(Z*alp)^2)))^2) );
	Ecalc:=mc2 + eps0:
	for i from 1 to (2*finalloop+1) do
		Ecalc:=Ecalc + (alp^(2*i))*E.i:
	od:
	RH:=(1/2)*mc2*alp^2:
#
#  Final Verification
#
#  Exact energy
#
	far:=6+2*(2*finalloop+1);
	fix(taylor(Edirac,alp,far));
#
#  Calculated energy
#
	fix(taylor(Ecalc,alp,far));
#
save n,l,k,A,A2,lam,eps0,E.(1..(2*loop+1)),Ecalc,Edirac,RH,
		R1.(0..loop),R2.(0..loop), Dirac.n.l.label.`.m`;
quit;