/* 量子力学の演算子法 水素原子 続
量子力学の演算子法(調和振動子 ∼ 水素原子)
緒方秀教 電気通信大学

と云うPDFに書いてあった事から波動関数の計算をMaxima Script で書いてみました。

*/
depends([Rw, f, fw, phi,psi,dum ] , r) ;
declare([l , n ] , integer) ;
assume(kappa>0) ;
assume(kappa0 >0) ;
assume(n>0) ;
Schrodinger: -diff(f,r,2) -2/r*diff(f,r) +l*(l+1)/r^2*f - kappa/r*f-epsilon*f = 0 ;
subst(f=(fw/r) , Schrodinger);
ev(%,diff) ;
Schrodinger_w : expand(%*r) ;
D(l,f):= diff(f,r)+l/r*f-kappa/(2*l)*f ;
U(l,f):= -diff(f,r)+l/r*f-kappa/(2*l)*f ;
U(l,D(l,f)) ;
[expand(U(l,D(l,f))),Schrodinger] ;
H(l,f):=(-kappa*f)/r+(l^2*f)/r^2+(l*f)/r^2-diff(f,r,2) ;
fwn : kappa^(n+1/2)*sqrt(n^((-2*n)-1)/gamma(2*n+1))*r^n*%e^-((kappa*r)/(2*n)) ;
subst([l=(n-1),f=fwn/r,epsilon=(-kappa^2/(4*n^2))],Schrodinger) ;
factor(expand(ev(%,diff))) ;

for n_:1 unless n_ > 10 do (dum:subst(n = n_,fwn),Rw[n_,n_-1]:dum,display(Rw[n_,n_-1]),for l_ : n_-1 step -1 unless l_ < 1 do (Rw[n_,l_-1]:factor(expand(D(l_,dum))),dum:Rw[n_,l_-1],display(Rw[n_,l_-1]))) $
for n:1 while n<=10 do ( for l:n-1 while l>=0 step -1 do ( dummy:integrate(Rw[n,l]^2,r,0,inf) , Rw[n,l]:Rw[n,l]/sqrt(dummy) ) ) ;
for n:1 while n<=10 do ( for l:n-1 while l>=0 step -1 do ( display(integrate(Rw[n,l]^2,r,0,inf) ) ) )$
disp("H(n-l,Rw[n,n-l]) = epsilon*Rw[n,n-l]") ;
for n unless n > 5 do ( for l:1 while l<=n do disp('Rw[n,n-l],ev(factor(expand(ev(H(n-l,Rw[n,n-l]) = epsilon*Rw[n,n-l],epsilon = -kappa^2/(4*n^2)))))) ) $

disp("probability vs. r") $
Const: [ ab = 5.292*10^(-11) , kappa = 1.203037e+10*%pi ] ;
ev(plot2d(subst(r=x*ab,[Rw[10,9]^2,Rw[10,8]^2,Rw[10,7]^2,Rw[10,6]^2,Rw[10,5]^2,Rw[10,4]^2,Rw[10,3]^2,Rw[10,2]^2,Rw[10,1]^2,Rw[10,0]^2 ]),[x,0,300]),Const ) ;


/* コメント化

コメント終了 */