射电天文偏振线的绘制

#读取数据I是流量强度，PA是偏振角，PI是偏振强度
I_data =  pf.getdata('9000.I.fits')     
PA_hdr = pf.getheader('9000.PA.fits')
PI_hdr = pf.getheader('9000.PI.fits')
I_hdr = pf.getheader('9000.I.fits')

#print(PA_data.shape, PI_data.shape, I_data.shape)

Nx = 4   #这两个值为取值，越小对偏振的取值次数就会越频繁，反之亦然
Ny = 4
scale = 0.0007  #需要对PI进行比例的缩放，具体缩放程度看PI和I的数值差距
sigma = 9.89*10**-5 #根据自己数据的大小设置sigma，并且设置不绘制小于3倍sigma的值偏振信号
M0, N0 =PA_data.shape
#print(M0,N0)

PA_w = WCS(PA_hdr)
PI_w = WCS(PI_hdr)
I_w = WCS(I_hdr)

for i in range(0,M0,Nx):
    for j in range(0,N0,Ny):
    	#没有偏振角或者偏振小于sigma
        if np.isnan(PA_data[i,j]) or PI_data[i,j]< sigma: continue    
        pa0 = np.radians(PA_data[i, j]) + np.pi/2.
        pi0 = np.sqrt(PI_data[i, j]* scale) / 2.

        l, b = PA_w.wcs_pix2world([[j+1,i+1,0]], 0)[0][0:2]
        y1 = b+pi0*np.cos(pa0)
        x1 = l+pi0*np.sin(pa0)/np.cos(np.radians(y1))
        y2 = b-pi0*np.cos(pa0)
        x2 = l-pi0*np.sin(pa0)/np.cos(np.radians(y2))

        px1, py1 = I_w.wcs_world2pix([[x1, y1, 0, 0]], 0)[0][0:2]
        px2, py2 = I_w.wcs_world2pix([[x2, y2, 0, 0]], 0)[0][0:2]
        plt.plot([px1,px2],[py1,py2],'green') 
        
'''
#画偏振的箭头指向，拿的别人的代码，还没有尝试过，先这里吧
ang0 = 252.
d_ang = 59.2208573466
ang = np.radians(ang0 + d_ang)
x1 = 68.77
y1 = 2.82
y2 = y1 + 0.5*np.cos(ang)
x2 = x1 + 0.5*np.sin(ang)/np.cos(np.radians(y2))
px1, py1 = I_data.convproj.topixel((x1,y1))
px2, py2 = I_data.convproj.topixel((x2,y2))
plt.arrow(px1, py1, px2-px1, py2-py1, width=3., fc='c')
'''

#plt.show()