import matplotlib.pyplot as plt
import os
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astropy.io import fits
import warnings
warnings.filterwarnings('ignore')
T = Table.read('atoms_B3_12M_hotcore_candidate_modified.tbl',format='ascii')
#spws = ['spw23','spw25','spw27','spw29','spw31','spw33','spw35','spw37','spw39']
lines = ['H$^{13}$CN','H$^{13}$CO$^+$','CCH','SiO','HCN','HCO$^+$']
freqs = [86339.9,86754.3,87316.9,86846.96,88631.847,89188.526]
def seek_file(path,sourcename):
files = []
for filename in os.listdir(path):
fp = os.path.join(filename)
fmts = '.txt'
if os.path.isfile(path+'/'+fp) and fmts in filename:
files.append(fp)
files.sort()
return files
def show_spec(path,sourcename,ID,files,ra,dec):
a = 0
coord = SkyCoord(ra,dec,unit='deg',frame='fk5')
hh,hm,hs = coord.ra.hms
dd,dm,ds = coord.dec.dms
dm,ds = abs(dm),abs(ds)
for j in range(len(files)):
filename = files[j]
x,y = np.loadtxt(path+'/'+filename,usecols=(0,1),unpack=True)
fitsfile = str(filename)[:-4]+'.fits'
hdu = fits.open(path+'/'+fitsfile)
windows = hdu[0].header['FILNAM04']
plt.subplots_adjust(hspace=0.3)
if len(x)>2000:
a = a+1
ax1 = plt.subplot(5,2,a)
ax1.tick_params(labelsize=4,direction='in', length=1, width=0.2,pad=0.7)
plt.title('%s_%s_%s' %(sourcename,ID,windows),loc='center',fontsize=4,pad=0.7)
if a==9:
plt.xlabel('Frequency (MHz)',fontsize=4)
plt.ylabel('Temperature (K)',fontsize=4)
a = a+1
ax2 = plt.subplot(5,2,a)
ax2.tick_params(labelsize=4,direction='in', length=1, width=0.2,pad=0.7)
ax1.plot(x[:int(len(x)/2)],y[:int(len(x)/2)],c='k',linewidth=0.2)
ax1.set_xlim(xmin=np.nanmin(x[:int(len(x)/2)]),xmax=np.nanmax(x[:int(len(x)/2)]))
ax2.plot(x[int(len(x)/2):],y[int(len(x)/2):],c='k',linewidth=0.2)
ax2.set_xlim(xmin=np.nanmin(x[int(len(x)/2):]),xmax=np.nanmax(x[int(len(x)/2):]))
if a==10:
plt.xlabel('Frequency (MHz)',fontsize=4)
plt.title('%s_%s_%s' %(sourcename,ID,windows),loc='center',fontsize=4,pad=0.7)
plt.ylabel('Temperature (K)',fontsize=4)
else:
a = a+1
ax = plt.subplot(5,2,a)
ax.tick_params(labelsize=4,direction='in', length=1, width=0.2,pad=0.7)
ax.plot(x,y,c='k',linewidth=0.2)
ax.set_xlim(xmin=np.nanmin(x),xmax=np.nanmax(x))
#plt.xlabel('Frequency (MHz)',fontsize=2)
plt.ylabel('Temperature (K)',fontsize=4)
for n in range(len(lines)):
if np.nanmin(x)<freqs[n]<np.nanmax(x):
plt.title('%s_%s_%s %s' %(sourcename,ID,windows,lines[n]),loc='center',fontsize=4,pad=0.7)
#plt.tight_layout()
plt.savefig('spec/figure/%s_%s.pdf' %(sourcename,ID))
plt.close('all')
for i in range(len(T['sname'])):
path = 'spec/'+T['sname'][i]+'/'+str(T['ID'][i])
files = seek_file(path,T['sname'][i])
show_spec(path,T['sname'][i],str(T['ID'][i]),files,T['ra'][i],T['dec'][i])