In [1]:
%matplotlib inline
import matplotlib 
#matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from astropy.io import fits
import numpy as np
import glob 

### this is an case of 3C286
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))
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)
        
        
x=np.linspace(1, 1024, num=1024)

Iarray = np.add(EX, EY)
Qarray = np.subtract(EX, EY)

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)


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,:]


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.)

allI = np.average(calI, axis=0)
allQ = np.average(calQ, axis=0)
allU = np.average(calU, axis=0)
allV = np.average(calV, axis=0)

## check cal intensity vary with the frequency
plt.plot(x, allI, 'r-', x, allQ, 'g-', x, allU, 'b-', x, allV,'k-')
#plt.ylim(-8.e7, 8.e7)
plt.xlim(0, 1024)
Out[1]:
(0, 1024)
In [2]:
#check cal signals vary with the frequency
parray = np.divide(np.power(np.add(np.power(allQ, 2), np.power(allU, 2)),0.5) , allI)
plt.plot(x, parray, 'r-')
plt.ylim(0, 2)
plt.xlim(0, 1024)

plt.xlabel('channel')
plt.ylabel('linear polarization degree')
Out[2]:
<matplotlib.text.Text at 0x10428c750>
In [3]:
## check rf signals (targets) linear polarization degrees vary with the frequency
rfI = np.divide(np.subtract(np.add(np.add(p3I, p4I),  np.add(p1I, p2I)), np.multiply(calI, 2.)), 4.)
rfQ = np.divide(np.subtract(np.add(np.add(p3Q, p4Q),  np.add(p1Q, p2Q)), np.multiply(calI, 2.)), 4.)
rfU = np.divide(np.subtract(np.add(np.add(p3U, p4U),  np.add(p1U, p2U)), np.multiply(calI, 2.)), 4.)
rfV = np.divide(np.subtract(np.add(np.add(p3V, p4V),  np.add(p1V, p2V)), np.multiply(calI, 2.)), 4.)


### all signals have been averaged.
allrfI = np.average(rfI, axis=0)
allrfQ = np.average(rfQ, axis=0)
allrfU = np.average(rfU, axis=0)
allrfV = np.average(rfV, axis=0)

rfparray = np.divide(np.power(np.add(np.power(allrfQ, 2), np.power(allrfU, 2)),0.5) , allrfI)


plt.plot(x, rfparray, 'r-')
plt.ylim(0, 1)
plt.xlim(0, 1024)

plt.xlabel('channel')
plt.ylabel('linear polarization degree')
Out[3]:
<matplotlib.text.Text at 0x10a4c2b90>
In [ ]: