ridge/scripts/01_crabmodel_plot_poly.py

79 lines
2.3 KiB
Python
Executable File

#!/usr/bin/env python
__author__ = "Roman Krivonos"
__copyright__ = "Space Research Institute (IKI)"
from astropy.table import Table, Column
import matplotlib.pyplot as plt
import sys, os
from ridge.utils import *
from ridge.config import *
scale = 1e-3
fn="detcnts.A01.crabmodel.fits"
dat = Table.read(proddir+fn)
df1 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.A02.crabmodel.fits"
dat = Table.read(proddir+fn)
df2 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.A03.crabmodel.fits"
dat = Table.read(proddir+fn)
df3 = dat.to_pandas().sort_values(by=['REV'])
fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True, figsize=(9, 7), dpi=100)
#fig.suptitle('Vertically stacked subplots')
#plt.figure(figsize=(8, 6), dpi=80)
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(1)
ax2.spines[axis].set_linewidth(1)
ax3.spines[axis].set_linewidth(1)
ax1.tick_params(axis="both", width=1, labelsize=14)
ax2.tick_params(axis="both", width=1, labelsize=14)
ax3.tick_params(axis="both", width=1, labelsize=14)
ax1.set_title("25-60 keV")
ax2.set_title("60-80 keV")
ax3.set_title("80-200 keV")
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax3.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax1.errorbar(df1['REV'], df1['B_EST']/scale, yerr=df1['ERR']/scale, fmt='o' )
ax1.plot(df1['REV'], df1['B_POLY']/scale, color='r', linewidth=2)
ax1.grid(visible=True)
ax2.errorbar(df2['REV'], df2['B_EST']/scale, yerr=df2['ERR']/scale, fmt='o' )
ax2.plot(df2['REV'], df2['B_POLY']/scale, color='r', linewidth=2)
ax2.grid(visible=True)
ax3.errorbar(df3['REV'], df3['B_EST']/scale, yerr=df3['ERR']/scale, fmt='o' )
ax3.plot(df3['REV'], df3['B_POLY']/scale, color='r', linewidth=2)
ax3.grid(visible=True)
#x = np.arange(-10, 10, 0.001)
#plot normal distribution with mean 0 and standard deviation 1
#plt.plot(x, norm.pdf(x, 0, 1), color='red', linewidth=2)
plt.xlabel('Revolution',fontsize=14, fontweight='normal')
ax2.set_ylabel('Count rate, x10$^{-3}$ cts s$^{-1}$ pix$^{-1}$',fontsize=14, fontweight='normal')
#plt.xscale('linear')
#plt.yscale('linear')
if not os.path.exists(figdir):
os.makedirs(figdir)
filename=figdir+'crabmodel_poly.png'
plt.savefig(filename, bbox_inches='tight')
plt.close(fig)
print("Result is saved as {}".format(filename))