利用python的matplotlib包绘制氢原子的径向分布函数

# to get the  angular momentum of a oneparticle system

import matplotlib.pyplot as plt
import numpy as np
from scipy import  integrate

# input the quantum number
# l = int(eval(input("please enter the quantum number l : ")))
# n = abs(int(eval(input('please enter the quantum number n : '))))
# z=int(input('please enter the atom number Z :'))
# m1,m2=input("please enter the molar mass of each atom:").split()
# m1,m2=(float(m1),float(m2))
l=1
n=3
z=1
m1=1.6726e-27
m2=9.109e-31
relamass=m1*m2/((m1+m2))
a=4*3.1415926*8.8542e-12*6.63*6.63e-68/(4*3.1415926*3.1415926*relamass*1.6*1.6e-38)
if n < 1+l:
    print("l should be less than n-1.")
    exit()

#to generate the wavefunction
jmax =n-l-1
b0=1
s1='0'

for j in range(0,jmax + 1):
    #print(j)
    exec(f'b{j + 1}=2*z*(j+l+1-n)/(n*a*(j+1)*(j+2*l+2))*b{j}')
    exec(f's1+="+"+str(b{j})+"*"+"x**{j}"')
#        exec(f'print(a{j+2})')
print(s1)

s2=f'x**{l}'
s3 = f'np.exp(-z*x/(n*a))'
s = '('+s1+')'+'*'+'('+s2+')'+'*'+s3
print(s)


#the radial function value for explicit x 
def radial(x):
    return eval(s)
#
def radialintegrate(x):
    return (eval(s)**2)*(x**2)
#normalization number
kk=integrate.quad(radialintegrate,0,1e-8)


#to draw the figure
plt.figure(figsize=(10, 8))
plt.title("the wave function picture of radial function.")
plt.xlim(0,5e-9)
x = np.linspace(0,5e-9, 100000)
y1 = radial(x)   #y1 is the radial function
#print(y1)
y2=radialintegrate(x)/kk[0]

max1 = max(abs(y2))
plt.xlabel("r", loc='left')
plt.ylim(-1.1*max1,  1.1*max1)
plt.ylabel('Wave function', loc='top')
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
# ax.xaxis.set_ticks_position('bottom')
# ax.yaxis.set_ticks_position('left')
ax.spines['bottom'].set_position(('data', 0))
ax.spines['left'].set_position(('data', 0))

l1, = plt.plot(x, y2, color='blue', linewidth=1.0, linestyle='-', label='down')
#plt.text(x=1e-6, y=0.9e-7, s=f'quantum number:{u}' + '\n' + 'E=(u+0.5)*h*v={:.4}'.format(
#    E) + '\n' + 'frequency={:.2f}(cm-1)\nmass={:.3}(g/mol)'.format(v, m))

plt.show()