记录“TDE”工作常用代码（光变分析）

asu552

已于 2023-11-15 14:14:50 修改

阅读量202

点赞数

文章标签： python 开发语言

于 2023-11-15 11:39:48 首次发布

本文链接：https://blog.csdn.net/asu552/article/details/134412837

版权

TDE光学光变工作，主要从光变数据的下载，处理作图，统计分析来展开。之后是其TDE-type的统计分析，以及光谱得到。下面是本人工作的步骤，并附上处理代码，以防止遗失和缺漏。

1.TDE光学数据的处理（以ZTF数据为例子，在irsa-ztf网址中，选取一个给定坐标一个角秒内的观测波段，以zr,zg,zi,将不同观测id的观测结果折叠在一起，得到对应波段的fits文件），实际上是对fits文件的处理。

######ztf 作图
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
import os
import math

ax = plt.gca()
ax.invert_yaxis() #反转y轴，因为星等跟光度反比
name='AT2021yzv' #TNS名称
num=59511.5 #peak处mjd
num1='59511.5'
path1 = 'R/lc(39).fits' #zr
path2 = 'G/lc(39).fits' #zg
#for file_name in os.listdir(path1):
    #print(file_name)

hdul1 = fits.open(path1) #
data1 = hdul1[1].data
chi1 = data1['chi']
catf1 = data1['catflags']
MAG1 = data1['mag'] 
MJD1 = data1['mjd']
MAGERR1 = data1['magerr']
r = []
mjd1 = []
yerr1 = []
yerr2 = []
yerr4 = []
l1 = len(MAG1)
v = 0

hdul2 = fits.open(path2)
data2 = hdul2[1].data
chi2 = data2['chi']
catf2 =data2['catflags']
MAG2 = data2['mag'] 
MJD2 = data2['mjd']
MAGERR2 = data2['magerr']
g = []
mjd2 = []
magerr1 = []
magerr2 = []
j = 0
l2 = len(MAG2)
for v in range (0, l1): #筛选条件 chi<4
    if chi1[v] < 4:
        if catf1[v]==0:
            r.append(MAG1[v])
            mjd1 = np.append(mjd1, MJD1[v]) 
            yerr1.append(MAGERR1[v])
            magerr1.append(MAGERR1[v])
            v=v+1


for j in range(0,l2):
    if chi2[j] < 4:
        if catf2[j]==0 :
            g.append(MAG2[j])
            mjd2 = np.append(mjd2,MJD2[j]) 
            yerr2.append(MAGERR2[j])
            magerr2.append(MAGERR2[j])
            j=j+1
            
print('done')
plt.xlabel('mjd-'+num1)
plt.ylabel('mag')
plt.errorbar(mjd1,r,yerr1,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='r',elinewidth=1,color='k',ms=2)
plt.errorbar(mjd2,g,yerr2,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='b',elinewidth=1,color='k',ms=2)
#plt.errorbar(mjd3,I,yerr4,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='m',elinewidth=1,color='k',ms=2)
label=['r','g']
plt.legend(label,loc=1)
plt.axvline(num,color="grey",linewidth=1, linestyle='--')
plt.title(name+' light curve')
plt.show()

gnew = []
rnew = []
mjd1new = []
mjd2new = []
rerr = []
gerr = []
yerr3 = []
m=0
z=0
for m in range (0, len(g)):  
    for z in range (0, len(r)-1):
        delta = (mjd1[z] - mjd2[m]) **2
        delta_1 = (mjd1[z-1] - mjd2[m]) **2
        delta1 = (mjd1[z+1] - mjd2[m]) **2
        if delta<0.25:  #筛选0.5天的两个波段数据，作g-r
            if delta_1>delta<delta1:
                mjd1new.append(mjd1[z])
                mjd2new.append(mjd2[m])
                rnew.append(r[z])
                gnew.append(g[m])
                rerr.append(magerr1[z])
                gerr.append(magerr2[m])
                break;
                z=z+1
            
            else:
                continue;
                m=m+1
n=0
g_r =[]
for n in range(0,len(rnew)):
    g__r = gnew[n] - rnew[n]
    g_r.append(g__r)
    yerr3.append(((rerr[n])**2 + (gerr[n])**2)**0.5)

#plt.plot(mjd2new,g_r,'b.')
#plt.plot(mjd2,g,'r.')
plt.errorbar(mjd2new,g_r,yerr3,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='m',elinewidth=1,color='k',ms=1.5)
#plt.xlim(58400,59000)
#plt.ylim(0.3,1)
plt.title(name+'  g-r')
plt.xlabel('mjd')
plt.ylabel('g-r')
plt.axvline(num,color="grey",linewidth=1,linestyle='--')
plt.show()

运行展示：

附：多波段的光学数据绘图：

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import leastsq
import os
from astropy.io import fits
import pandas as pd


#from astropy.io import fits
ax = plt.gca()
ax.invert_yaxis() 
path1 = 'R/lc(33).fits'
path2 = 'G/lc(31).fits'

filename='iPTF16axa light-curve'
f2 = np.genfromtxt("V-iPTF-16axa.csv", delimiter=",")
f1 = np.genfromtxt("iPTF-16axa.csv", delimiter="")

####ZTF
hdul1 = fits.open(path1)
data1 = hdul1[1].data
chi1 = data1['chi']
catf1 =data1['catflags']
MAG1 = data1['mag'] 
MJD1 = data1['mjd']
MAGERR1 = data1['magerr']
r = []
mjd1 = []
l1 = len(MAG1)
i = 0
hdul2 = fits.open(path2)
data2 = hdul2[1].data
chi2 = data2['chi']
catf2 =data2['catflags']
MAG2 = data2['mag'] 
MJD2 = data2['mjd']
MAGERR2 = data2['magerr']
g = []
mjd2 = []
magerr1 = []
magerr2 = []
j = 0
l2 = len(MAG2)
for i in range (0, l1):
    if chi1[i] < 4:
        if catf1[i]==0:
            r.append(MAG1[i])
            mjd1 = np.append(mjd1, MJD1[i]) 
            magerr1.append(MAGERR1[i])
            i=i+1

print(i)
for j in range (0, l2):
    if chi2[j] < 4:
        if catf2[j]==0:
            g.append(MAG2[j])
            mjd2 = np.append(mjd2, MJD2[j]) 
            magerr2.append(MAGERR2[j])
            j=j+1
print(j)
####

#print(f1.to_string())##to_string() 用于返回 DataFrame 类型的数据，避免中间部分以 ... 代替。

####ATLAS
#MJD = f1[:,0]
#Band = f1[:,5]
#Ra = f1[:,8]
#Dec = f1[:,9]
#MAG = f1[:,1]
#MAGerr = f1[:,2]
#uJy = f1[:,3]
#duJy = f1[:,4]

#a = 0
#mjd3=[]
#mjd4=[]
#mag3=[]
#yerr3=[]
#mag4=[]
#yerr4=[]
#l = len(MJD)
#for a in range(l):
    #if uJy[a]>=30:
        #if 0<=duJy[a]<=70:
            #if Band[a] ==0:
                #mjd3.append(MJD[a])
                #mag3.append(MAG[a])
                #yerr3.append(MAGerr[a])
            #if Band[a]==1:
                #mjd4.append(MJD[a])
                #mag4.append(MAG[a])
                #yerr4.append(MAGerr[a])
f1 = np.genfromtxt("iPTF-16axa-1day.csv",delimiter=",")  ##与上面相比，加了一天的bin
MJD = f1[:,0]
Band = f1[:,3]
uJy = f1[:,1]
duJy = f1[:,2]
num = f1[:,4]
a = 0

mjd3=[]
mjd4=[]
flux3=[]
yerr3=[]
flux4=[]
yerr4=[]
l = len(MJD)
for a in range(l): #限制条件
    if uJy[a]>=30:
        if 0<=duJy[a]<=70:
            if Band[a] ==0:
                mjd3.append(MJD[a])
                flux3.append(uJy[a])
                yerr3.append(duJy[a])
            if Band[a]==1:
                mjd4.append(MJD[a])
                flux4.append(uJy[a])
                yerr4.append(duJy[a])

