#! /usr/bin/env python3 # -*- coding: iso-8859-15 -*- from glob import glob import numpy as np from numpy import pi,linspace,sqrt,append,arange,array,zeros,ones,dot,angle,exp,conj,cos,sin from astropy.io import fits from astropy.time import Time as astroTime import time from scipy.linalg import pinv as pinv2 badfilesList=[ "2018-09-22T07:56:25",# beta pictoris "2019-09-12T05:47:38",# beta pictoris "2020-01-07T02:07:02",#: "Beta Pictorisb", "2020-01-07T02:11:20",#: "Beta Pictorisb", "2020-01-07T02:29:01",#: "Beta Pictorisb", "2020-01-07T02:32:33",#: "Beta Pictorisb", "2020-01-07T02:46:50",#: "Beta Pictorisb", "2020-01-07T02:49:20",#: "Beta Pictorisb", "2020-01-07T02:52:05",#: "Beta Pictorisb", "2020-01-07T03:07:10",#: "Beta Pictorisb", "2020-01-07T03:09:22",#: "Beta Pictorisb", "2020-01-07T03:29:16",#: "Beta Pictorisb", "2020-01-07T03:31:00",#: "Beta Pictorisb", "2020-01-07T03:37:22",#: "Beta Pictorisb", "2020-01-07T03:45:50",#: "Beta Pictorisb", "2020-01-07T03:56:14",#: "Beta Pictorisb", "2020-01-07T03:58:25",#: "Beta Pictorisb", "2020-01-07T04:05:28",#: "Beta Pictorisb", "2020-01-07T04:09:28",#: "Beta Pictorisb", "2020-01-07T04:15:35",#: "Beta Pictorisb", "2020-01-07T04:16:45",#: "Beta Pictorisb", "2020-01-07T04:22:56",#: "Beta Pictorisb", "2020-01-07T04:26:44",#: "Beta Pictorisb", "2020-01-07T04:33:12",#: "Beta Pictorisb", "2024-02-11T05:25:34", # HD 115470: "2024-02-11T05:40:59", # HD 115470: "2024-03-19T04:13:09", # HR 2562: "2024-06-04T07:37:32", # HD206 "2024-05-31T09:39:45", # HD206 "2024-05-31T09:46:59", # HD206 "2024-06-04T07:37:32", # HD206 "2024-06-02T00:59:02", # HD135344B "2024-06-02T01:03:55",# HD135344B "2024-06-02T01:10:44", # HD135344B ] Bad_Tel_1=[] Bad_Tel_2=[ "2023-05-08T06:33:30",#StKM 1-1494 "2023-05-08T06:40:48",#StKM 1-1494 "2023-05-08T06:48:14",#StKM 1-1494 "2023-05-08T06:50:21",#StKM 1-1494 "2023-05-08T06:57:38",#StKM 1-1494 "2023-05-08T07:12:40",#StKM 1-1494 "2023-05-08T07:19:57",#StKM 1-1494 "2023-05-08T07:27:25",#StKM 1-1494 ] Bad_Tel_3=[ "2022-10-24T01:12:26", #HD984 ] Bad_Tel_4=[] class bTextColors: HEADER = '\033[95m' OKBLUE = '\033[94m' OKGREEN = '\033[92m' WARNING = '\033[93m' FAIL = '\033[91m' ENDC = '\033[0m' BOLD = '\033[1m' UNDERLINE = '\033[4m' # Code Mathias def calculatePhaseRefArclength(phaseref_coeffs, wav): """ Calculate the arclength for the polynomial fit to phaseRef args: phaseref_coeffs: array of size (ndit, nchannel, 4) containing the coeffs of the poly fit to phase ref wav: wavelength axis (in m) """ # calculate normalized wavelength grid used to fit phaseRef w0 = wav.mean() wref = (w0/wav - 1)/(np.max(wav) - np.min(wav))*w0 # prepare array for values ndit = np.shape(phaseref_coeffs)[0] nchannel = np.shape(phaseref_coeffs)[1] phaseRefArclength = np.zeros([ndit, nchannel]) # calculate arclength for dit in range(ndit): for c in range(nchannel): coeffs = phaseref_coeffs[dit, c, :] # calculate the coeffs of 1+(dy/dx)**2 arclengthpolycoeff = np.array([1+coeffs[1]**2, 4*coeffs[1]*coeffs[2], 4*coeffs[2]**2+6*coeffs[1]*coeffs[3], 12*coeffs[2]*coeffs[3], 9*coeffs[3]**2]) # integrate sqrt(1+(dx/dy)**2) do get arclength. Minus sign because wref is decreasing phaseRefArclength[dit, c] = -np.trapz(np.sqrt(np.polyval(arclengthpolycoeff[::-1], wref)), wref) return phaseRefArclength def angle_unwrap(C): Cmean=C.sum(axis=0) P=angle(C*conj(Cmean)) return P+angle(Cmean) def factorize_number(n): """Factorise un nombre en nombres premiers""" i = 2 factors = [] while i * i <= n: if n % i: i += 1 else: n //= i factors.append(i) if n > 1: factors.append(n) return factors class GravDataSet: """GRAVITY class""" def __init__(self,file,fheader,ext): self.file=file self.ext=ext self.fheader=fheader self.swap=0 self.swap_ratio=1 self.centred=False self.target=0 self.date=fheader["DATE-OBS"] # self.dit=fheader["HIERARCH ESO DET2 SEQ1 DIT"] self.ndit_compact=fheader["HIERARCH ESO DET2 NDIT"] # self.ndit_compact=self.ndit self.ftstart=astroTime([fheader["HIERARCH ESO PCR ACQ START"]],format='isot',scale='utc').mjd[0]*24*3600 self.Ref_Amp=[] self.Ref_Phase=[] self.ext_day=str(ext)+str(int(self.ftstart/24/3600+0.5)) self.res_axis_pola=fheader["HIERARCH ESO INS SPEC RES"]+fheader["HIERARCH ESO INS OPTI11 ID"]+fheader["HIERARCH ESO INS POLA MODE"] try: self.mjd=(fheader["MJD-OBS"]+fheader["MJD-END"])/2 except: self.mjd=fheader["MJD-OBS"] self.wave=fits.getdata(self.file,'OI_WAVELENGTH',self.ext).field('EFF_WAVE') h_v={} h_v["DATE"]=fheader["DATE-OBS"] h_v["DIT"]=fheader["HIERARCH ESO DET2 SEQ1 DIT"] h_v["GAIN"]=fheader["HIERARCH ESO INS DET2 GAIN"] h_v["DITFT"]=fheader["HIERARCH ESO DET3 SEQ1 DIT"] h_v["GAINFT"]=fheader["HIERARCH ESO INS DET3 GAIN"] h_v["NDIT"]=fheader["HIERARCH ESO DET2 NDIT"] h_v["RES"]=fheader["HIERARCH ESO INS SPEC RES"] h_v["POLA"]=fheader["HIERARCH ESO INS POLA MODE"] h_v["AXIS"]=fheader["HIERARCH ESO INS OPTI11 ID"] h_v["PROG_ID"]=fheader["HIERARCH ESO OBS PROG ID"] h_v["SEPTRK"]=fheader["HIERARCH ESO FT KAL SEPTRK"] h_v["SXY_CTU"]=array([(fheader["HIERARCH ESO FT KAL SEPCTUX%i"%i],fheader["HIERARCH ESO FT KAL SEPCTUY%i"%i]) for i in 1+arange(4)]) h_v["MET_ROW"]=fheader["HIERARCH ESO QC MET REF ROW"] h_v["MET_DIFF"]=array([fheader["HIERARCH ESO QC MET DRS DIFF"+str(i+1)] for i in range(4)]) h_v["REJ_RATIO"]=array([fheader["HIERARCH ESO QC REJECTED RATIO FT"],fheader["HIERARCH ESO QC REJECTED RATIO SC"]]) try: h_v["SXY"]=array([fheader["HIERARCH ESO INS SOBJ X"],fheader["HIERARCH ESO INS SOBJ Y"]]) except: h_v["SXY"]=array([0.0,0.0]) try: h_v["NAME"]=fheader["HIERARCH ESO OBS TARG NAME"] h_v["NAME_SC"]=fheader["HIERARCH ESO INS SOBJ NAME"] h_v["NAME_FT"]=fheader["HIERARCH ESO FT ROBJ NAME"] h_v["RA"]=fheader["RA"] h_v["DEC"]=fheader["DEC"] except: h_v["NAME"]="Noname" h_v["NAME_SC"]="Noname_SC" h_v["NAME_FT"]="Noname_FT" h_v["RA"]=0.0 h_v["DEC"]=0.0 try: h_v["UT"]=fheader["HIERARCH ESO ISS CONF T1NAME"] h_v["STATIONS"]=fheader["HIERARCH ESO ISS CONF STATION1"]+fheader["HIERARCH ESO ISS CONF STATION2"]+fheader["HIERARCH ESO ISS CONF STATION3"]+fheader["HIERARCH ESO ISS CONF STATION4"] h_v["AIRM_START"]=fheader["HIERARCH ESO ISS AIRM START"] h_v["SEEING_START"]=fheader["HIERARCH ESO ISS AMBI FWHM START"] h_v["TAU0_START"]=fheader["HIERARCH ESO ISS AMBI TAU0 START"] h_v["AIRM_END"]=fheader["HIERARCH ESO ISS AIRM END"] h_v["SEEING_END"]=fheader["HIERARCH ESO ISS AMBI FWHM END"] h_v["TAU0_END"]=fheader["HIERARCH ESO ISS AMBI TAU0 END"] h_v["DATE_END"]=fheader["HIERARCH ESO PCR ACQ END"] except: h_v["UT"]="None" h_v["STATIONS"]="None" h_v["AIRM_START"]=0 h_v["SEEING_START"]=0 h_v["TAU0_START"]=0 h_v["AIRM_END"]=0 h_v["SEEING_END"]=0 h_v["TAU0_END"]=0 h_v["DATE_END"]="None" if np.any(np.array(badfilesList)==self.date): h_v["flag"]=True else: h_v["flag"]=False self.h_v=h_v def get_xy_and_comment(self,verbose): h_v=self.h_v self.sxySCFTName=str(h_v['SXY'][0])+"_"+str(h_v['SXY'][1])+h_v['NAME_FT'] self.sxyFTSCName=str(-h_v['SXY'][0])+"_"+str(-h_v['SXY'][1])+h_v['NAME_SC'] self.sxySCFTName+=h_v["POLA"]+h_v["AXIS"]+h_v["RES"]+self.ext_day self.sxyFTSCName+=h_v["POLA"]+h_v["AXIS"]+h_v["RES"]+self.ext_day # np.array([str(x)+"_"+str(y) for x,y in h_v['SXY']]) # self.sxySCFT_name=xy+name for xy,name in zip(sxySCFT,h_v['NAME_FT']) # sxyFTSC=np.array([str(-x)+"_"+str(-y) for x,y in h_v['SXY']]) # self.sxyFTSC_name=np.array([xy+name for xy,name in zip(sxyFTSC,h_v['NAME_SC'])]) # sxySCFT[h_v['sep']<1]='None' if h_v["SEPTRK"]>0.5: self.sxy=h_v["SXY"]+0 else: self.sxy=h_v["SXY_CTU"].mean(axis=0) self.sep=np.sqrt(self.sxy[0]**2+self.sxy[1]**2) if self.sep<10: self.centred=True self.sxySCFTName='None' h_v['comment']=[] h_v['comment_color']=[] data,dit,ndit,red,pola,axis,sxy,name = (h_v["DATE"], h_v["DIT"], h_v["NDIT"], h_v["RES"], h_v["POLA"], h_v["AXIS"], h_v["SXY"], h_v["NAME"]) h_v['comment']=data+" "+str(dit)+" "+str(ndit)+" "+red+" "+pola+" "+axis+" "+str(sxy)+" "+name h_v['comment_color']=h_v['comment'] if h_v["flag"]: h_v['comment_color']=bTextColors.FAIL+h_v['comment_color']+bTextColors.ENDC h_v['verification']=[] h_v['verification_color']=[] data,rej, met = (h_v["DATE"],h_v["REJ_RATIO"],h_v["MET_DIFF"]) string_rej="reject FT=%0.2i%%,SC=%0.2i%%"%(rej[0],rej[1]) string_rej_color=string_rej if rej[1]>50: string_rej_color=bTextColors.FAIL+string_rej+bTextColors.ENDC elif rej[0]>20: string_rej_color=bTextColors.WARNING+string_rej+bTextColors.ENDC string_met=" DRS/TAC-MET=%.3f,%.3f,%.3f,%.3f"%(met[0],met[1],met[2],met[3]) string_met_color=string_met if max(met)>0.1: string_met_color=bTextColors.FAIL+string_met+bTextColors.ENDC elif max(met)>0.01: string_met_color=bTextColors.WARNING+string_met+bTextColors.ENDC h_v['verification']+=[data+" "+string_rej+string_met] h_v['verification_color']+=[data+" "+string_rej_color+string_met_color] if verbose==True: print(h_v["comment_color"]) print(h_v["verification_color"]) def read_oi_flux(self): oi_flux_data={ "TIME":[], 'FLUX':[], 'FDDL_FT':[], 'FDDL_SC':[], 'PHASE_MET_FC':[], 'PHASE_MET_FCFT':[], 'PHASE_MET_FCSC':[], } N=self.h_v["NDIT"] FT_Start=self.ftstart hdul=fits.open(self.file,memmap=True,mode="readonly") data=hdul['OI_FLUX',self.ext].data for k in oi_flux_data.keys(): data_column=data.field(k).reshape((N,4,-1)) if k == "TIME": oi_flux_data[k]=data_column[:,0,0]*1e-6+FT_Start elif data_column.shape[2] == 1: oi_flux_data[k]=data_column[:,:,0] else: oi_flux_data[k]=data_column hdul.close() self.oi_flux=oi_flux_data self.flux=self.oi_flux["FLUX"]/self.h_v["DIT"] def read_oi_vis(self,remove_telluric=False): if remove_telluric: tellurics=(self.wave<2.0799139e-06)|(self.wave>2.4096553e-06) # tellurics=(self.wave<2.0799139e-06)|(self.wave>2.3796553e-06) self.wave=self.wave[~tellurics] else: tellurics=self.wave>200 oi_vis_data={ "TIME":[], 'VISDATA_FT':[], 'VISDATA':[], 'VISERR':[], 'PHASE_REF':[], 'PHASE_REF_COEFF':[], 'FLAG':[], 'REJECTION_FLAG':[], 'UCOORD':[], 'VCOORD':[], 'MJD':[], 'OPD_DISP':[], 'PHASE_MET_TELFC':[], 'SELF_REF':[], 'SELF_REF_COEFF':[], } N=self.h_v["NDIT"] FT_Start=self.ftstart hdul=fits.open(self.file,memmap=True,mode="readonly") data=hdul['OI_VIS',self.ext].data for k in oi_vis_data.keys(): data_column=data.field(k).reshape((N,6,-1)) if k == "TIME": oi_vis_data[k]=data_column[:,0,0]*1e-6+FT_Start elif data_column.shape[2] == 1: oi_vis_data[k]=data_column[:,:,0] elif data_column.shape[2] < 8 : oi_vis_data[k]=data_column[:,:,:] else: oi_vis_data[k]=data_column[:,:,~tellurics] hdul.close() # remove metrology glitches for e in range(len(oi_vis_data["FLAG"])): oi_vis_data["FLAG"][e]|=np.isnan(oi_vis_data["PHASE_MET_TELFC"][e]) oi_vis_data["PHASE_MET_TELFC"][e][np.isnan(oi_vis_data["PHASE_MET_TELFC"][e])]=0 oi_vis_data['phaseMet'] = 2*pi/self.wave*oi_vis_data["OPD_DISP"] + oi_vis_data["PHASE_MET_TELFC"] oi_vis_data['arcLength']=calculatePhaseRefArclength(oi_vis_data["PHASE_REF_COEFF"], self.wave) self.flag = oi_vis_data["FLAG"] | ((oi_vis_data["REJECTION_FLAG"]&19) >0)[:,:,None] self.flag |= self.h_v["flag"] self.flag |= oi_vis_data['arcLength'][:,:,None]>5 # 5 self.flag |= oi_vis_data["FLAG"].mean(axis=2)[:,:,None]>.3 if np.any(np.array(Bad_Tel_1)==self.date): print("Bad files on telescope 1:",self.date) self.flag[:,0]=True self.flag[:,1]=True self.flag[:,2]=True if np.any(np.array(Bad_Tel_2)==self.date): print("Bad files on telescope 2:",self.date) self.flag[:,0]=True self.flag[:,3]=True self.flag[:,4]=True if np.any(np.array(Bad_Tel_3)==self.date): print("Bad files on telescope 3:",self.date) self.flag[:,1]=True self.flag[:,3]=True self.flag[:,5]=True if np.any(np.array(Bad_Tel_4)==self.date): print("Bad files on telescope 4:",self.date) self.flag[:,2]=True self.flag[:,4]=True self.flag[:,5]=True self.oi_vis=oi_vis_data self.ucoord =oi_vis_data["UCOORD"][:,:,None]*2*pi/self.wave/1000/(180/pi*3600) self.vcoord =oi_vis_data["VCOORD"][:,:,None]*2*pi/self.wave/1000/(180/pi*3600) self.uvmax = np.max(np.sqrt( self.ucoord**2 + self.vcoord**2 )) self.ucoord_diff = np.diff(self.ucoord,axis=0).mean(axis=0) self.vcoord_diff = np.diff(self.vcoord,axis=0).mean(axis=0) self.visData=self.oi_vis["VISDATA"]/self.h_v["DIT"] self.visData[self.flag]=0 self.visData*=exp(1j*(self.oi_vis["PHASE_REF"]-self.oi_vis["phaseMet"])) self.visData_ft=self.oi_vis["VISDATA_FT"]/self.h_v["DIT"] self.var=abs(self.oi_vis["VISERR"]/self.h_v["DIT"])**2/2 self.weight=1/self.var try: self.weight[self.flag]=1e-3*np.min(self.weight[~self.flag]) except: self.weight[self.flag]=1e-15 self.visData_weighted=self.visData*self.weight self.visData_ftnormed=self.visData/abs(self.visData_ft).mean(axis=2)[:,:,None] self.weight_ftnormed=self.weight*abs(self.visData_ft).mean(axis=2)[:,:,None]**2 self.visData_ftnormed_weighted=self.visData_ftnormed*self.weight_ftnormed def compact_data(self,full=False): # self.visData_raw=self.visData # self.visData_ftnormed_raw=self.visData_ftnormed if full: fact=self.ndit_compact else: fn=factorize_number(self.ndit_compact) fact=1 for fi in fn: if self.ndit_compact/fact>11: fact*=fi self.ndit_compact/=fact phase=self.ucoord*self.sxy[0] + self.vcoord*self.sxy[1] def vec_col(vec,fact=fact,mean=False): if mean: return vec.reshape(fact,-1,*vec.shape[1:]).mean(axis=0) else: return vec.reshape(fact,-1,*vec.shape[1:]).sum(axis=0) self.visData=vec_col(self.visData_weighted*exp(1j*phase))/vec_col(self.weight) self.visData_ftnormed=vec_col(self.visData_ftnormed_weighted*exp(1j*phase))/vec_col(self.weight_ftnormed) self.ucoord_diff*=fact self.vcoord_diff*=fact self.ucoord=vec_col(self.ucoord,mean=True) self.vcoord=vec_col(self.vcoord,mean=True) self.weight=vec_col(self.weight) self.weight_ftnormed=vec_col(self.weight_ftnormed) phase=self.ucoord*self.sxy[0] + self.vcoord*self.sxy[1] self.visData=self.visData*exp(-1j*phase) self.visData_ftnormed=self.visData_ftnormed*exp(-1j*phase) self.visData_weighted=self.visData*self.weight self.visData_ftnormed_weighted=self.visData_ftnormed*self.weight_ftnormed # print(self.ndit_compact*fact,self.ndit_compact) class bTextColors: HEADER = '\033[95m' OKBLUE = '\033[94m' OKGREEN = '\033[92m' WARNING = '\033[93m' FAIL = '\033[91m' ENDC = '\033[0m' BOLD = '\033[1m' UNDERLINE = '\033[4m' # Code Mathias def calculatePhaseRefArclength(phaseref_coeffs, wav): """ Calculate the arclength for the polynomial fit to phaseRef args: phaseref_coeffs: array of size (ndit, nchannel, 4) containing the coeffs of the poly fit to phase ref wav: wavelength axis (in m) """ # calculate normalized wavelength grid used to fit phaseRef w0 = wav.mean() wref = (w0/wav - 1)/(np.max(wav) - np.min(wav))*w0 # prepare array for values ndit = np.shape(phaseref_coeffs)[0] nchannel = np.shape(phaseref_coeffs)[1] phaseRefArclength = np.zeros([ndit, nchannel]) # calculate arclength for dit in range(ndit): for c in range(nchannel): coeffs = phaseref_coeffs[dit, c, :] # calculate the coeffs of 1+(dy/dx)**2 arclengthpolycoeff = np.array([1+coeffs[1]**2, 4*coeffs[1]*coeffs[2], 4*coeffs[2]**2+6*coeffs[1]*coeffs[3], 12*coeffs[2]*coeffs[3], 9*coeffs[3]**2]) # integrate sqrt(1+(dx/dy)**2) do get arclength. Minus sign because wref is decreasing phaseRefArclength[dit, c] = -np.trapz(np.sqrt(np.polyval(arclengthpolycoeff[::-1], wref)), wref) return phaseRefArclength def angle_unwrap(C): Cmean=C.sum(axis=0) P=angle(C*conj(Cmean)) return P+angle(Cmean) #%% OI_header def get_header(files,verbose=False): h_v={ "file":[], "file_i":[], "extension":[], "DATE":[], "RA":[], "DEC":[], "DIT":[], "DITFT":[], "NDIT":[], "RES":[], "AXIS":[], "POLA":[], "SXY":[], "SXY_CTU":[], "SXY_KAL":[], "AST_KAL":[], "AST_MOD":[], "LKDT":[], "FT_START":[], "WAVELENG":[], "TAU0_START":[], "TAU0_END":[], "SEEING_START":[], "SEEING_END":[], "AIRM_START":[], "AIRM_END":[], "DATE_END":[], "NAME":[], "NAME_FT":[], "NAME_SC":[], "MET_ROW":[], "MET_DIFF":[], "REJ_RATIO":[], "MET_ROW":[], "UT":[], "ADDMET":[], "SEPTRK":[], "STATIONS":[], "PROG_ID":[], "flag":[], } for file_i,f in enumerate(files): fheader=fits.getheader(f) if fheader['HIERARCH ESO INS POLA MODE'] == "COMBINED": extension_list=[10] else: extension_list=[11,12] for extension in extension_list: h_v["file"]+=[f] h_v["file_i"]+=[file_i] h_v["extension"]+=[extension] h_v["DATE"]+=[fheader["DATE-OBS"]] h_v["DIT"]+=[fheader["HIERARCH ESO DET2 SEQ1 DIT"]] h_v["DITFT"]+=[fheader["HIERARCH ESO DET3 SEQ1 DIT"]] h_v["NDIT"]+=[fheader["HIERARCH ESO DET2 NDIT"]] h_v["RES"]+=[fheader["HIERARCH ESO INS SPEC RES"]] h_v["AXIS"]+=[fheader["HIERARCH ESO INS OPTI11 ID"]] h_v["POLA"]+=[fheader["HIERARCH ESO INS POLA MODE"]] h_v["PROG_ID"]+=[fheader["HIERARCH ESO OBS PROG ID"]] try: h_v["SXY"]+=[array([fheader["HIERARCH ESO INS SOBJ X"],fheader["HIERARCH ESO INS SOBJ Y"]])] except: h_v["SXY"]+=[array([0.0,0.0])] try: h_v["SXY_CTU"]+=[array([(fheader["HIERARCH ESO FT KAL SEPCTUX%i"%i],fheader["HIERARCH ESO FT KAL SEPCTUY%i"%i]) for i in 1+arange(4)])] h_v["SXY_KAL"]+=[array([(fheader["HIERARCH ESO FT KAL SEPX%i"%i],fheader["HIERARCH ESO FT KAL SEPY%i"%i]) for i in 1+arange(4)])] h_v["AST_KAL"]+=[array([fheader["HIERARCH ESO FT KAL AST"+str(i+1)] for i in range(4)])] h_v["AST_MOD"]+=[array([fheader["HIERARCH ESO FT KAL ASTMOD"+str(i+1)] for i in range(4)])] except: h_v["SXY_CTU"]+=[array([h_v["SXY"][-1] for i in 1+arange(4)])] h_v["SXY_KAL"]+=[array([h_v["SXY"][-1] for i in 1+arange(4)])] h_v["AST_KAL"]+=[array([0.0,0.0])] h_v["AST_MOD"]+=[array([0.0,0.0])] print("OLD DATA???? !!!! No KALCTU...") h_v["LKDT"]+=[array([fheader["HIERARCH ESO OCS MET LKDT_FC"+str(i+1)] for i in range(4)])] h_v["FT_START"]+=[astroTime([fheader["HIERARCH ESO PCR ACQ START"]],format='isot',scale='utc').mjd[0]*24*3600] h_v["WAVELENG"]+=[fheader["HIERARCH ESO QC MET LAMBDA MEAN"]] try: h_v["NAME"]+=[fheader["HIERARCH ESO OBS TARG NAME"]] h_v["NAME_SC"]+=[fheader["HIERARCH ESO INS SOBJ NAME"]] h_v["NAME_FT"]+=[fheader["HIERARCH ESO FT ROBJ NAME"]] h_v["RA"]+=[fheader["RA"]] h_v["DEC"]+=[fheader["DEC"]] except: h_v["NAME"]+=["Noname"] h_v["NAME_SC"]+=["Noname_SC"] h_v["NAME_FT"]+=["Noname_FT"] h_v["RA"]+=[0.0] h_v["DEC"]+=[0.0] h_v["MET_ROW"]+=[fheader["HIERARCH ESO QC MET REF ROW"]] h_v["MET_DIFF"]+=[array([fheader["HIERARCH ESO QC MET DRS DIFF"+str(i+1)] for i in range(4)])] h_v["REJ_RATIO"]+=[array([fheader["HIERARCH ESO QC REJECTED RATIO FT"],fheader["HIERARCH ESO QC REJECTED RATIO SC"]])] try: h_v["UT"]+=[fheader["HIERARCH ESO ISS CONF T1NAME"]] h_v["STATIONS"]+=[fheader["HIERARCH ESO ISS CONF STATION1"]+ fheader["HIERARCH ESO ISS CONF STATION2"]+ fheader["HIERARCH ESO ISS CONF STATION3"]+ fheader["HIERARCH ESO ISS CONF STATION4"]] h_v["AIRM_START"]+=[fheader["HIERARCH ESO ISS AIRM START"]] h_v["SEEING_START"]+=[fheader["HIERARCH ESO ISS AMBI FWHM START"]] h_v["TAU0_START"]+=[fheader["HIERARCH ESO ISS AMBI TAU0 START"]] h_v["AIRM_END"]+=[fheader["HIERARCH ESO ISS AIRM END"]] h_v["SEEING_END"]+=[fheader["HIERARCH ESO ISS AMBI FWHM END"]] h_v["TAU0_END"]+=[fheader["HIERARCH ESO ISS AMBI TAU0 END"]] h_v["DATE_END"]+=[fheader["HIERARCH ESO PCR ACQ END"]] except: h_v["UT"]+=["None"] h_v["STATIONS"]+=["None"] h_v["AIRM_START"]+=[0] h_v["SEEING_START"]+=[0] h_v["TAU0_START"]+=[0] h_v["AIRM_END"]+=[0] h_v["SEEING_END"]+=[0] h_v["TAU0_END"]+=[0] h_v["DATE_END"]+=["None"] try: h_v["SEPTRK"]+=[fheader["HIERARCH ESO FT KAL SEPTRK"]] except: h_v["SEPTRK"]+=[0] h_v["ADDMET"]+=[np.ones(4)] for tel in range(4): try: off=fheader["HIERARCH ESO ADD MET PH_FC%i_FT"%(tel+1)] h_v["ADDMET"][-1][tel]+=off*2*pi except: pass try: off=fheader["HIERARCH ESO ADD MET PH_FC%i_SC"%(tel+1)] h_v["ADDMET"][-1][tel]-=off*2*pi except: pass if np.any(np.array(badfilesList)==h_v["DATE"][-1]): h_v["flag"]+=[True] else: h_v["flag"]+=[False] h_v["sxy_best"]=[] for sxy,sxy_ctu,septrk in zip(h_v["SXY"],h_v["SXY_CTU"],h_v["SEPTRK"]): if septrk>0.5: h_v["sxy_best"]+=[sxy] else: h_v["sxy_best"]+=[sxy_ctu.mean(axis=0)] for k in h_v.keys(): h_v[k]=np.array(h_v[k]) if len(h_v["SXY"])>0: h_v['sep']=sqrt(h_v["SXY"][:,0]**2+h_v["SXY"][:,1]**2) h_v['comment']=[] h_v['comment_color']=[] for data,dit,ndit,red,pola,axis,sxy,name in zip(h_v["DATE"], h_v["DIT"], h_v["NDIT"], h_v["RES"], h_v["POLA"], h_v["AXIS"], h_v["SXY"], h_v["NAME"]): h_v['comment']+=[data+" "+str(dit)+" "+str(ndit)+" "+red+" "+pola+" "+axis+" "+str(sxy)+" "+name] for comment,flag in zip(h_v["comment"],h_v["flag"]): h_v['comment_color']+=[comment] if flag: h_v['comment_color'][-1]=bTextColors.FAIL+h_v['comment_color'][-1]+bTextColors.ENDC h_v['extensionANDdate']=[] for ext,date in zip(h_v["extension"],h_v["FT_START"]): h_v['extensionANDdate']+=[ext+100*int(date/24/3600+0.5)] h_v['verification']=[] h_v['verification_color']=[] for data,rej, met in zip(h_v["DATE"],h_v["REJ_RATIO"],h_v["MET_DIFF"]): string_rej="reject FT=%0.2i%%,SC=%0.2i%%"%(rej[0],rej[1]) string_rej_color=string_rej if rej[1]>50: string_rej_color=bTextColors.FAIL+string_rej+bTextColors.ENDC elif rej[0]>20: string_rej_color=bTextColors.WARNING+string_rej+bTextColors.ENDC string_met=" DRS/TAC-MET=%.3f,%.3f,%.3f,%.3f"%(met[0],met[1],met[2],met[3]) string_met_color=string_met if max(met)>0.1: string_met_color=bTextColors.FAIL+string_met+bTextColors.ENDC elif max(met)>0.01: string_met_color=bTextColors.WARNING+string_met+bTextColors.ENDC h_v['verification']+=[data+" "+string_rej+string_met] h_v['verification_color']+=[data+" "+string_rej_color+string_met_color] if verbose==True: for i,e in enumerate(h_v["extension"]): if (e==11)|(e==10): print(h_v["comment_color"][i]) for i,e in enumerate(h_v["extension"]): if (e==11)|(e==10): print(h_v["verification_color"][i]) return h_v #%% OI_VIS DATA def get_oi_vis(h_v,verbose=False,files_dir=None): oi_vis_data={ "TIME":[], 'VISDATA_FT':[], 'VISDATA':[], 'VISERR':[], 'PHASE_REF':[], 'PHASE_REF_COEFF':[], 'FLAG':[], 'REJECTION_FLAG':[], 'UCOORD':[], 'VCOORD':[], 'MJD':[], 'OPD_DISP':[], 'PHASE_MET_TELFC':[], 'SELF_REF':[], 'SELF_REF_COEFF':[], } for f,extension,N,FT_Start in zip(h_v["file"],h_v["extension"],h_v["NDIT"],h_v["FT_START"]): if files_dir: f=files_dir+f if verbose==True: print(f) hdul=fits.open(f,memmap=True,mode="readonly") data=hdul['OI_VIS',extension].data for k in oi_vis_data.keys(): data_column=data.field(k).reshape((N,6,-1)) if k == "TIME": oi_vis_data[k]+=[data_column[:,0,0]*1e-6+FT_Start] elif data_column.shape[2] == 1: oi_vis_data[k]+=[data_column[:,:,0]] else: oi_vis_data[k]+=[data_column] hdul.close() # remove metrology glitches for e in range(len(oi_vis_data["FLAG"])): oi_vis_data["FLAG"][e]|=np.isnan(oi_vis_data["PHASE_MET_TELFC"][e]) oi_vis_data["PHASE_MET_TELFC"][e][np.isnan(oi_vis_data["PHASE_MET_TELFC"][e])]=0 oi_vis_data['wave']=[] oi_vis_data['phaseMet']=[] oi_vis_data['ucoord']=[] oi_vis_data['vcoord']=[] oi_vis_data['arcLength']=[] oi_vis_data['flag']=[] for f,extension in zip(h_v["file"],h_v["extension"]): if files_dir: f=files_dir+f oi_vis_data['wave']+=[fits.getdata(f,'OI_WAVELENGTH',extension).field('EFF_WAVE')] for OPD_DISP,PHASE_MET_TELFC,wave in zip(oi_vis_data["OPD_DISP"],oi_vis_data["PHASE_MET_TELFC"],oi_vis_data["wave"]): oi_vis_data['phaseMet'] += [2*pi/wave*OPD_DISP + PHASE_MET_TELFC] for UCOORD,VCOORD,wave in zip(oi_vis_data["UCOORD"],oi_vis_data["VCOORD"],oi_vis_data["wave"]): oi_vis_data['ucoord'] +=[UCOORD[:,:,None]*2*pi/wave] oi_vis_data['vcoord'] +=[VCOORD[:,:,None]*2*pi/wave] for phaseref_coeffs, wave in zip(oi_vis_data["PHASE_REF_COEFF"],oi_vis_data["wave"]): oi_vis_data['arcLength']+=[calculatePhaseRefArclength(phaseref_coeffs, wave)] for flag,flag2,arcLength,flagBadFile in zip(oi_vis_data["FLAG"],oi_vis_data["REJECTION_FLAG"],oi_vis_data["arcLength"],h_v["flag"]): flag3 = flag | ((flag2&19) >0)[:,:,None] flag3|=flagBadFile flag3|=arcLength[:,:,None]>5 # 5 flag3|=flag.mean(axis=2)[:,:,None]>.3 oi_vis_data['flag'] +=[flag3] return oi_vis_data #%% OI_FLUX DATA def get_oi_flux(h_v): oi_flux_data={ "TIME":[], 'FLUX':[], 'FDDL_FT':[], 'FDDL_SC':[], 'PHASE_MET_FC':[], 'PHASE_MET_FCFT':[], 'PHASE_MET_FCSC':[], } for f,extension,N,FT_Start in zip(h_v["file"],h_v["extension"],h_v["NDIT"],h_v["FT_START"]): hdul=fits.open(f,memmap=True,mode="readonly") data=hdul['OI_FLUX',extension].data for k in oi_flux_data.keys(): data_column=data.field(k).reshape((N,4,-1)) if k == "TIME": oi_flux_data[k]+=[data_column[:,0,0]*1e-6+FT_Start] elif data_column.shape[2] == 1: oi_flux_data[k]+=[data_column[:,:,0]] else: oi_flux_data[k]+=[data_column] hdul.close() return oi_flux_data #%% OI_VIS_MET DATA def get_oi_met(h_v): oi_met_data={ "TIME":[], 'PHASE_FC_DRS':[], 'PHASE_TELFC_CORR':[], 'OPD_TELFC_CORR_XY':[], 'OPD_FC_CORR':[], 'OPD_TELFC_MCORR':[] } for f,FT_Start,extension in zip(h_v["file"],h_v["FT_START"],h_v["extension"]): if extension < 12: hdul=fits.open(f,memmap=True,mode="readonly") data=hdul['OI_VIS_MET'].data N=int(len(data)/4) for k in oi_met_data.keys(): data_column=data.field(k).reshape((N,4,-1)) if k == "TIME": oi_met_data[k]+=[data_column[:,0,0]*1e-6+FT_Start] elif data_column.shape[2] == 1: oi_met_data[k]+=[data_column[:,:,0]] else: oi_met_data[k]+=[data_column] hdul.close() oi_met_data['astRes']=[] oi_met_data['sepXRes']=[] oi_met_data['sepYRes']=[] for OPD_TELFC_CORR_XY,wave in zip(oi_met_data["OPD_TELFC_CORR_XY"],h_v["WAVELENG"]): PHASE_TELFC_CORR_XY=OPD_TELFC_CORR_XY/wave*2*pi oi_met_data['astRes']+=[PHASE_TELFC_CORR_XY[:,:,3]-PHASE_TELFC_CORR_XY[:,:,2]+PHASE_TELFC_CORR_XY[:,:,1]-PHASE_TELFC_CORR_XY[:,:,0]] oi_met_data['sepXRes']+=[PHASE_TELFC_CORR_XY[:,:,3]-PHASE_TELFC_CORR_XY[:,:,2]-PHASE_TELFC_CORR_XY[:,:,1]+PHASE_TELFC_CORR_XY[:,:,0]] oi_met_data['sepYRes']+=[PHASE_TELFC_CORR_XY[:,:,3]+PHASE_TELFC_CORR_XY[:,:,2]-PHASE_TELFC_CORR_XY[:,:,1]-PHASE_TELFC_CORR_XY[:,:,0]] return oi_met_data #%% OI_ACQ def get_oi_acq(h_v): oi_acq_data={ "TIME":[], 'PUPIL_NSPOT':[], 'PUPIL_X':[], 'PUPIL_Y':[], "PUPIL_ACQ":[], "PUPIL_FIELD":[] } for f,FT_Start,extension in zip(h_v["file"],h_v["FT_START"],h_v["extension"]): if extension < 12: hdul=fits.open(f,memmap=True,mode="readonly") try: data=hdul['OI_VIS_ACQ'].data N=int(len(data)/4) X=data.field("PUPIL_X").reshape((N,4)) Y=data.field("PUPIL_Y").reshape((N,4)) G=(X**2+Y**2).sum(axis=1)>0.001 for k in oi_acq_data.keys(): data_column=data.field(k).reshape((N,4,-1))[G] if k == "TIME": oi_acq_data[k]+=[data_column[:,0,0]*1e-6+FT_Start] elif data_column.shape[2] == 1: oi_acq_data[k]+=[data_column[:,:,0]] else: oi_acq_data[k]+=[data_column] except: pass try: data=hdul['IMAGING_DATA_ACQ'].data[0] Ns=55 oi_acq_data["PUPIL_FIELD"]+=[np.array([data[120-Ns:120+Ns,120+250*i-Ns:120+250*i+Ns] for i in range(4)])] oi_acq_data["PUPIL_ACQ"]+=[np.array([data[870-Ns:870+Ns,120+250*i-Ns:120+250*i+Ns] for i in range(4)])] except: pass hdul.close() return oi_acq_data #%%