from scipy import * from scipy import weave from scipy import integrate from scipy import optimize from pylab import * code_Numerov=""" // void Numerov(const container& F, int Nmax, double dh, container& Solution) double dx = dh; double h2 = dx*dx; double h12 = h2/12; double w0 = (1-h12*F(0))*Solution(0); double Fx = F(1); double w1 = (1-h12*Fx)*Solution(1); double Phi = Solution(1); double w2; for (int i=2; i=(n_max-l): break u0=u1 dE = dE/1.091 Eb += Ebl E0 = Ebl[0][1] Eb.sort(comp) return Eb Z_screened=1. Z=82 n_max=5 lmax=3 #R = logspace(-7,2.3,1000) #R0 = linspace(1e-7,50,2**10+1) R0 = linspace(1e-7,90,2**12+1) R=R0[::-1] Eb = FindBoundStates(Z_screened,R,n_max,lmax) rho = zeros(len(R),dtype=float) Nelec=0 for i,(l,Ene) in enumerate(Eb): print 'state=', l, Ene u = ComputeSchrod(Ene,R,l,Z_screened,True) dN = 2*(2*l+1) if Z >= Nelec+dN: ferm=1. else: ferm= (Z-Nelec)/(dN+0.0) Nelec += dN*ferm rho += u*u*dN*ferm / (4*pi*R**2) if Nelec>=Z: break print 'Adding ', dN*ferm, 'electrons of l=', l, 'Ene=', Ene print 'Nelec=', Nelec print 'Norm=', integrate.romb(rho*4*pi*R**2,R[0]-R[1]) plot(R,rho*4*pi*R**2, 'o-') show()