mag_ao = -2.5 * np.log10(np.abs(flux3)) + 23.9
m_ao_err_up = -2.5 * np.log10(np.abs(flux3) + yerr3) + 23.9
m_ao_err_down = -2.5 * np.log10(np.abs(flux3) - yerr3) + 23.9
magerr_ao = np.abs(m_ao_err_up - m_ao_err_down) / 2

mag_ac = -2.5 * np.log10(np.abs(flux4)) + 23.9
m_ac_err_up = -2.5 * np.log10(np.abs(flux4)+yerr4) + 23.9
m_ac_err_down = -2.5 * np.log10(np.abs(flux4)-yerr4) + 23.9
magerr_ac = np.abs(m_ac_err_up - m_ac_err_down) / 2
                                
                
####

####CATALINA

VMAG = f2[:,1]
VMJD = f2[:,5]
VRA = f2[:,3]
VDEC = f2[:,4]
VMAGerr = f2[:,2]
Blend = f2[:,6]
b=0
mjd5 =[]
mag5=[]
yerr5=[]
l3 = len(VMAG)
for b in range(l3):
    if Blend[b]==0:
        mjd5.append(VMJD[b])
        mag5.append(VMAG[b])
        yerr5.append(VMAGerr[b])

###ASASSN
data = np.genfromtxt("AS15-oi.csv",delimiter=",")
umjd=data[:,0];
umag=data[:,1];
umagerr=data[:,2]
bmjd=data[:,3];
bmag=data[:,4];
bmagerr=data[:,5]
vmjd=data[:,6];
vmag=data[:,7];
vmagerr=data[:,8]

###WISE
f3 = np.genfromtxt("iPTF16axa-wise.csv",delimiter=",")
mjd6 = f3[:,0]
mag6 = f3[:,1]
magerr6 = f3[:,2]


###
r_2=[a-3 for a in r] 
g_2=[a-3 for a in g]
mag4_=[a+1 for a in mag4]
plt.errorbar(mjd1,r_2,magerr1,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='r',elinewidth=1,color='k',ms=1)#mjd1=r
plt.errorbar(mjd2,g_2,magerr2,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='b',elinewidth=1,color='k',ms=1)#mjd2=g
#plt.plot(mjd3,mag3,'y3')
#plt.plot(mjd4,mag4,'c4')
#plt.errorbar(mjd3,mag3,yerr3,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='orange',elinewidth=1,color='k',ms=1)#mjd3=o
#plt.errorbar(mjd4,mag4_,yerr4,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='green',elinewidth=1,color='k',ms=1)#mjd4=c
plt.errorbar(mjd3,mag_ao-1,magerr_ao,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='orange',elinewidth=1,color='orange',ms=1.5)
plt.errorbar(mjd4,mag_ac,magerr_ac,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='green',elinewidth=1,color='green',ms=1.5)
plt.errorbar(mjd5,mag5,yerr5,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='m',elinewidth=1,color='k',ms=1)#mjd5=v
plt.errorbar(mjd6,mag6+1,magerr6,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='r',elinewidth=1,color='k',ms=1)
#plt.errorbar(umjd,umag,umagerr,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='violet',elinewidth=1,color='k',ms=1)#U
#plt.errorbar(bmjd,bmag,bmagerr,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='turquoise',elinewidth=1,color='k',ms=1)#B
#plt.errorbar(vmjd,vmag,vmagerr,xerr=None,fmt='o',capthick=1,capsize=1.5,ecolor='springgreen',elinewidth=1,color='k',ms=1)#V


plt.title(filename)
#plt.ylim(23,17)
#plt.xlim(59000,60000)
plt.xlabel('mjd(57414)')
plt.ylabel('mag')
label=['r-3','g-3','o','c+1','v','W1+1']
plt.legend(label,loc=0)
plt.axvline(57414, color="grey", linewidth=1,  linestyle='--')
plt.show()

效果大概这样：