TT-plots 从温度定标获得区域内的平均谱指数

针对某个源内得多个pixe点。在两个频率上分别对应的温度数据
可以做他们的TT-plots得到的斜率加上2即为流量定标下的谱指数
提供了一种其他方式计算谱指数，平均的好处是可以剔除某个不好的点，当然坏处也是，如果谱指数不均匀，可能带来比较大的偏差。

#!/usr/bin/env python
from astropy.io import fits as pf
import numpy as np
import os
import sys
from scipy import odr
from scipy import signal
from astropy.modeling import models, fitting
import matplotlib.pyplot as plt

def cal_alpha(data1,data2):
    n = data1.size
    x = np.zeros(n) 
    y = np.zeros(n)
    kk = 0
    for i in range(n):
        if np.isnan(data1[i]) or np.isnan(data2[i]): continue
        x[kk] = data1[i]
        y[kk] = data2[i]
        kk += 1

    if kk<10: return np.nan,np.nan,np.nan,np.nan

    p_init = models.Polynomial1D(1)
    fit_p = fitting.LinearLSQFitter()
    p = fit_p(p_init, x[:kk], y[:kk])
    a, b = p.c0.value, p.c1.value
    return a, b


def main():
	
	#设定基础频率
    freq0 = 1000.00357627869
    df = 0.00762939453125
    chan1 = 4258
    chan2 = 60916
    #前面四行参数可以不用管，直接设定的基础频率就可以了
    freq1 = (freq0 + chan1 * df)/1.e3  # in GHz   
    freq2 = (freq0 + chan2 * df)/1.e3  # in GHz
	
	#对指定区域进行拟合，例如本文就是设定mask大于0.5区域是目标源
    maskfile = 'mask.fits'
    mask = pf.getdata(maskfile)  

	#导入两个频率下的温度数据
    data1 = pf.getdata(inf1.fits)[mask>0.5]
    data2 = pf.getdata(inf2.fits)[mask>0.5]
	
	#以下不需要修改
    a1, b1 = cal_alpha(data1, data2)
    a2, b2 = cal_alpha(data2, data1)
    alpha1 = np.log10(b1)/np.log10(freq2/freq1)
    alpha2 = np.log10(b2)/np.log10(freq1/freq2)
    print(alpha1, alpha2)
    alpha_1 = (alpha1 + alpha2) / 2.
    d_alpha_1 = abs(alpha2 - alpha1)
    a1_p = a1 * 1.
    b1_p = b1 * 1.
    a2_p = -a2/b2
    b2_p = 1./b2
	
	#画图
    fig = plt.figure(figsize=(15,10))
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(data1, data2, 'o')
    ax.plot(data1, a1_p + b1_p * data1)
    ax.plot(data1, a2_p + b2_p * data1)
    ss = r'$-$' + '%4.2f' % abs(alpha_1) + r'$\pm$' + '%4.2f' % d_alpha_1
    ax.text(0.1, 0.9, ss, horizontalalignment='left', verticalalignment='center', transform=ax.transAxes)
    ax.set_xlabel('T (K, at %5.3f GHz)' % freq1)
    ax.set_ylabel('T (K, at %5.3f GHz)' % freq2)
        
    plt.show()
    #plt.savefig('tt.pdf', bbox_inches='tight')
	#最后得到的斜率加上2即为谱指数
	
if __name__ == "__main__":
    main()