h2p_solver.py 文件源码

def R_xi(E_vec,L0):

   traj = np.zeros([x_rad.shape[0]+1,E_vec.shape[0]])

   n=0
   for E in E_vec:

      c2 = D*D*E/4.

      L = newton(newton_ang_func,L0,args=(c2,m),tol=1e-8,maxiter=200)

      slope = -(-L+m*(m+1.)+2.*D+c2)/(2.*(m+1.))
      z0 = [1+step*slope,slope]

      z = odeint(g,z0,x_rad,args=(c2,L,D,m))

      temp=pow(x_rad,2.0)-1.
      temp=pow(temp,m/2.)

      zz=temp*z[:,0]

      first_zz = np.array([1])
      zz=np.append(first_zz, zz)

      traj[:,n] = zz

      n += 1

   xt = np.append(np.array([1]),x_rad)

   figure = plt.figure(figsize=(12, 11))
   plt.plot(xt,traj, linewidth=2.0, label = '')
   plt.ylabel('R($\\xi$)')#, fontdict=font)
   plt.xlabel('$\\xi$')#, fontdict=font)
   #plt.xlim(1.0,10.0)
   #plt.ylim(-1.0,1.0)
   plt.locator_params(axis='x', nbins=10)
   plt.locator_params(axis='y', nbins=10)
   plt.tick_params(axis='x', pad=15)
   #plt.legend(loc=1)

   plt.show()
   plt.close()