#
#	This program calculates the matrix elements for the Stark effect
#
#	by T. C. Scott					(Fall 1988)

#	Since only the first few coefficients are of interest, this
#       program is designed to calculate only the first few coefficients 
#	and is not intended for higher-order expansions. 

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 `Int`:
#
Delta:=proc(i,j) if i=j then 1; else 0; fi; end;
Select:=proc(n,l,choice,k,label)
	if (choice<>1) and (l<>0) then 
		k:= -l; label:=m.(abs(k)); 
	else 
		k:=l+1; label:=k;
	fi;
end;
prep:=proc(subord,ord,a)
local i,sm;
	sm:=0:
	for i from 0 to subord do
		sm := sm+ alp^(2*i)*R.a.i:
	od:
	sm :=subs(t=lam*r,A*sm):
	sm :=taylor(sm,alp,2*ord+1):
sm:
end:
prepare:=proc(ord,R1,R2)
	R1:= prep(ord,ord,1):
	R2:= prep(ord-1,ord,2):
end:
#
#   Define Zero-order Eigenfunction (upper and lower components)
#
#
	ord:=1;
	if ord > 1 then ord:=1 fi;  # limit "ord" at one
	a0:=1;
	Z:=1;
#
	n1:=1; l1:=0; mm1:=0; choice1:=1;
	n2:=2; l2:=1; mm2:=0; choice2:=1;
#
	Select(n1,l1,choice1,'k1','label1');
	label1;
	read Dirac.n1.l1.label1.`.m`;
	prepare(ord,'R1i','R2i'):
#	
	Select(n2,l2,choice2,'k2','label2');
	label2;
	read Dirac.n2.l2.label2.`.m`;
	prepare(ord,'R1f','R2f'):
#
	arg:=taylor(R1i*R1f+R2i*R2f,alp,2*ord+1):
	arg:=convert(arg,polynom):
	arg:=subs(r=t,(r^3)*arg):
	arg:=expand(simplify(normal(arg))):
	stark:=expand(map(IIntd,arg,t)):
	Almk:=Delta(mm2,mm1)*(signum(k1)*((sqrt((k1-mm1-1)*(k1+mm1)
	)*Delta(k2,k1-1)/(2*k1-1))+(sqrt((k1-mm1)*(k1+mm1+1)
	)*Delta(k2,k1+1)/(2*k1+1)))+((2*mm1+1)/((2*k1-1)*(2*k1+1))
	)*Delta(k2,-k1));
	stark0:=simplify(Almk*coeff(stark,alp,0));
	stark1:=simplify(Almk*coeff(stark,alp,2));
	evalf(stark0);
	evalf(stark1);
	save stark0,stark1,`stark.m`;
quit;
