畫rotational-diagram

import numpy as np
import astropy as a
import astropy.units as u
import astropy.constants as c
from scipy import optimize
import matplotlib.pyplot as plt
plt.ion()

B = 24690.2
h = 6.626e-27     #erg/s
k = 1.38e-16      #erg/K
Tbg = 2.73        #K
nu = [356626.667, 356874.537, 345903.916, 345919.26, 344443.433, 344109.039, 342729.796, 343215.87]
mu = 0.899e-18
Eu = [727.82556, 717.88089, 332.6488, 459.43023, 451.22691, 419.39868, 227.47312, 1098.9015]
M = [63.246, 76.107, 181.02, 108.12, 146.69, 187.26, 232.87, 3.5619]
M_err = [1,1,1,1,1,0.25,1,1]
def f(x,p,t):
    return -x/t+p

y = []
for i in range(8):
    _y = np.log(3*k*M[i]*1e41/(8*12.943*nu[i]*(np.pi)**3))
    y.append(_y)
    
x = Eu
x_m = x[0:7]#[x[2],x[3],x[5]]#[x[2],x[3],x[5]]
y_m = y[0:7]#[y[2],y[3],y[5]]#[y[2],y[3],y[5]]
a,b = optimize.curve_fit(f, x_m, y_m, method='lm')
ylim = [30,50]
x1 = np.arange(200,1100)
y1 = -x1/a[1]+a[0]
plt.text(800,max(y)-0.4,"CH$_3$OH",fontsize=12)
plt.text(800,max(y)-0.8,"T$_{rot}$=%.2f (K)" %(a[1]),fontsize=12)
plt.plot(x1,y1)
plt.xlabel('E$_u$/k')
plt.ylabel(r'ln($\frac{3kW}{8 \pi^3 \nu S_{ij} \mu^2}$)')
column_density = np.exp(y[4]+Eu[4]/a[1])/1e19
plt.text(800,np.max(y)-1.2,"N=%.2f (10$^{19}$ cm$^{-2}$)" %(column_density),fontsize=12)
plt.scatter(x,y)
#plt.plot(x1,c='k')
#plt.text(x[0])
plt.savefig('fitting_1234567.pdf')