In [3]:
### This is a script to check all data and bandpass ###
%matplotlib inline
import matplotlib
#matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from astropy.io import fits
import numpy as np
import glob
###### combine all data in a single file ########
myfile = glob.glob("EFFBG_2017-12-06_49*/*/S45mm-SPECPOL-ARRAYDATA*fits")
EX= np.empty(shape=(0,1024)) # already squared
EY= np.empty(shape=(0,1024)) # already squared
Uarray=np.empty(shape=(0,1024))
Varray=np.empty(shape=(0,1024))
obst = np.empty(shape=(0))
for i in range(np.size(myfile)):
tbdata = fits.getdata(myfile[i])
EX = np.concatenate((tbdata['data'][:,0,:], EX), axis=0)
EY = np.concatenate((tbdata['data'][:,1,:], EY), axis=0)
Uarray = np.concatenate((tbdata['data'][:,2,:], Uarray), axis=0)
Varray = np.concatenate((tbdata['data'][:,3,:], Varray), axis=0)
obst = np.concatenate((tbdata['MJD'], obst), axis=0)
#hdulist = fits.open('S45mm-SPECPOL-ARRAYDATA-1.fits')
#hdulist.info()
#hdulist[0].header
#tbdata = hdulist[1].data
x=np.linspace(1, 1024, num=1024)
#hdulist[1].header['FREQRES']
#EX = tbdata['data'][:,0,:]
#EY = tbdata['data'][:,1,:]
### I = EX^2 + EY^2
### Q = EX^2- EY^2
#Iarray = np.add(EX, EY)
#Qarray = np.subtract(EX, EY)
Iarray = EX + EY
Qarray = EX - EY
#Iarray = tbdata['data'][:,0,:]
#Qarray = tbdata['data'][:,1,:]
#Uarray = tbdata['data'][:,2,:]
#Varray = tbdata['data'][:,3,:]
### define four phase for the calibration and the time
arz = Iarray.shape[0]/4
p1I=np.ndarray(shape=(arz,1024), dtype=float)
p2I=np.ndarray(shape=(arz,1024), dtype=float)
p3I=np.ndarray(shape=(arz,1024), dtype=float)
p4I=np.ndarray(shape=(arz,1024), dtype=float)
p1Q=np.ndarray(shape=(arz,1024), dtype=float)
p2Q=np.ndarray(shape=(arz,1024), dtype=float)
p3Q=np.ndarray(shape=(arz,1024), dtype=float)
p4Q=np.ndarray(shape=(arz,1024), dtype=float)
p1U=np.ndarray(shape=(arz,1024), dtype=float)
p2U=np.ndarray(shape=(arz,1024), dtype=float)
p3U=np.ndarray(shape=(arz,1024), dtype=float)
p4U=np.ndarray(shape=(arz,1024), dtype=float)
p1V=np.ndarray(shape=(arz,1024), dtype=float)
p2V=np.ndarray(shape=(arz,1024), dtype=float)
p3V=np.ndarray(shape=(arz,1024), dtype=float)
p4V=np.ndarray(shape=(arz,1024), dtype=float)
phtim = np.ndarray(shape=(arz), dtype=float)
for i in range(arz):
p1I[i,:] = Iarray[4*i,:]
p2I[i,:] = Iarray[4*i+1,:]
p3I[i,:] = Iarray[4*i+2,:]
p4I[i,:] = Iarray[4*i+3,:]
p1Q[i,:] = Qarray[4*i,:]
p2Q[i,:] = Qarray[4*i+1,:]
p3Q[i,:] = Qarray[4*i+2,:]
p4Q[i,:] = Qarray[4*i+3,:]
p1U[i,:] = Uarray[4*i,:]
p2U[i,:] = Uarray[4*i+1,:]
p3U[i,:] = Uarray[4*i+2,:]
p4U[i,:] = Uarray[4*i+3,:]
p1V[i,:] = Varray[4*i,:]
p2V[i,:] = Varray[4*i+1,:]
p3V[i,:] = Varray[4*i+2,:]
p4V[i,:] = Varray[4*i+3,:]
phtim[i] = np.average(obst[4*i:4*i+3])
calI = (p3I + p4I - (p1I+p2I))/2
calQ = (p3Q + p4Q - (p1Q+p2Q))/2
calU = (p3U + p4U - (p1U+p2U))/2
calV = (p3V + p4V - (p1V+p2V))/2
#calI = np.divide(np.subtract(np.add(p3I, p4I), np.add(p1I, p2I)), 2.)
#calQ = np.divide(np.subtract(np.add(p3Q, p4Q), np.add(p1Q, p2Q)), 2.)
#calU = np.divide(np.subtract(np.add(p3U, p4U), np.add(p1U, p2U)), 2.)
#calV = np.divide(np.subtract(np.add(p3V, p4V), np.add(p1V, p2V)), 2.)
#### cal I, Q, U, V
############### This shows all calibrated data (input signals) in a glimpse ##################
fig1, ax = plt.subplots(2,2)
ax[0,0].imshow(calI, aspect =1./6, cmap='jet',vmin=calI.min()/1e4*0, vmax=calI.max()/5e4)
ax[0,0].text(900, 500,"I", color='w')
#ax1.set_xlabel('channel')
#ax1.set_ylabel('time')
ax[0,1].imshow(calQ, aspect =1./6, cmap='jet',vmin=calQ.min()/1e4*0, vmax=calQ.max()/5e4)
ax[0,1].text(900, 500,"Q", color='w')
ax[1,0].imshow(calU, aspect =1./6, cmap='jet',vmin=calU.min()/1e4*0, vmax=calU.max()/5e4)
ax[1,0].text(900, 500,"U", color='w')
ax[1,1].imshow(calV, aspect =1./6, cmap='jet',vmin=calV.min()/1e4*0, vmax=calV.max()/5e4)
ax[1,1].text(900, 500,"V", color='w')
fig1.text(0.5, 0.04, 'Channel (cal-signal)', ha='center')
fig1.text(0.04, 0.5, 'Time', va='center', rotation='vertical')
Out[3]:
In [2]:
rfI = (p3I+p4I + p1I +p2I- calI*2)/4
rfQ = (p3Q+p4Q + p1Q +p2Q- calQ*2)/4
rfU = (p3U+p4U + p1U +p2U- calU*2)/4
rfV = (p3V+p4V + p1V +p2V- calV*2)/4
############### This shows all observed data (targets) in a glimpse ##################
fig2, axes = plt.subplots(2,2)
cax1 = axes[0,0].imshow(rfI, aspect =1./6, cmap='jet',vmin=0, vmax=rfI.max()/1e3)
#cbar = fig.colorbar(cax = cbaxes1, ax = axes[0,0], fraction =0.046, pad =0.04, orientation='horizontal')
axes[0,0].text(900, 500,"I", color='w')
cax2 = axes[0,1].imshow(rfQ, aspect =1./6,cmap='jet',vmin=rfQ.min()/1e4, vmax=rfQ.max()/1e4)
#fig.colorbar(cax2, ax = axes[0,1])
axes[0,1].text(900, 500,"Q", color='w')
cax3 = axes[1,0].imshow(rfU, aspect =1./6,cmap='jet',vmin=rfU.min()/1e4, vmax=rfU.max()/1e4)
axes[1,0].text(900, 500,"U", color='k')
cax4 = axes[1,1].imshow(rfV, aspect =1./6,cmap='jet',vmin=rfV.min()/1e4, vmax=rfV.max()/1e4)
axes[1,1].text(900, 500,"V", color='k')
fig2.text(0.5, 0.04, 'Channel (sky)', ha='center')
fig2.text(0.04, 0.5, 'Time', va='center', rotation='vertical')
Out[2]:
In [ ]: