`help/text/SHM`:=TEXT(
`HELP FOR: SHM, a procedure which generates an analytical form`,
`          for the first n harmonic oscillator levels and their`,
`          energies`,
` `,
`CALLING SEQUENCE:`,
`   SHM();`,
` `,
`SYNOPSIS:`,
`-The procedure is called with empty brackets as no parameters`,
` are passed.`,
`-The number of eigenstates is then requested. Enter the number of`,
` states required followed by a semi-colon and RETURN.`,
`-The specified states are stored in an array called SHMwf, the`,
` ground state being SHMwf(1). The corresponding energies are`,
` stored in the array SHMen.`,
` `,
`EXAMPLES:`,
`The first 5 harmonic oscillator levels are obtained by `,
`typing`,
`   SHM();`,
`then, in response to "no of levels/states to be calculated",`,
`   5;`
):

#--------------------------------------------------------------------------




SHM:=proc()
     local n,iop,nf,a,pic,pikky;
     global SHMen, SHMwf;
#calculate the nth simple harmonic oscillator state
#using Hermite polynomials from `orthopoly'
#iop is a plot option;1 for eigenvalues, 2 for eigenvectors

n:=readstat(`no of levels/states to be calculated`);
#iop:=readstat(`plot option - 0=no plot,1=energy levels,2=eigenstates`);

SHMen:=array(1..n);
SHMwf:=array(1..n);
a:=1;#a^2=MWc/h from Eisberg/Rae
for i from 1 to n do
  nf:=sqrt(1/(sqrt(Pi)*2^(i-1)*(i-1)!));
  SHMwf[i]:=nf*exp(-(a*x)^2/2)*H((i-1),a*x):
#  SHMen[i]:=(i-1+1/2);#factor 1=hWc i.e. SHO frequency
od;
print(n);print(`eigenstates have been stored in SHMwf[/n/]`);
end:
#if iop<>0 then
#pic:={};

#if iop=1 then
#  pic:={plot(y^2,y=-4..4,0..SHMen[n]+2,title=`SHO energy levels`)};
#  for i from 1 to n do
#    pikky:=plot(SHMen[i]);
#    pic:=pic union {pikky};
#  od;
#fi;

#if iop=2 then
#  for i from 1 to n do
#    pikky:=plot(SHMwf[i],y,title=`SHO eigenstates`);
#    pic:=pic union {pikky};
#  od;
#fi;

#display(pic);

#fi;

#end:

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


`help/text/Hwave`:=TEXT(
`HELP FOR: Hwave, a procedure which generates an analytical form`,
`          for a specified hydrogen atom type wavefunction`,
` `,
`CALLING SEQUENCE:`,
`   Hwave(args);`,
` `,
`SYNOPSIS:`,
`-the first 3 parameters specify the wavefunction, they are the`,
` quantum numbers n, l, m`,
`-two more parameters are optional. If unspecified the wavefunctions`,
` contain the nuclear charge number Z and the radius a. They must`,
` be specified if the wavefunctions are to be integrated`,
` `,
`EXAMPLES:`,
`The (3,2,1) wavefunction of a single electron atom is obtained by `,
`typing`,
`   Hwave(3,2,1);`,
`The corresponding hydrogen atom (Z=1) wavefunction, setting a=1,`,
`is obtained by typing`,
`   Hwave(3,2,1,1,1);`
):

#--------------------------------------------------------------------------



Hwave:=proc()
       local nf1,nf2,ctt,x,Z,a;

n:=args[1];l:=args[2];m:=args[3];
if nargs=5 then Z:=args[4];a:=args[5]; fi;

#solve the time independent Schrodinger equation for the hydrogen atom
#take the known solutions in terms of Laguerre and Legendre polynomials
#[with(orthopoly)] have to convert these to associated polynomials
#n,l,m are the orbital quantum numbers.

#Z:=1;a:=1;

#associated Laguerre polynomial
lag:=diff(L(n+l,x),x$(2*l+1));#print(");
lag:=subs(x=2*Z*r/(n*a),lag)*exp(-Z*r/(n*a))*(2*Z*r/(n*a))**(l);

#associated Legendre polynomial
leg:=P(l,ctt);#print(leg);
if abs(m)>0 then
  leg:=(1-ctt**2)**(abs(m)/2)*diff(leg,ctt$abs(m));#print(");
fi;
leg:=subs(ctt=cos(t),leg);

nf1:=-sqrt((2*Z/(n*a))^3*(n-l-1)!/((2*n)*((n+l)!)^3));
#nf1:=1;
nf2:=(-1)^m*sqrt((2*l+1)*(l-abs(m))!/(4*Pi*(l+abs(m))!));
#nf2:=2;
wf:=nf1*nf2*lag*leg*exp(I*m*f);


end:

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Th:=proc(l,m)

leg:=P(l,ctt);
leg:=(1-ctt**2)**(abs(m)/2)*diff(leg,ctt$abs(m));
leg:=subs(ctt=cos(t),leg);

end:

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

R:=proc(n,l)

lag:=diff(L(n+l,r),r$(2*l+1));

end:

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

with(plots):with(orthopoly):