# -*- coding: utf-8 -*- """ Created on Mon Sep 14 11:41:26 2015 @author: kervella """ try: import pyfits except: from astropy.io import fits as pyfits import numpy as np #import os #import gc # =============================================================================== # Python classes to read GRAVITY files at different stages of the data processing # =============================================================================== # # Oifits, Oilist, P2vm, Wave, Preproc, P2vmreduced, Rawdata, Badpix, Profile, Dark, Argon # class Oifits: def __init__(self, filename): # open fits file self.filename = filename gravi_oifits=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_oifits[0].header.copy() self.exptime_sc = self.header['HIERARCH ESO DET2 SEQ1 DIT'] # Mode of the instrument if 'SINGLE' in str(self.header['HIERARCH ESO PRO CATG']): self.single = True else: self.single = False if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True # get the FDDL, OI_ARRAY, OI_WAVELENGTH and OI_VIS_MET extentions for hdu in gravi_oifits : # OI_TARGET if hdu.name == 'OI_TARGET': # test if OI_TARGET is present in header self.oi_target = gravi_oifits['OI_TARGET'].data if hdu.name == 'OI_ARRAY' : oi_array=hdu.data.copy() if hdu.name == 'ARRAY_GEOMETRY' : # TBC, not present (2015-09) oi_array=hdu.data.copy() # get the OI_VIS, OI_FLUX and T3 extensions in combined polarization mode if self.polarsplit == False: for hdu in gravi_oifits : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_vis_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') : oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_vis2_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS2') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') : oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_flux_sc=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc.field('TIME')[::4]/3600. # in decimal hours headerflux = hdu.header.copy() columnflux = hdu.columns.names # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') : oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_T3') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_t3_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_t3_ft=hdu.data.copy() if (hdu.name == 'OI_T3_FT') : oi_t3_ft=hdu.data.copy() # get the OI_VIS, OI_FLUX and T3 extensions in split polarization mode (S considered =P1, P considered =P2) if self.polarsplit == True: for hdu in gravi_oifits : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_vis_sc_p=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_vis_sc_s=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_vis2_sc_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_vis2_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_flux_sc_p=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc_p.field('TIME')[::4]/3600. # in decimal hours headerflux = hdu.header.copy() columnflux = hdu.columns.names if (hdu.name == 'OI_FLUX') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_flux_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_FLUX_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_t3_sc_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_t3_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() gravi_oifits.close() del gravi_oifits # gc.collect() # Wavelength scales self.wave_sc=oi_wave_sc.field('EFF_WAVE')*1.e6 # wavelength scale in microns self.minwave_sc=np.min(self.wave_sc) self.maxwave_sc=np.max(self.wave_sc) self.wave_ft=oi_wave_ft.field('EFF_WAVE')*1.e6 self.minwave_ft=np.min(self.wave_sc) self.maxwave_ft=np.max(self.wave_sc) self.nwave_sc = self.wave_sc.shape[0] self.nwave_ft = self.wave_ft.shape[0] # OI_ARRAY parameters self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') # Time scale of the different records in the OIFITS file self.time = self.oi_flux_sc_time # Number of records self.nrecords = self.time.shape[0] # OI_FLUX table parameters if self.polarsplit == False: if 'VISFLUX' in columnflux: self.oi_flux_sc = oi_flux_sc.field('VISFLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('VISFLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc = oi_flux_sc.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft = oi_flux_ft.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_ft) else: self.oi_flux_sc = oi_flux_sc.field('FLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('FLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc = oi_flux_sc.field('FLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft = oi_flux_ft.field('FLUXERR').reshape(self.nrecords,4,self.nwave_ft) if self.polarsplit == True: if 'VISFLUX' in columnflux: self.oi_flux_sc_s = oi_flux_sc_s.field('VISFLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('VISFLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('VISFLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('VISFLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_s = oi_flux_sc_s.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_s = oi_flux_ft_s.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_p = oi_flux_sc_p.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_p = oi_flux_ft_p.field('VISFLUXERR').reshape(self.nrecords,4,self.nwave_ft) else: self.oi_flux_sc_s = oi_flux_sc_s.field('FLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('FLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('FLUX').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('FLUX').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_s = oi_flux_sc_s.field('FLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_s = oi_flux_ft_s.field('FLUXERR').reshape(self.nrecords,4,self.nwave_ft) self.oi_flux_err_sc_p = oi_flux_sc_p.field('FLUXERR').reshape(self.nrecords,4,self.nwave_sc) self.oi_flux_err_ft_p = oi_flux_ft_p.field('FLUXERR').reshape(self.nrecords,4,self.nwave_ft) # OI_VIS and T3 table parameters if self.polarsplit == False: self.oi_vis_sc_time = oi_vis_sc.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_sc_mjd = oi_vis_sc.field('MJD')[::6] self.oi_vis_sc_target_id = oi_vis_sc.field('TARGET_ID')[::6] self.oi_vis_sc_int_time = oi_vis_sc.field('INT_TIME') # in seconds # Complex visibilities self.oi_vis_sc_visdata = oi_vis_sc.field('VISDATA') self.oi_vis_sc_viserr = oi_vis_sc.field('VISERR') self.oi_vis_sc_visamp = oi_vis_sc.field('VISAMP') self.oi_vis_sc_visamperr = oi_vis_sc.field('VISAMPERR') self.oi_vis_sc_visphi = oi_vis_sc.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr = oi_vis_sc.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_flag = oi_vis_sc.field('FLAG') self.oi_vis_ft_time = oi_vis_ft.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_ft_mjd = oi_vis_ft.field('MJD')[::6] self.oi_vis_ft_target_id = oi_vis_ft.field('TARGET_ID')[::6] self.oi_vis_ft_int_time = oi_vis_ft.field('INT_TIME') # in seconds self.oi_vis_ft_visdata = oi_vis_ft.field('VISDATA') self.oi_vis_ft_viserr = oi_vis_ft.field('VISERR') self.oi_vis_ft_visamp = oi_vis_ft.field('VISAMP') self.oi_vis_ft_visamperr = oi_vis_ft.field('VISAMPERR') self.oi_vis_ft_visphi = oi_vis_ft.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr = oi_vis_ft.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) # u,v coordinates different for SC and FT nbaseline= 6 self.ucoord_ft = oi_vis_ft.field('UCOORD') self.vcoord_ft = oi_vis_ft.field('VCOORD') self.projbase_ft = np.zeros(nbaseline) self.azimbase_ft = np.zeros(nbaseline) for baseline in range(0,nbaseline): self.projbase_ft[baseline] = np.sqrt(np.square(self.ucoord_ft[baseline])+np.square(self.vcoord_ft[baseline])) self.azimbase_ft[baseline] = np.arctan2(self.vcoord_ft[baseline], self.ucoord_ft[baseline]) * 180./np.pi self.ucoord_sc = oi_vis_sc.field('UCOORD') self.vcoord_sc = oi_vis_sc.field('VCOORD') self.projbase_sc = np.zeros(nbaseline) self.azimbase_sc = np.zeros(nbaseline) for baseline in range(0,nbaseline): self.projbase_sc[baseline] = np.sqrt(np.square(self.ucoord_sc[baseline])+np.square(self.vcoord_sc[baseline])) self.azimbase_sc[baseline] = np.arctan2(self.vcoord_sc[baseline], self.ucoord_sc[baseline]) * 180./np.pi # Station index common self.sta_index = oi_vis_sc.field('STA_INDEX') # Squared visibilities self.oi_vis2_sc_vis2data = oi_vis2_sc.field('VIS2DATA') self.oi_vis2_sc_vis2err = oi_vis2_sc.field('VIS2ERR') self.oi_vis2_sc_ucoord = oi_vis2_sc.field('UCOORD') self.oi_vis2_sc_vcoord = oi_vis2_sc.field('VCOORD') self.oi_vis2_sc_sta_index = oi_vis2_sc.field('STA_INDEX') self.oi_vis2_sc_flag = oi_vis2_sc.field('FLAG') self.oi_vis2_ft_vis2data = oi_vis2_ft.field('VIS2DATA') self.oi_vis2_ft_vis2err = oi_vis2_ft.field('VIS2ERR') self.oi_vis2_ft_ucoord = oi_vis2_ft.field('UCOORD') self.oi_vis2_ft_vcoord = oi_vis2_ft.field('VCOORD') self.oi_vis2_ft_sta_index = oi_vis2_ft.field('STA_INDEX') self.oi_vis2_ft_flag = oi_vis2_ft.field('FLAG') # Closure quantities self.t3_u1coord = oi_t3_sc.field('U1COORD') self.t3_v1coord = oi_t3_sc.field('V1COORD') self.t3_baseline1 = np.sqrt(np.square(self.t3_u1coord)+np.square(self.t3_v1coord)) self.t3_u2coord = oi_t3_sc.field('U2COORD') self.t3_v2coord = oi_t3_sc.field('V2COORD') self.t3_baseline2 = np.sqrt(np.square(self.t3_u2coord)+np.square(self.t3_v2coord)) self.t3_u3coord = self.t3_u1coord+self.t3_u2coord self.t3_v3coord = self.t3_v1coord+self.t3_v2coord self.t3_baseline3 = np.sqrt(np.square(self.t3_u3coord)+np.square(self.t3_v3coord)) self.t3_baseavg = np.mean([self.t3_baseline1,self.t3_baseline2,self.t3_baseline3],axis=0) self.t3phi_sc = oi_t3_sc.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_sc = oi_t3_sc.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_ft = oi_t3_ft.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_ft = oi_t3_ft.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3amp_sc = oi_t3_sc.field('T3AMP') self.t3amperr_sc = oi_t3_sc.field('T3AMPERR') self.t3amp_ft = oi_t3_ft.field('T3AMP') self.t3amperr_ft = oi_t3_ft.field('T3AMPERR') self.t3_flag = oi_t3_ft.field('FLAG') # Station indices for the closure quantities self.t3_sc_staindex = oi_t3_sc.field('STA_INDEX') self.t3_ft_staindex = oi_t3_ft.field('STA_INDEX') # Astrometric quantities (in SC OIVIS table) self.gdelay_astrometry = oi_vis_sc.field('GDELAY_ASTROMETRY') self.phase_astrometry = oi_vis_sc.field('PHASE_ASTROMETRY') if self.polarsplit == True: self.oi_vis_sc_time = oi_vis_sc_s.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_sc_mjd = oi_vis_sc_s.field('MJD')[::6] self.oi_vis_sc_target_id = oi_vis_sc_s.field('TARGET_ID') # same for s and p polar self.oi_vis_sc_int_time = oi_vis_sc_s.field('INT_TIME') # in seconds # Complex visibilities self.oi_vis_sc_visdata_s = oi_vis_sc_s.field('VISDATA') self.oi_vis_sc_viserr_s = oi_vis_sc_s.field('VISERR') self.oi_vis_sc_visamp_s = oi_vis_sc_s.field('VISAMP') self.oi_vis_sc_visamperr_s = oi_vis_sc_s.field('VISAMPERR') self.oi_vis_sc_visphi_s = oi_vis_sc_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr_s = oi_vis_sc_s.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_flag_s = oi_vis_sc_s.field('FLAG') self.oi_vis_sc_visdata_p = oi_vis_sc_p.field('VISDATA') self.oi_vis_sc_viserr_p = oi_vis_sc_p.field('VISERR') self.oi_vis_sc_visamp_p = oi_vis_sc_p.field('VISAMP') self.oi_vis_sc_visamperr_p = oi_vis_sc_p.field('VISAMPERR') self.oi_vis_sc_visphi_p = oi_vis_sc_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr_p = oi_vis_sc_p.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_time = oi_vis_ft_s.field('TIME')[::6]/3600. # in decimal hours self.oi_vis_ft_mjd = oi_vis_ft_s.field('MJD')[::6] self.oi_vis_ft_target_id = oi_vis_ft_s.field('TARGET_ID') # same for s and p polar self.oi_vis_ft_int_time = oi_vis_ft_s.field('INT_TIME') # in seconds self.oi_vis_ft_visdata_s = oi_vis_ft_s.field('VISDATA') self.oi_vis_ft_viserr_s = oi_vis_ft_s.field('VISERR') self.oi_vis_ft_visamp_s = oi_vis_ft_s.field('VISAMP') self.oi_vis_ft_visamperr_s = oi_vis_ft_s.field('VISAMPERR') self.oi_vis_ft_visphi_s = oi_vis_ft_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr_s = oi_vis_ft_s.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_sta_index_s = oi_vis_ft_s.field('STA_INDEX') self.oi_vis_ft_visdata_p = oi_vis_ft_p.field('VISDATA') self.oi_vis_ft_viserr_p = oi_vis_ft_p.field('VISERR') self.oi_vis_ft_visamp_p = oi_vis_ft_p.field('VISAMP') self.oi_vis_ft_visamperr_p = oi_vis_ft_p.field('VISAMPERR') self.oi_vis_ft_visphi_p = oi_vis_ft_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr_p = oi_vis_ft_p.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) # u,v coordinates common to both polarizations, different for SC and FT nbaseline= 6 self.ucoord_ft = oi_vis_ft_p.field('UCOORD') self.vcoord_ft = oi_vis_ft_p.field('VCOORD') self.projbase_ft = np.zeros(nbaseline) self.azimbase_ft = np.zeros(nbaseline) for baseline in range(0,nbaseline): self.projbase_ft[baseline] = np.sqrt(np.square(self.ucoord_ft[baseline])+np.square(self.vcoord_ft[baseline])) self.azimbase_ft[baseline] = np.arctan2(self.vcoord_ft[baseline], self.ucoord_ft[baseline]) * 180./np.pi self.ucoord_sc = oi_vis_sc_p.field('UCOORD') self.vcoord_sc = oi_vis_sc_p.field('VCOORD') self.projbase_sc = np.zeros(nbaseline) self.azimbase_sc = np.zeros(nbaseline) for baseline in range(0,nbaseline): self.projbase_sc[baseline] = np.sqrt(np.square(self.ucoord_sc[baseline])+np.square(self.vcoord_sc[baseline])) self.azimbase_sc[baseline] = np.arctan2(self.vcoord_sc[baseline], self.ucoord_sc[baseline]) * 180./np.pi # Station index common to both polarizations, SC and FT self.sta_index = oi_vis_sc_s.field('STA_INDEX') # Squared visibilities self.oi_vis2_sc_vis2data_s = oi_vis2_sc_s.field('VIS2DATA') self.oi_vis2_sc_vis2err_s = oi_vis2_sc_s.field('VIS2ERR') self.oi_vis2_sc_sta_index_s = oi_vis2_sc_s.field('STA_INDEX') self.oi_vis2_sc_flag_s = oi_vis2_sc_s.field('FLAG') self.oi_vis2_sc_vis2data_p = oi_vis2_sc_p.field('VIS2DATA') self.oi_vis2_sc_vis2err_p = oi_vis2_sc_p.field('VIS2ERR') self.oi_vis2_sc_sta_index_p = oi_vis2_sc_p.field('STA_INDEX') self.oi_vis2_sc_flag_p = oi_vis2_sc_p.field('FLAG') self.oi_vis2_ft_vis2data_s = oi_vis2_ft_s.field('VIS2DATA') self.oi_vis2_ft_vis2err_s = oi_vis2_ft_s.field('VIS2ERR') self.oi_vis2_ft_flag_s = oi_vis2_ft_s.field('FLAG') self.oi_vis2_ft_vis2data_p = oi_vis2_ft_p.field('VIS2DATA') self.oi_vis2_ft_vis2err_p = oi_vis2_ft_p.field('VIS2ERR') self.oi_vis2_ft_flag_p = oi_vis2_ft_p.field('FLAG') # Closure quantities self.t3_u1coord = oi_t3_sc_s.field('U1COORD') self.t3_v1coord = oi_t3_sc_s.field('V1COORD') self.t3_baseline1 = np.sqrt(np.square(self.t3_u1coord)+np.square(self.t3_v1coord)) self.t3_u2coord = oi_t3_sc_s.field('U2COORD') self.t3_v2coord = oi_t3_sc_s.field('V2COORD') self.t3_baseline2 = np.sqrt(np.square(self.t3_u2coord)+np.square(self.t3_v2coord)) self.t3_u3coord = self.t3_u1coord+self.t3_u2coord self.t3_v3coord = self.t3_v1coord+self.t3_v2coord self.t3_baseline3 = np.sqrt(np.square(self.t3_u3coord)+np.square(self.t3_v3coord)) self.t3_baseavg = np.mean([self.t3_baseline1,self.t3_baseline2,self.t3_baseline3],axis=0) self.t3phi_sc_s = oi_t3_sc_s.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_sc_s = oi_t3_sc_s.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_ft_s = oi_t3_ft_s.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_ft_s = oi_t3_ft_s.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_sc_p = oi_t3_sc_p.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_sc_p = oi_t3_sc_p.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phi_ft_p = oi_t3_ft_p.field('T3PHI')*np.pi/180. # conversion in radians (after Sept. 16) self.t3phierr_ft_p = oi_t3_ft_p.field('T3PHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.t3amp_sc_s = oi_t3_sc_s.field('T3AMP') self.t3amperr_sc_s = oi_t3_sc_s.field('T3AMPERR') self.t3amp_ft_s = oi_t3_ft_s.field('T3AMP') self.t3amperr_ft_s = oi_t3_ft_s.field('T3AMPERR') self.t3_flag_s = oi_t3_ft_s.field('FLAG') self.t3amp_sc_p = oi_t3_sc_p.field('T3AMP') self.t3amperr_sc_p = oi_t3_sc_p.field('T3AMPERR') self.t3amp_ft_p = oi_t3_ft_p.field('T3AMP') self.t3amperr_ft_p = oi_t3_ft_p.field('T3AMPERR') self.t3_flag_p = oi_t3_ft_p.field('FLAG') # Station indices for closure quantities self.t3_sc_staindex = oi_t3_sc_s.field('STA_INDEX') self.t3_ft_staindex = oi_t3_ft_s.field('STA_INDEX') # Astrometric quantities (in SC OIVIS table) self.gdelay_astrometry_s = oi_vis_sc_s.field('GDELAY_ASTROMETRY') self.gdelay_astrometry_p = oi_vis_sc_p.field('GDELAY_ASTROMETRY') self.phase_astrometry_s = oi_vis_sc_s.field('PHASE_ASTROMETRY') self.phase_astrometry_p = oi_vis_sc_p.field('PHASE_ASTROMETRY') class Oilist: def __init__(self, filelist, fits_keywords): self.oi = [] filelist_out = [] timescale_out = [] oi_out = [] nfiles_in = len(filelist) nfiles = 0 for filename in filelist: gravi_file = pyfits.open(filename+'.fits',memmap=True) header = gravi_file[0].header.copy() gravi_file.close() del gravi_file # Final OIFITS data files only key_names = list(fits_keywords.keys()) type_ok = True for strname in key_names: type_ok *= (strname in header) if type_ok: keys_ok = True for name in key_names: keys_ok *= (header[name] in fits_keywords[name]) if keys_ok == True: nfiles = nfiles + 1 filelist_out.append(filename) print(' Loading OIFITS file %i/%i: %s.fits'%(nfiles,nfiles_in,filename)) a = Oifits(filename+'.fits') oi_out.append(a) timescale_out.append(float(a.header['MJD-OBS'])) # timescale of the measurements # TO BE CHANGED: (SEVERAL CAN BE PRESENT IN EACH FILE) del a del header # gc.collect() # Sort the output structure as a function of MJD of observation sorted_timescale = np.array(np.argsort(timescale_out)) # indices of sorted timescale self.filelist = [] self.timescale = [] for i in range(0,sorted_timescale.shape[0]): self.oi.append(oi_out[sorted_timescale[i]]) self.filelist.append(filelist_out[sorted_timescale[i]]) self.timescale.append(timescale_out[sorted_timescale[i]]) self.nfiles = nfiles class P2vm: def __init__(self, filename): # open fits file self.filename = filename gravi_p2vm=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_p2vm[0].header #print self.header['HIERARCH ESO PRO REC1 PIPE ID'] if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True #print self.polarsplit self.ports_sc = np.array(gravi_p2vm['IMAGING_DETECTOR_SC'].data['PORTS'])-1 # telescopes from 0 to 3 self.ports_ft = np.array(gravi_p2vm['IMAGING_DETECTOR_FT'].data['PORTS'])-1 self.regname_sc=gravi_p2vm['P2VM_SC'].data['REGNAME'] self.regname_ft=gravi_p2vm['P2VM_FT'].data['REGNAME'] self.transmission_sc=gravi_p2vm['P2VM_SC'].data['TRANSMISSION'] self.transmission_ft=gravi_p2vm['P2VM_FT'].data['TRANSMISSION'] self.transmission_met=gravi_p2vm['P2VM_MET'].data['TRANSMISSION'] self.nwave_sc=self.transmission_sc.shape[2] self.nwave_ft=self.transmission_ft.shape[2] self.nregion_sc=self.transmission_sc.shape[0] self.nregion_ft=self.transmission_ft.shape[0] self.coherence_sc=gravi_p2vm['P2VM_SC'].data['COHERENCE'] self.coherence_ft=gravi_p2vm['P2VM_FT'].data['COHERENCE'] self.coherence_met=gravi_p2vm['P2VM_MET'].data['COHERENCE'] self.phase_sc=gravi_p2vm['P2VM_SC'].data['PHASE'] self.phase_ft=gravi_p2vm['P2VM_FT'].data['PHASE'] self.phase_met=gravi_p2vm['P2VM_MET'].data['PHASE'] self.cmatrix_sc=np.array(gravi_p2vm['P2VM_SC'].data['C MATRIX']).reshape(self.nregion_sc,6,self.nwave_sc) self.cmatrix_ft=np.array(gravi_p2vm['P2VM_FT'].data['C MATRIX']).reshape(self.nregion_ft,6,self.nwave_ft) # get the OI_WAVELENGTH extention # these tables are present twice in the P2VM files (one for each polarization) but contain the same data if (self.polarsplit == False): for hdu in gravi_p2vm : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): self.wave_sc=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_SC') : self.wave_sc=hdu.data['EFF_WAVE']*1.e6 # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): self.wave_ft=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_FT') : self.wave_ft=hdu.data['EFF_WAVE']*1.e6 if (self.polarsplit == True): for hdu in gravi_p2vm : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): self.wave_sc=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_SC') : self.wave_sc=hdu.data['EFF_WAVE']*1.e6 # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): self.wave_ft=hdu.data['EFF_WAVE']*1.e6 if (hdu.name == 'OI_WAVELENGTH_FT') : self.wave_ft=hdu.data['EFF_WAVE']*1.e6 self.minwave_sc = np.min(self.wave_sc) self.maxwave_sc = np.max(self.wave_sc) self.minwave_ft = np.min(self.wave_ft) self.maxwave_ft = np.max(self.wave_ft) gravi_p2vm.close() class P2vmlist: def __init__(self, filelist, fits_keywords): self.p2vm = [] filelist_out = [] timescale_out = [] p2vm_out = [] nfiles_in = len(filelist) nfiles = 0 for filename in filelist: gravi_file = pyfits.open(filename+'.fits',memmap=True) header = gravi_file[0].header.copy() gravi_file.close() del gravi_file # P2VM data files only with proper keywords key_names = list(fits_keywords.keys()) type_ok = True for strname in key_names: type_ok *= (strname in header) if type_ok: keys_ok = True for name in key_names: keys_ok *= (header[name] in fits_keywords[name]) if keys_ok == True: nfiles = nfiles + 1 filelist_out.append(filename) print(' Loading P2VM file %i/%i: %s.fits'%(nfiles,nfiles_in,filename)) a = P2vm(filename+'.fits') p2vm_out.append(a) timescale_out.append(float(a.header['MJD-OBS'])) # timescale of the measurements # TO BE CHANGED: (SEVERAL CAN BE PRESENT IN EACH FILE) del a del header # gc.collect() # Sort the output structure as a function of MJD of observation sorted_timescale = np.array(np.argsort(timescale_out)) # indices of sorted timescale self.filelist = [] self.timescale = [] for i in range(0,sorted_timescale.shape[0]): self.p2vm.append(p2vm_out[sorted_timescale[i]]) self.filelist.append(filelist_out[sorted_timescale[i]]) self.timescale.append(timescale_out[sorted_timescale[i]]) self.nfiles = nfiles class Wave: def __init__(self, filename): # open fits file self.filename = filename gravi_wave=pyfits.open(filename, mode='readonly') # get the WAVE_DATA_SC extention self.valid = 0 for hdu in gravi_wave : #print hdu.name if hdu.name=='WAVE_DATA_SC' : self.valid=1 if self.valid == 0 : print("No WAVE_DATA_SC table in this file") return None self.valid = 0 for hdu in gravi_wave : #print hdu.name if hdu.name=='WAVE_DATA_FT' : self.valid=1 if self.valid == 0 : print("No WAVE_DATA_FT table in this file") return None # Main header self.header = gravi_wave[0].header if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True self.wavefiber_in = False for hdu in gravi_wave : if hdu.name=='WAVE_DATA_SC' : h = hdu.header self.startx = h['STARTX'] self.nwave_sc= h['NWAVE'] self.nregion_sc= h['TFIELDS'] if hdu.name=='WAVE_DATA_FT' : h = hdu.header self.nwave_ft= h['NWAVE'] self.nregion_ft= h['TFIELDS'] if hdu.name=='WAVE_FIBRE' : self.wavefiber_in = True if self.polarsplit == True: wfs = np.zeros((6,self.nwave_sc)) wfp = np.zeros((6,self.nwave_sc)) wfs[0,:]=hdu.data['BASE_12_S']*1E6 wfs[1,:]=hdu.data['BASE_13_S']*1E6 wfs[2,:]=hdu.data['BASE_14_S']*1E6 wfs[3,:]=hdu.data['BASE_23_S']*1E6 wfs[4,:]=hdu.data['BASE_24_S']*1E6 wfs[5,:]=hdu.data['BASE_34_S']*1E6 wfp[0,:]=hdu.data['BASE_12_P']*1E6 wfp[1,:]=hdu.data['BASE_13_P']*1E6 wfp[2,:]=hdu.data['BASE_14_P']*1E6 wfp[3,:]=hdu.data['BASE_23_P']*1E6 wfp[4,:]=hdu.data['BASE_24_P']*1E6 wfp[5,:]=hdu.data['BASE_34_P']*1E6 self.wavefiber_s = wfs self.wavefiber_p = wfp if self.polarsplit == False: wfc = np.zeros((6,self.nwave_sc)) wfc[0,:]=hdu.data['BASE_12_C']*1E6 wfc[1,:]=hdu.data['BASE_13_C']*1E6 wfc[2,:]=hdu.data['BASE_14_C']*1E6 wfc[3,:]=hdu.data['BASE_23_C']*1E6 wfc[4,:]=hdu.data['BASE_24_C']*1E6 wfc[5,:]=hdu.data['BASE_34_C']*1E6 self.wavefiber_c = wfc #self.wave_sc = wave_sc.reshape(self.nregion_sc,self.nwave_sc) #self.wave_ft = wave_ft.reshape(self.nregion_ft,self.nwave_ft) # Reshaping of wavelength vectors (output in microns) wavetable_sc = np.zeros([self.nregion_sc,self.nwave_sc],dtype='d') wave_sc = gravi_wave['WAVE_DATA_SC'].data for output in range(1,self.nregion_sc+1): wavetable_sc[output-1,:] = wave_sc.field('DATA'+str(output))[0].reshape(1,self.nwave_sc) self.wave_sc = wavetable_sc*1E6 wavetable_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') wave_ft = gravi_wave['WAVE_DATA_FT'].data for output in range(1,self.nregion_ft+1): wavetable_ft[output-1,:] = wave_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft) self.wave_ft = wavetable_ft*1E6 self.minwave_ft = np.min(self.wave_ft) self.maxwave_ft = np.max(self.wave_ft) imaging_det_sc = gravi_wave['IMAGING_DETECTOR_SC'].data self.regname_sc = imaging_det_sc.field('REGNAME') self.center_sc = (imaging_det_sc.field('CENTER')) imaging_det_ft = gravi_wave['IMAGING_DETECTOR_FT'].data self.regname_ft = imaging_det_ft.field('REGNAME') self.center_ft = (imaging_det_ft.field('CENTER')) if self.polarsplit == True: self.base_list_sc = [self.regname_sc[i*8][0:2] for i in range(0,6)] self.base_list_ft = [self.regname_ft[i*8][0:2] for i in range(0,6)] if self.polarsplit == False: self.base_list_sc = [self.regname_sc[i*4][0:2] for i in range(0,6)] self.base_list_ft = [self.regname_ft[i*4][0:2] for i in range(0,6)] # get the TEST_WAVE engineering image cube extention: map of the SC spectral extraction + map of wavelengths # 0: interpolated wavelength map # 1: profile map # 2: non-interpolated wavelength map (not available for older processed data) self.test_wave = gravi_wave['TEST_WAVE'].data self.test_wave[0,:,:] *= 1E6 # to get microns in the wavelength tables self.test_wave[2,:,:] *= 1E6 # to get microns in the wavelength tables ind=np.where(self.test_wave[0,:,:] > 1) self.minwave_sc = np.min(self.test_wave[0,ind[0],ind[1]]) self.maxwave_sc = np.max(self.test_wave[0,ind[0],ind[1]]) #self.minwave_sc = np.min(np.select(self.test_wave[0,:,:] > 1.0,self.test_wave[0,:,:])) #waself.maxwave_sc = np.max(np.select(self.test_wave[0,:,:] > 1.0,self.test_wave[0,:,:])) gravi_wave.close() class Preproc: def __init__(self, filename): # open fits file self.filename = filename gravi_preproc=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_preproc[0].header # Single mode of the instrument (no metrology data) #if 'SINGLE' in str(self.header['HIERARCH ESO DO CATG']): if 'SINGLE' in str(self.header['HIERARCH ESO PRO REC1 RAW1 CATG']): self.single = True else: self.single = False if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False self.nregion = 24 else: self.polarsplit = True self.nregion = 48 fddl = None # in case the FDDL data are not in the file detector_sc = None detector_ft = None oi_array = None spectrum_sc = None spectrum_ft = None metrology = None self.fddlpresent = False self.gdelay_in = False # get the FDDL, OI_ARRAY, OI_WAVELENGTH, OI_VIS_MET and OPDC extentions for hdu in gravi_preproc : if hdu.name == 'ARRAY_GEOMETRY' : # TBC oi_array=hdu.data if hdu.name == 'FDDL': fddl = hdu.data self.fddlpresent = True if hdu.name == 'IMAGING_DETECTOR_SC' : detector_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_FT' : detector_ft=hdu.data if hdu.name == 'METROLOGY' : metrology=hdu.data if hdu.name == 'SPECTRUM_DATA_FT' : spectrum_ft=hdu.data if hdu.name == 'SPECTRUM_DATA_SC' : spectrum_sc=hdu.data if hdu.name == 'OPDC' : opdc=hdu.data self.gdelay_in = True if self.polarsplit == False: for hdu in gravi_preproc : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if self.polarsplit == True: for hdu in gravi_preproc : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() # Wavelength scales self.wave_sc=oi_wave_sc.field('EFF_WAVE')*1.e6 # wavelength scale in microns self.minwave_sc=np.min(self.wave_sc) self.maxwave_sc=np.max(self.wave_sc) self.wave_ft=oi_wave_ft.field('EFF_WAVE')*1.e6 self.minwave_ft=np.min(self.wave_sc) self.maxwave_ft=np.max(self.wave_sc) self.nwave_sc = self.wave_sc.shape[0] self.nwave_ft = self.wave_ft.shape[0] # OPDC parameters if self.gdelay_in == True: self.opdc_time = opdc.field('TIME')/1e6 # in seconds self.opdc_piezo_offset = opdc.field('PIEZO_DL_OFFSET') # in VOLTS self.opdc_fddl_offset = opdc.field('VLTI_DL_OFFSET')*1e3 # FDDL (NOT VLTI) offsets converted in microns self.opdc_state = opdc.field('STATE') # IMAGING_DETECTOR parameters if detector_sc is not None: self.regname_sc = detector_sc.field('REGNAME') self.ports_sc = detector_sc.field('PORTS') else: self.regname_sc = '' if detector_ft is not None: self.regname_ft = detector_ft.field('REGNAME') self.ports_ft = detector_ft.field('PORTS') else: self.regname_ft = '' # OI_ARRAY parameters if oi_array is not None: self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1e6 # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # Output spectra if spectrum_sc is not None: self.spectrum_sc = [] # initialization self.time_sc = spectrum_sc.field('TIME')/1E6 # in seconds self.spectrum_sc = np.zeros((self.nregion, self.time_sc.shape[0], self.nwave_sc)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_sc[i-1,:,:] = np.mean(spectrum_sc.field('DATA'+str(i)), axis=2) if spectrum_ft is not None: self.spectrum_ft = [] # initialization self.time_ft = spectrum_ft.field('TIME')/1E6 # in seconds self.spectrum_ft = np.zeros((self.nregion , self.time_ft.shape[0], self.nwave_ft)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_ft[i-1,:,:] = np.mean(spectrum_ft.field('DATA'+str(i)), axis=2) # METROLOGY parameters self.time_met = None self.metrology_volt = None if metrology is not None: self.time_met = metrology.field('TIME')/1e6 # in seconds self.metrology_volt = metrology.field('VOLT') self.nframe_met = self.time_met.shape[0] # Number of frames if 'HIERARCH ESO DET2 NDIT' in self.header: self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] else: self.nframe_sc = 0 self.nframe_ft = self.spectrum_ft.shape[1] gravi_preproc.close() class P2vmreduced: #------------------------------- #Columns in the OI_VIS of the SC # #TARGET_ID : id listed in OI_TARGET #TIME [usec] : timestamp #MJD [day] : #INT_TIME [s] : integration time of this frame #VISDATA [e,e] : complex coherent flux spectra of SC in this frame #VISERR [e,e] : complex error on the coherent flux spectra #UCOORD [m] : uv-plane of this SC frame #VCOORD [m] : uv-plane of this SC frame #STA_INDEX : station index in the OI_ARRAY # #PHASE_MET_FC [rad] : unwrap FT-SC phase as computed by the DRS algorithm #PHASE_MET_TEL [rad] : unwrap FT-SC phase as computed by the DRS algorithm, mean of 4 diodes #OPD_MET_FC [m] : unwrap SC-FT delay as computed by the TAC algorihm #OPD_MET_TEL [m] : unwrap SC-FT delay as computed by the TAC algorihm, 4 diodes # #VISDATA_FT [e,e] : VISDATA spectra of FT (integrated in this SC frame) #PHASE_REF [rad] : reference phase, actually -1*arg{VISDATA_FT}, re-interpolated in the SC wavelength. # #NFRAME_FT : number of averaged FT frame in this SC frame #NFRAME_MET : number of averaged MET frame in this SC frame # #GDELAY [m] : real-time GD #SNR : real-time SNR #GDELAY_BOOT [m] : best GD estimate, accounting closing triangles #SNR_BOOT : best SNR estimate, accounting closing triangles # #V_FACTOR_FT : measured visibility loss on the FT #V_FACTOR : predicted visibility loss of this SC frame (re-interpolation of V_FACTOR_FT on the SC wavelengths) # #FRINGEDET_RATIO : fraction of FT frame accepted in this SC frame #REJECTION_FLAG : this frame is accepted/rejected # #-------------------------------- #Columns in the OI_FLUX of the SC # #TARGET_ID : id listed in OI_TARGET #TIME [usec] : timestamp #MJD [day] : #INT_TIME [s] : integration time of this frame #FLUX [e] : flux spectra #FLUXERR [e] : error on flux spectra #STA_INDEX : station index in the OI_ARRAY #TOTALFLUX_SC [e] : total flux of SC in this SC frame (integrated over spectrum) #TOTALFLUX_FT [e] : total flux of FT in this SC frame (integrated over spectrum) # #------------------------- #Columns in the OI_VIS_MET # #PHASE_TEL [rad] : 4 diodes x 4 beams phases at telescope, unwrap by pipeline algorithm (FT-SC) #PHASE_FC [rad] : 4 beams phases at combiner, unwrap by pipeline algorithm (FT-SC) #OPD_TEL [m] : 4 diodes x 4 beams OPD at telescope, unwrap by TAC algorithm (SC-FT) #OPD_FC [m] : 4 beams OPDs at telescope, unwrap by TAC algorithm (SC-FT) #FLAG_FC, FLAG_TEL: flags computed by TAC algorithm #VAMP_FC_FT, VAMP_FC_SC, VAMP_TEL_FT, VAMP_TEL_SC: Volt amplitude # #------------------------- #Columns in the OPDC table # #TIME [usec] : timestamp #PIEZO_DL_OFFSET #VLTI_DL_OFFSET #STATE : global fringe tracking state #STEPS : target phase modulation per baseline (scrambled), in units of pi/8 # #---------------------------------------- #HEADER information about target position # #FT.ROBJ.ALPHA #FT.ROBJ.DELTA #FT.ROBJ.EPOCH #FT.ROBJ.EQUINOX #FT.ROBJ.PARALLAX #FT.ROBJ.PMA #FT.ROBJ.PMD #FT.ROBJ.MAG #FT.ROBJ.DIAMETER #FT.ROBJ.NAME # #INS.SOBJ.X #INS.SOBJ.Y #INS.SOBJ.DIAMETER #INS.SOBJ.MAG #INS.SOBJ.NAME # #And during a SWAP: # #FT.ROBJ.ALPHA/DELTA += INS.SOBJ.X/Y #INS.SOBJ.XY *= -1 #INS.SOBJ.DIAMETER <=> INS.ROBJ.DIAMETER #INS.SOBJ.MAG <=> INS.ROBJ.MAG #INS.SOBJ.NAME <=> FT.ROBJ.NAME # #------------------- #Changes 2016-02-01: # #OI_VIS_MET #rename: PHI -> PHASE_FC [rad], 4 beams #rename: PHI_TEL -> PHASE_TEL, 4 diodes x 4 beams #change: OPD_TEL and OPD_FC to add the zero # #OI_VIS of the SC #rename: PHI_MET -> PHASE_MET_FC [rad], scalar #rename: PHI_MET_TEL -> PHASE_MET_TEL [rad], scalar, mean of four diodes #add: OPD_MET_FC [m] scalar #add: OPD_MET_TEL [m] four diodes #rename: EXP_PHASE -> PHASE_REF #change: GDELAY and GDELAY_BOOT -> now in [m] # def __init__(self, filename): # open fits file self.filename = filename gravi_p2vmreduced=pyfits.open(filename, mode='readonly') nbase = 6 # Main FITS header self.header = gravi_p2vmreduced[0].header # Single mode of the instrument (no OI_VIS_MET data) #if 'SINGLE' in str(self.header['HIERARCH ESO DO CATG']): if 'SINGLE' in str(self.header['HIERARCH ESO PRO REC1 RAW1 CATG']): self.single = True else: self.single = False if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True fddl = None # in case the FDDL data are not in the file oi_vis_met = None oi_flux_sc = None oi_flux_sc_s = None oi_vis_sc = None oi_vis_sc_s = None self.gdelay_in=False self.vfactor_in=False # get the FDDL, OI_ARRAY extentions for hdu in gravi_p2vmreduced : if hdu.name == 'OI_ARRAY' : oi_array=hdu.data if hdu.name == 'ARRAY_GEOMETRY' : # TBC oi_array=hdu.data if hdu.name == 'FDDL': fddl = hdu.data if hdu.name == 'OI_VIS_MET' : oi_vis_met=hdu.data if hdu.name == 'OPDC' : opdc=hdu.data self.gdelay_in = True # get the OI_VIS, OI_FLUX and T3 extensions in combined polarization mode if self.polarsplit == False: for hdu in gravi_p2vmreduced : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_vis_sc=hdu.data.copy() header_oi_vis_sc = hdu.header.copy() column_oi_vis_sc = hdu.columns.names # New table names 2015-11-03 if (hdu.name == 'OI_VIS') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') : oi_vis_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_vis2_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS2') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') : oi_vis2_ft=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_flux_sc=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc.field('TIME')[::4]/3600. # in decimal hours headerflux = hdu.header.copy() columnflux = hdu.columns.names # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') : oi_flux_ft=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_T3') & ('SPECTRO_SC' in str(hdu.header['INSNAME'])): oi_t3_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('SPECTRO_FT' in str(hdu.header['INSNAME'])): oi_t3_ft=hdu.data.copy() if (hdu.name == 'OI_T3_FT') : oi_t3_ft=hdu.data.copy() # get the OI_VIS, OI_FLUX and T3 extensions in split polarization mode (S considered =P1, P considered =P2) if self.polarsplit == True: for hdu in gravi_p2vmreduced : header = hdu.header.copy() if 'INSNAME' in header: # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_wave_sc=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_SC') : oi_wave_sc=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_WAVELENGTH') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_WAVELENGTH_FT') : oi_wave_ft=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_vis_sc_p=hdu.data.copy() header_oi_vis_sc_p = hdu.header.copy() column_oi_vis_sc_p = hdu.columns.names if (hdu.name == 'OI_VIS') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_vis_sc_s=hdu.data.copy() header_oi_vis_sc_s = hdu.header.copy() column_oi_vis_sc_s = hdu.columns.names if (hdu.name == 'OI_VIS2') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_vis2_sc_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_vis2_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_VIS') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis_ft_s=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_vis2_ft_p=hdu.data.copy() if (hdu.name == 'OI_VIS2_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_vis2_ft_s=hdu.data.copy() if (hdu.name == 'OI_FLUX') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_flux_sc_p=hdu.data.copy() self.oi_flux_sc_time = oi_flux_sc_p.field('TIME')[::4]/3600. # in decimal hours headerflux = hdu.header.copy() columnflux = hdu.columns.names if (hdu.name == 'OI_FLUX') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_flux_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_FLUX') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_FLUX_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_flux_ft_p=hdu.data.copy() self.oi_flux_ft_time = oi_flux_ft_p.field('TIME')[::4]/3600. # in decimal hours if (hdu.name == 'OI_FLUX_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_flux_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3') & ('SPECTRO_SC_P1' in str(hdu.header['INSNAME'])): oi_t3_sc_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('SPECTRO_SC_P2' in str(hdu.header['INSNAME'])): oi_t3_sc_s=hdu.data.copy() # New table names 2015-11-03 if (hdu.name == 'OI_T3') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('SPECTRO_FT_P1' in str(hdu.header['INSNAME'])): oi_t3_ft_p=hdu.data.copy() if (hdu.name == 'OI_T3_FT') & ('SPECTRO_FT_P2' in str(hdu.header['INSNAME'])): oi_t3_ft_s=hdu.data.copy() # Wavelength scales self.wave_sc = oi_wave_sc.field('EFF_WAVE')*1.e6 # wavelength scale in microns self.minwave_sc=np.min(self.wave_sc) self.maxwave_sc=np.max(self.wave_sc) self.nwave_sc = self.wave_sc.shape[0] self.wave_ft = oi_wave_ft.field('EFF_WAVE')*1.e6 self.minwave_ft=np.min(self.wave_sc) self.maxwave_ft=np.max(self.wave_sc) self.nwave_ft = self.wave_ft.shape[0] # OI_ARRAY parameters self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1e6 # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # OPDC parameters if self.gdelay_in == True: # print("Read OPDC data") self.opdc_time = opdc.field('TIME')/1e6 # in seconds self.opdc_piezo_offset = opdc.field('PIEZO_DL_OFFSET') # in VOLTS self.opdc_fddl_offset = opdc.field('VLTI_DL_OFFSET')*1e3 # FDDL (NOT VLTI) offsets converted in microns self.opdc_state = opdc.field('STATE') # OI_VIS_MET parameters if oi_vis_met is not None: print("Load the MET and TAC data") self.oi_vis_met_time = oi_vis_met.field('TIME')/1e6 # in seconds if 'PHI' in oi_vis_met.columns.names: self.oi_vis_met_phase = oi_vis_met.field('PHI') if 'PHASE_FC' in oi_vis_met.columns.names: self.oi_vis_met_phase = oi_vis_met.field('PHASE_FC') if 'PHI_TEL' in oi_vis_met.columns.names: self.oi_vis_met_phase_tel = oi_vis_met.field('PHI_TEL') if 'PHASE_TEL' in oi_vis_met.columns.names: self.oi_vis_met_phase_tel = oi_vis_met.field('PHASE_TEL') if 'OPD_TEL' in oi_vis_met.columns.names: self.met_opd_tel = oi_vis_met.field('OPD_TEL') if 'OPD_FC' in oi_vis_met.columns.names: self.met_opd_fc = oi_vis_met.field('OPD_FC') if 'FLAG_TEL' in oi_vis_met.columns.names: self.met_flag_tel = oi_vis_met.field('FLAG_TEL') if 'FLAG_FC' in oi_vis_met.columns.names: self.met_flag_fc = oi_vis_met.field('FLAG_FC') if 'VAMP_TEL_SC' in oi_vis_met.columns.names: self.met_vamp_tel_sc = oi_vis_met.field('VAMP_TEL_SC') if 'VAMP_TEL_FT' in oi_vis_met.columns.names: self.met_vamp_tel_ft = oi_vis_met.field('VAMP_TEL_FT') if 'VAMP_FC_SC' in oi_vis_met.columns.names: self.met_vamp_fc_sc = oi_vis_met.field('VAMP_FC_SC') if 'VAMP_FC_FT' in oi_vis_met.columns.names: self.met_vamp_fc_ft = oi_vis_met.field('VAMP_FC_FT') # DELAY to account for the missing first frame of the SC # SCdelay = self.header['HIERARCH ESO INS TIM1 PERIOD'] > now included in DRS # OI_FLUX table parameters if (oi_flux_sc is not None) & (self.polarsplit == False): print("Load FLUX data") self.oi_flux_sc_time = oi_flux_sc.field('TIME')[::4]/1e6 self.oi_flux_ft_time = oi_flux_ft.field('TIME')[::4]/1e6 if 'VISFLUX' in columnflux: self.oi_flux_sc = oi_flux_sc.field('VISFLUX').reshape(self.oi_flux_sc_time.shape[0],4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('VISFLUX').reshape(self.oi_flux_ft_time.shape[0],4,self.nwave_ft) else: self.oi_flux_sc = oi_flux_sc.field('FLUX').reshape(self.oi_flux_sc_time.shape[0],4,self.nwave_sc) self.oi_flux_ft = oi_flux_ft.field('FLUX').reshape(self.oi_flux_ft_time.shape[0],4,self.nwave_ft) if 'TOTALFLUX_SC' in columnflux: self.totalflux_sc = oi_flux_sc.field('TOTALFLUX_SC').reshape(self.oi_flux_sc_time.shape[0],4) self.totalflux_ft = oi_flux_sc.field('TOTALFLUX_FT').reshape(self.oi_flux_sc_time.shape[0],4) if (oi_flux_sc_s is not None) & (self.polarsplit == True): self.oi_flux_sc_time_s = oi_flux_sc_s.field('TIME')[::4]/1e6 self.oi_flux_ft_time_s = oi_flux_ft_s.field('TIME')[::4]/1e6 self.oi_flux_sc_time_p = oi_flux_sc_p.field('TIME')[::4]/1e6 self.oi_flux_ft_time_p = oi_flux_ft_p.field('TIME')[::4]/1e6 if 'VISFLUX' in columnflux: self.oi_flux_sc_s = oi_flux_sc_s.field('VISFLUX').reshape(self.oi_flux_sc_time_s.shape[0],4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('VISFLUX').reshape(self.oi_flux_ft_time_s.shape[0],4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('VISFLUX').reshape(self.oi_flux_sc_time_p.shape[0],4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('VISFLUX').reshape(self.oi_flux_ft_time_p.shape[0],4,self.nwave_ft) else: self.oi_flux_sc_s = oi_flux_sc_s.field('FLUX').reshape(self.oi_flux_sc_time_s.shape[0],4,self.nwave_sc) self.oi_flux_ft_s = oi_flux_ft_s.field('FLUX').reshape(self.oi_flux_ft_time_s.shape[0],4,self.nwave_ft) self.oi_flux_sc_p = oi_flux_sc_p.field('FLUX').reshape(self.oi_flux_sc_time_p.shape[0],4,self.nwave_sc) self.oi_flux_ft_p = oi_flux_ft_p.field('FLUX').reshape(self.oi_flux_ft_time_p.shape[0],4,self.nwave_ft) if 'TOTALFLUX_SC' in columnflux: self.totalflux_sc_s = oi_flux_sc_s.field('TOTALFLUX_SC').reshape(self.oi_flux_sc_time_s.shape[0],4) self.totalflux_ft_s = oi_flux_sc_s.field('TOTALFLUX_FT').reshape(self.oi_flux_sc_time_s.shape[0],4) self.totalflux_sc_p = oi_flux_sc_p.field('TOTALFLUX_SC').reshape(self.oi_flux_sc_time_p.shape[0],4) self.totalflux_ft_p = oi_flux_sc_p.field('TOTALFLUX_FT').reshape(self.oi_flux_sc_time_p.shape[0],4) # Number of frames self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] if self.polarsplit == False: self.nframe_ft = self.oi_flux_ft_time.shape[0] if self.polarsplit == True: self.nframe_ft = self.oi_flux_ft_time_s.shape[0] self.nframe_met = self.oi_vis_met_time.shape[0] # OI_VIS table parameters if (oi_vis_sc is not None) & (self.polarsplit == False): print("Load VIS data") self.time_sc = oi_vis_sc.field('TIME')[::nbase]/1e6 if 'EXP_PHASE' in oi_vis_sc.columns.names: self.oi_vis_sc_exp_phase = oi_vis_sc.field('EXP_PHASE') if 'PHASE_REF' in oi_vis_sc.columns.names: self.oi_vis_sc_exp_phase = oi_vis_sc.field('PHASE_REF') if 'PHI_MET' in oi_vis_sc.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc.field('PHI_MET') if 'PHASE_MET_FC' in oi_vis_sc.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc.field('PHASE_MET_FC') if 'PHASE_MET_TEL' in oi_vis_sc.columns.names: self.oi_vis_sc_phi_met_tel = oi_vis_sc.field('PHASE_MET_TEL') if 'OPD_MET_FC' in oi_vis_sc.columns.names: self.oi_vis_sc_opd_met = oi_vis_sc.field('OPD_MET_FC').reshape(self.nframe_sc,nbase).T*1e6 # in microns if 'OPD_MET_TEL' in oi_vis_sc.columns.names: self.oi_vis_sc_opd_met_tel = oi_vis_sc.field('OPD_MET_TEL')*1e6 # in microns self.oi_vis_sc_visdata = oi_vis_sc.field('VISDATA') self.oi_vis_sc_viserr = oi_vis_sc.field('VISERR') self.oi_vis_sc_visamp = oi_vis_sc.field('VISAMP') self.oi_vis_sc_visamperr = oi_vis_sc.field('VISAMPERR') self.oi_vis_sc_visphi = oi_vis_sc.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr = oi_vis_sc.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_ucoord = oi_vis_sc.field('UCOORD') self.oi_vis_sc_vcoord = oi_vis_sc.field('VCOORD') self.oi_vis_sc_sta_index = oi_vis_sc.field('STA_INDEX') if 'VISDATA_FT' in oi_vis_sc.columns.names: self.oi_vis_sc_visdataft = oi_vis_sc.field('VISDATA_FT') if 'VISVAR_FT' in oi_vis_sc.columns.names: self.oi_vis_sc_visvarft = oi_vis_sc.field('VISVAR_FT') self.oi_vis_sc_vispowerft = oi_vis_sc.field('VISPOWER_FT') self.oi_vis2_sc = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) for baseline in range(0,nbase): # self.oi_vis2_sc[baseline,:,:] = oi_vis2_sc.field('VIS2DATA')[baseline::nbase,:] self.oi_vis2_sc[baseline,:,:] = oi_vis_sc.field('VISAMP')[baseline::nbase,:]**2 * np.sign(oi_vis_sc.field('VISAMP')[baseline::nbase,:]) self.time_ft = oi_vis_ft.field('TIME')[::nbase]/1e6 self.oi_vis_ft_visdata = oi_vis_ft.field('VISDATA') self.oi_vis_ft_viserr = oi_vis_ft.field('VISERR') self.oi_vis_ft_visamp = oi_vis_ft.field('VISAMP') self.oi_vis_ft_visamperr = oi_vis_ft.field('VISAMPERR') self.oi_vis_ft_visphi = oi_vis_ft.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr = oi_vis_ft.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_ucoord = oi_vis_ft.field('UCOORD') self.oi_vis_ft_vcoord = oi_vis_ft.field('VCOORD') self.oi_vis_ft_sta_index = oi_vis_ft.field('STA_INDEX') if 'TARGET_PHASE' in oi_vis_ft.columns.names: self.oi_vis_ft_target_phase = oi_vis_ft.field('TARGET_PHASE') else: self.oi_vis_ft_target_phase = self.oi_vis_ft_visphi * 0.0 self.oi_vis2_ft = np.zeros((nbase,self.nframe_ft,self.nwave_ft)) for baseline in range(0,nbase): # self.oi_vis2_ft[baseline,:,:] = oi_vis2_ft.field('VIS2DATA')[baseline::nbase,:] # baseline, time, wavelength self.oi_vis2_ft[baseline,:,:] = oi_vis_ft.field('VISAMP')[baseline::nbase,:]**2 * np.sign(oi_vis_ft.field('VISAMP')[baseline::nbase,:]) # baseline, time, wavelength if self.gdelay_in == True: # Group delay in microns self.gdelay_sc = oi_vis_sc.field('GDELAY').reshape(self.nframe_sc,nbase).T*1e6 if 'GDELAY_FT' in oi_vis_sc.columns.names: self.gdelay_ft_resamp = oi_vis_sc.field('GDELAY_FT').reshape(self.nframe_sc,nbase).T*1e6 # GD of FT at SC frequency self.gdelay_ft = oi_vis_ft.field('GDELAY').reshape(self.nframe_ft,nbase).T*1e6 self.gdelay_boot_ft = oi_vis_ft.field('GDELAY_BOOT').reshape(self.nframe_ft,nbase).T*1e6 # Signal to noise of fringes self.snr_sc = oi_vis_sc.field('SNR').reshape(self.nframe_sc,nbase).T self.snr_ft = oi_vis_ft.field('SNR').reshape(self.nframe_ft,nbase).T self.snr_boot_ft = oi_vis_ft.field('SNR_BOOT').reshape(self.nframe_ft,nbase).T if 'V_FACTOR' in column_oi_vis_sc: print("Load VFACTOR data") self.vfactor_in = True self.oi_vis_vfactor_ft = oi_vis_sc.field('V_FACTOR_FT').reshape(self.nframe_sc,nbase,self.nwave_ft) self.oi_vis_vfactor = oi_vis_sc.field('V_FACTOR').reshape(self.nframe_sc,nbase,self.nwave_sc) self.oi_vis_reject_flag = oi_vis_sc.field('REJECTION_FLAG').reshape(self.nframe_sc,nbase) self.oi_vis_fringe_det_ratio = oi_vis_sc.field('FRINGEDET_RATIO').reshape(self.nframe_sc,nbase) if (oi_vis_sc_s is not None) & (self.polarsplit == True): self.time_sc = oi_vis_sc_s.field('TIME')[::nbase]/1e6 if 'EXP_PHASE' in oi_vis_sc_s.columns.names: self.oi_vis_sc_exp_phase_s = oi_vis_sc_s.field('EXP_PHASE') self.oi_vis_sc_exp_phase_p = oi_vis_sc_p.field('EXP_PHASE') if 'PHASE_REF' in oi_vis_sc_s.columns.names: self.oi_vis_sc_exp_phase_s = oi_vis_sc_s.field('PHASE_REF') self.oi_vis_sc_exp_phase_p = oi_vis_sc_p.field('PHASE_REF') if 'PHI_MET' in oi_vis_sc_s.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc_s.field('PHI_MET') # self.oi_vis_sc_phi_met_p = oi_vis_sc_p.field('PHI_MET') if 'PHASE_MET_FC' in oi_vis_sc_s.columns.names: self.oi_vis_sc_phi_met = oi_vis_sc_s.field('PHASE_MET_FC') # self.oi_vis_sc_phi_met_p = oi_vis_sc_p.field('PHASE_MET_FC') if 'PHASE_MET_TEL' in oi_vis_sc_s.columns.names: self.oi_vis_sc_phi_met_tel = oi_vis_sc_s.field('PHASE_MET_TEL') #self.oi_vis_sc_phi_met_tel_p = oi_vis_sc_p.field('PHASE_MET_TEL') if 'OPD_MET_FC' in oi_vis_sc_s.columns.names: self.oi_vis_sc_opd_met = oi_vis_sc_s.field('OPD_MET_FC').reshape(self.nframe_sc,nbase).T*1e6 # in microns #self.oi_vis_sc_opd_met_p = oi_vis_sc_p.field('OPD_MET_FC').reshape(self.nframe_sc,nbase).T*1e6 # in microns if 'OPD_MET_TEL' in oi_vis_sc_s.columns.names: self.oi_vis_sc_opd_met_tel = oi_vis_sc_s.field('OPD_MET_TEL')*1e6 # in microns #self.oi_vis_sc_opd_met_tel_p = oi_vis_sc_p.field('OPD_MET_TEL') self.oi_vis_sc_visdata_s = oi_vis_sc_s.field('VISDATA') self.oi_vis_sc_viserr_s = oi_vis_sc_s.field('VISERR') self.oi_vis_sc_visamp_s = oi_vis_sc_s.field('VISAMP') self.oi_vis_sc_visamperr_s = oi_vis_sc_s.field('VISAMPERR') self.oi_vis_sc_visphi_s = oi_vis_sc_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr_s = oi_vis_sc_s.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_ucoord_s = oi_vis_sc_s.field('UCOORD') self.oi_vis_sc_vcoord_s = oi_vis_sc_s.field('VCOORD') self.oi_vis_sc_sta_index_s = oi_vis_sc_s.field('STA_INDEX') if 'VISDATA_FT' in oi_vis_sc_s.columns.names: self.oi_vis_sc_visdataft_s = oi_vis_sc_s.field('VISDATA_FT') self.oi_vis_sc_visdataft_p = oi_vis_sc_p.field('VISDATA_FT') if 'VISVAR_FT' in oi_vis_sc_s.columns.names: self.oi_vis_sc_visvarft_s = oi_vis_sc_s.field('VISVAR_FT') self.oi_vis_sc_visvarft_p = oi_vis_sc_p.field('VISVAR_FT') self.oi_vis_sc_vispowerft_s = oi_vis_sc_s.field('VISPOWER_FT') self.oi_vis_sc_vispowerft_p = oi_vis_sc_p.field('VISPOWER_FT') self.oi_vis_sc_visdata_p = oi_vis_sc_p.field('VISDATA') self.oi_vis_sc_viserr_p = oi_vis_sc_p.field('VISERR') self.oi_vis_sc_visamp_p = oi_vis_sc_p.field('VISAMP') self.oi_vis_sc_visamperr_p = oi_vis_sc_p.field('VISAMPERR') self.oi_vis_sc_visphi_p = oi_vis_sc_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_visphierr_p = oi_vis_sc_p.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_sc_ucoord_p = oi_vis_sc_p.field('UCOORD') self.oi_vis_sc_vcoord_p = oi_vis_sc_p.field('VCOORD') self.oi_vis_sc_sta_index_p = oi_vis_sc_p.field('STA_INDEX') self.oi_vis2_sc_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) self.oi_vis2_sc_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) for baseline in range(0,nbase): # self.oi_vis2_sc_s[baseline,:,:] = oi_vis2_sc_s.field('VIS2DATA')[baseline::nbase,:] # baseline, time, wavelength # self.oi_vis2_sc_p[baseline,:,:] = oi_vis2_sc_p.field('VIS2DATA')[baseline::nbase,:] self.oi_vis2_sc_s[baseline,:,:] = oi_vis_sc_s.field('VISAMP')[baseline::nbase,:]**2 * np.sign(oi_vis_sc_s.field('VISAMP')[baseline::nbase,:]) # baseline, time, wavelength self.oi_vis2_sc_p[baseline,:,:] = oi_vis_sc_p.field('VISAMP')[baseline::nbase,:]**2 * np.sign(oi_vis_sc_p.field('VISAMP')[baseline::nbase,:]) self.time_ft = oi_vis_ft_s.field('TIME')[::nbase]/1e6 self.oi_vis_ft_visdata_s = oi_vis_ft_s.field('VISDATA') self.oi_vis_ft_viserr_s = oi_vis_ft_s.field('VISERR') self.oi_vis_ft_visamp_s = oi_vis_ft_s.field('VISAMP') self.oi_vis_ft_visamperr_s = oi_vis_ft_s.field('VISAMPERR') self.oi_vis_ft_visphi_s = oi_vis_ft_s.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr_s = oi_vis_ft_s.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_ucoord_s = oi_vis_ft_s.field('UCOORD') self.oi_vis_ft_vcoord_s = oi_vis_ft_s.field('VCOORD') self.oi_vis_ft_sta_index_s = oi_vis_ft_s.field('STA_INDEX') self.oi_vis_ft_visdata_p = oi_vis_ft_p.field('VISDATA') self.oi_vis_ft_viserr_p = oi_vis_ft_p.field('VISERR') self.oi_vis_ft_visamp_p = oi_vis_ft_p.field('VISAMP') self.oi_vis_ft_visamperr_p = oi_vis_ft_p.field('VISAMPERR') self.oi_vis_ft_visphi_p = oi_vis_ft_p.field('VISPHI')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_visphierr_p = oi_vis_ft_p.field('VISPHIERR')*np.pi/180. # conversion in radians (after Sept. 16) self.oi_vis_ft_ucoord_p = oi_vis_ft_p.field('UCOORD') self.oi_vis_ft_vcoord_p = oi_vis_ft_p.field('VCOORD') self.oi_vis_ft_sta_index_p = oi_vis_ft_p.field('STA_INDEX') if 'TARGET_PHASE' in oi_vis_ft_s.columns.names: self.oi_vis_ft_target_phase_s = oi_vis_ft_s.field('TARGET_PHASE') self.oi_vis_ft_target_phase_p = oi_vis_ft_p.field('TARGET_PHASE') else: self.oi_vis_ft_target_phase_s = self.oi_vis_ft_visphi_s * 0.0 self.oi_vis_ft_target_phase_p = self.oi_vis_ft_visphi_p * 0.0 self.oi_vis2_ft_s = np.zeros((nbase,self.nframe_ft,self.nwave_ft)) self.oi_vis2_ft_p = np.zeros((nbase,self.nframe_ft,self.nwave_ft)) for baseline in range(0,nbase): # self.oi_vis2_ft_s[baseline,:,:] = oi_vis2_ft_s.field('VIS2DATA')[baseline::nbase,:] # self.oi_vis2_ft_p[baseline,:,:] = oi_vis2_ft_p.field('VIS2DATA')[baseline::nbase,:] self.oi_vis2_ft_s[baseline,:,:] = oi_vis_ft_s.field('VISAMP')[baseline::nbase,:]**2 * np.sign(oi_vis_ft_s.field('VISAMP')[baseline::nbase,:]) self.oi_vis2_ft_p[baseline,:,:] = oi_vis_ft_p.field('VISAMP')[baseline::nbase,:]**2 * np.sign(oi_vis_ft_s.field('VISAMP')[baseline::nbase,:]) if self.gdelay_in == True: # Group delay in microns self.gdelay_sc_s = oi_vis_sc_s.field('GDELAY').reshape(self.nframe_sc,nbase).T*1e6 self.gdelay_sc_p = oi_vis_sc_p.field('GDELAY').reshape(self.nframe_sc,nbase).T*1e6 if 'GDELAY_FT' in oi_vis_sc_s.columns.names: self.gdelay_ft_resamp_s = oi_vis_sc_s.field('GDELAY_FT').reshape(self.nframe_sc,nbase).T*1e6 # GD of FT at SC frequency self.gdelay_ft_resamp_p = oi_vis_sc_p.field('GDELAY_FT').reshape(self.nframe_sc,nbase).T*1e6 # GD of FT at SC frequency self.gdelay_ft_s = oi_vis_ft_s.field('GDELAY').reshape(self.nframe_ft,nbase).T*1e6 # at FT frequency self.gdelay_ft_p = oi_vis_ft_p.field('GDELAY').reshape(self.nframe_ft,nbase).T*1e6 # at FT frequency # self.gdelay_boot_sc_s = oi_vis_sc_s.field('GDELAY_BOOT').reshape(self.nframe_sc,nbase).T*coh_length_sc # self.gdelay_boot_sc_p = oi_vis_sc_p.field('GDELAY_BOOT').reshape(self.nframe_sc,nbase).T*coh_length_sc self.gdelay_boot_ft_s = oi_vis_ft_s.field('GDELAY_BOOT').reshape(self.nframe_ft,nbase).T*1e6 self.gdelay_boot_ft_p = oi_vis_ft_p.field('GDELAY_BOOT').reshape(self.nframe_ft,nbase).T*1e6 # Signal to noise of fringes self.snr_sc_s = oi_vis_sc_s.field('SNR').reshape(self.nframe_sc,nbase).T self.snr_sc_p = oi_vis_sc_p.field('SNR').reshape(self.nframe_sc,nbase).T self.snr_ft_s = oi_vis_ft_s.field('SNR').reshape(self.nframe_ft,nbase).T self.snr_ft_p = oi_vis_ft_p.field('SNR').reshape(self.nframe_ft,nbase).T # self.snr_boot_sc_s = oi_vis_sc_s.field('SNR_BOOT').reshape(self.nframe_sc,nbase).T # self.snr_boot_sc_p = oi_vis_sc_p.field('SNR_BOOT').reshape(self.nframe_sc,nbase).T self.snr_boot_ft_s = oi_vis_ft_s.field('SNR_BOOT').reshape(self.nframe_ft,nbase).T self.snr_boot_ft_p = oi_vis_ft_p.field('SNR_BOOT').reshape(self.nframe_ft,nbase).T if 'V_FACTOR' in column_oi_vis_sc_p: self.vfactor_in = True self.oi_vis_vfactor_ft_s = oi_vis_sc_s.field('V_FACTOR_FT').reshape(self.nframe_sc,nbase,self.nwave_ft) # time_sc, baseline, wavelength_ft self.oi_vis_vfactor_s = oi_vis_sc_s.field('V_FACTOR').reshape(self.nframe_sc,nbase,self.nwave_sc) self.oi_vis_vfactor_ft_p = oi_vis_sc_p.field('V_FACTOR_FT').reshape(self.nframe_sc,nbase,self.nwave_ft) self.oi_vis_vfactor_p = oi_vis_sc_p.field('V_FACTOR').reshape(self.nframe_sc,nbase,self.nwave_sc) self.oi_vis_reject_flag_s = oi_vis_sc_s.field('REJECTION_FLAG').reshape(self.nframe_sc,nbase) self.oi_vis_reject_flag_p = oi_vis_sc_p.field('REJECTION_FLAG').reshape(self.nframe_sc,nbase) self.oi_vis_fringe_det_ratio_s = oi_vis_sc_s.field('FRINGEDET_RATIO').reshape(self.nframe_sc,nbase) self.oi_vis_fringe_det_ratio_p = oi_vis_sc_p.field('FRINGEDET_RATIO').reshape(self.nframe_sc,nbase) #============================================================================== # Computation of differential OPD quantities #============================================================================== if (self.gdelay_in == True) & (self.polarsplit == False): self.mean_gd_sc = np.zeros(nbase) self.mean_gd_ft = np.zeros(nbase) self.gd_diff = np.zeros((nbase,self.nframe_sc)) self.opd_diff_avg = np.zeros((nbase,self.nframe_sc)) self.opd_diff = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.total_dopd_avg = np.zeros((nbase,self.nframe_sc)) self.total_dopd = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC for baseline in range(0,nbase): # time average of GD over the sequence self.mean_gd_sc[baseline] = np.nanmean(self.gdelay_sc[baseline,:]) self.mean_gd_ft[baseline] = np.nanmean(self.gdelay_ft_resamp[baseline,:]) # Differential group delay in microns self.gd_diff[baseline,:] = self.gdelay_sc[baseline,:] - self.gdelay_ft_resamp[baseline,:] - (self.mean_gd_sc[baseline] - self.mean_gd_ft[baseline]) # Wavelength average of differential phase, in microns, once the mean GD has been subtracted visdata_sc = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_sc[baseline]/self.wave_sc)) for data in self.oi_vis_sc_visdata[baseline::nbase,:] ]) visdata_ft = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_ft[baseline]/self.wave_ft)) for data in self.oi_vis_sc_visdataft[baseline::nbase,:] ]) ph_diff_avg = visdata_sc * np.conj(visdata_ft) self.opd_diff_avg[baseline,:] = np.angle ( ph_diff_avg * np.conj(np.mean(ph_diff_avg)) ) / (2.*np.pi) * np.mean(self.wave_sc) # Differential phase in microns, function of wavelength visdata_sc = self.oi_vis_sc_visdata[baseline::nbase,:] * np.exp(1.j*self.oi_vis_sc_exp_phase[baseline::nbase,:]) for wave in range(0,self.nwave_sc): # self.opd_diff[baseline,:,wave] = np.angle( visdata_sc[:,wave] * np.conj(visdata_sc[:,wave]) ) / (2.*np.pi) * self.wave_sc[wave] self.opd_diff[baseline,:,wave] = np.angle( visdata_sc[:,wave] * np.conj(np.mean(visdata_sc[:,wave])) ) / (2.*np.pi) * self.wave_sc[wave] # Construction of the total dOPD SC-FT including metrology and GD self.total_dopd_avg[baseline,:] = self.opd_diff_avg[baseline,:] + (self.mean_gd_sc[baseline] - self.mean_gd_ft[baseline]) + self.oi_vis_sc_opd_met[baseline,:] # Construction of the total dOPD SC-FT including metrology and GD as a function of wavelength for wave in range(0,self.nwave_sc): self.total_dopd[baseline,:, wave] = self.opd_diff[baseline,:,wave] + (self.mean_gd_sc[baseline] - self.mean_gd_ft[baseline]) + self.oi_vis_sc_opd_met[baseline,:] if (self.gdelay_in == True) & (self.polarsplit == True): self.mean_gd_sc_s = np.zeros(nbase) self.mean_gd_sc_p = np.zeros(nbase) self.mean_gd_ft_s = np.zeros(nbase) self.mean_gd_ft_p = np.zeros(nbase) self.gd_diff_s = np.zeros((nbase,self.nframe_sc)) self.gd_diff_p = np.zeros((nbase,self.nframe_sc)) self.opd_diff_avg_s = np.zeros((nbase,self.nframe_sc)) self.opd_diff_avg_p = np.zeros((nbase,self.nframe_sc)) self.opd_diff_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.opd_diff_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.total_dopd_avg_s = np.zeros((nbase,self.nframe_sc)) self.total_dopd_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC self.total_dopd_avg_p = np.zeros((nbase,self.nframe_sc)) self.total_dopd_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc)) # as a function of wavelength channel of SC for baseline in range(0,nbase): self.mean_gd_sc_s[baseline] = np.nanmean(self.gdelay_sc_s[baseline,:]) self.mean_gd_sc_p[baseline] = np.nanmean(self.gdelay_sc_p[baseline,:]) self.mean_gd_ft_s[baseline] = np.nanmean(self.gdelay_ft_resamp_s[baseline,:]) self.mean_gd_ft_p[baseline] = np.nanmean(self.gdelay_ft_resamp_p[baseline,:]) # Differential group delay in microns self.gd_diff_s[baseline,:] = self.gdelay_sc_s[baseline,:] - self.gdelay_ft_resamp_s[baseline,:] - (self.mean_gd_sc_s[baseline] - self.mean_gd_ft_s[baseline]) self.gd_diff_p[baseline,:] = self.gdelay_sc_p[baseline,:] - self.gdelay_ft_resamp_p[baseline,:] - (self.mean_gd_sc_p[baseline] - self.mean_gd_ft_p[baseline]) # Differential phase average in the band, in microns, once the mean_gd has been removed visdata_sc_s = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_sc_s[baseline]/self.wave_sc)) for data in self.oi_vis_sc_visdata_s[baseline::nbase,:] ]) visdata_sc_p = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_sc_p[baseline]/self.wave_sc)) for data in self.oi_vis_sc_visdata_p[baseline::nbase,:] ]) visdata_ft_s = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_ft_s[baseline]/self.wave_ft)) for data in self.oi_vis_sc_visdataft_s[baseline::nbase,:] ]) visdata_ft_p = np.array([ np.nanmean(data * np.exp(-2.j*np.pi*self.mean_gd_ft_p[baseline]/self.wave_ft)) for data in self.oi_vis_sc_visdataft_p[baseline::nbase,:] ]) ph_diff_avg_s = visdata_sc_s * np.conj(visdata_ft_s); ph_diff_avg_p = visdata_sc_p * np.conj(visdata_ft_p); self.opd_diff_avg_s[baseline,:] = np.angle ( ph_diff_avg_s * np.conj(np.mean(ph_diff_avg_s)) ) / (2.*np.pi) * np.mean(self.wave_sc) self.opd_diff_avg_p[baseline,:] = np.angle ( ph_diff_avg_p * np.conj(np.mean(ph_diff_avg_p)) ) / (2.*np.pi) * np.mean(self.wave_sc) # Differential phase in microns, function of wavelength visdata_sc_s = self.oi_vis_sc_visdata_s[baseline::nbase,:] * np.exp (1.j*self.oi_vis_sc_exp_phase_s[baseline::nbase,:]) visdata_sc_p = self.oi_vis_sc_visdata_p[baseline::nbase,:] * np.exp (1.j*self.oi_vis_sc_exp_phase_p[baseline::nbase,:]) for wave in range(0,self.nwave_sc): self.opd_diff_s[baseline,:,wave] = np.angle ( visdata_sc_s[:,wave] * np.conj(np.mean(visdata_sc_s[:,wave])) ) / (2.*np.pi) * self.wave_sc[wave] self.opd_diff_p[baseline,:,wave] = np.angle ( visdata_sc_p[:,wave] * np.conj(np.mean(visdata_sc_p[:,wave])) ) / (2.*np.pi) * self.wave_sc[wave] # Construction of the total dOPD SC-FT including metrology and GD self.total_dopd_avg_s[baseline,:] = self.opd_diff_avg_s[baseline,:] + (self.mean_gd_sc_s[baseline] - self.mean_gd_ft_s[baseline]) + self.oi_vis_sc_opd_met[baseline,:] self.total_dopd_avg_p[baseline,:] = self.opd_diff_avg_p[baseline,:] + (self.mean_gd_sc_p[baseline] - self.mean_gd_ft_p[baseline]) + self.oi_vis_sc_opd_met[baseline,:] # Construction of the total dOPD SC-FT including metrology and GD as a function of wavelength for wave in range(0,self.nwave_sc): self.total_dopd_s[baseline,:, wave] = self.opd_diff_s[baseline,:,wave] + (self.mean_gd_sc_s[baseline] - self.mean_gd_ft_s[baseline]) + self.oi_vis_sc_opd_met[baseline,:] self.total_dopd_p[baseline,:, wave] = self.opd_diff_p[baseline,:,wave] + (self.mean_gd_sc_p[baseline] - self.mean_gd_ft_p[baseline]) + self.oi_vis_sc_opd_met[baseline,:] # #============================================================================== # # Computation of the group delay signals # #============================================================================== # # # Smoothing of the FT group delay signals # smooth_gd_ft = 30 # number of samples for smoothing # # if (self.polarsplit == False): # nwave_sc = self.nwave_sc # nwave_ft = self.nwave_ft # # # SC group delay signals # visdata_sc = np.zeros((nbase,self.nframe_sc,self.nwave_sc),dtype=complex) # for base in range(0,nbase): # for wave in range(0,self.nwave_sc): # visdata_sc[base,:,wave] = self.oi_vis_sc_visdata[base::nbase,wave] # # lambda_mean_sc = np.mean(self.wave_sc) # delta_lambda_sc = self.wave_sc[1]-self.wave_sc[0] # self.group_delay_sc = np.zeros((nbase,self.nframe_sc)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_sc): # incstep_sc = visdata_sc[base,frame,0:nwave_sc-1]*np.conjugate(visdata_sc[base,frame,1:nwave_sc]) # self.group_delay_sc[base,frame] = np.angle(np.sum(incstep_sc))/(2*np.pi)*(lambda_mean_sc**2)/delta_lambda_sc # # # FT group delay signals # visdata_ft = np.zeros((nbase,self.nframe_ft,self.nwave_ft),dtype=complex) # for base in range(0,nbase): # for wave in range(0,self.nwave_ft): # visdata_ft[base,:,wave] = self.oi_vis_ft_visdata[base::nbase,wave] # # lambda_mean_ft = np.mean(self.wave_ft) # delta_lambda_ft = self.wave_ft[1]-self.wave_ft[0] # incstep_ft = np.zeros((nbase,self.nframe_ft),dtype=complex) # self.group_delay_ft = np.zeros((nbase,self.nframe_ft)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # incstep_ft[base,frame] = np.sum(visdata_ft[base,frame,0:nwave_ft-1]*np.conjugate(visdata_ft[base,frame,1:nwave_ft])) # # for base in range(0,nbase): # smoothing of the inter-spectra incstep # incstep_ft[base,:] = np.convolve(incstep_ft[base,:], np.ones((smooth_gd_ft,))/smooth_gd_ft, mode='same') # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # self.group_delay_ft[base,frame] = np.angle(incstep_ft[base,frame])/(2*np.pi)*(lambda_mean_ft**2)/delta_lambda_ft # # if (self.polarsplit == True): # nwave_sc = self.nwave_sc # nwave_ft = self.nwave_ft # # # SC group delay signals # visdata_sc_s = np.zeros((nbase,self.nframe_sc,self.nwave_sc),dtype=complex) # visdata_sc_p = np.zeros((nbase,self.nframe_sc,self.nwave_sc),dtype=complex) # for base in range(0,nbase): # for wave in range(0,self.nwave_sc): # visdata_sc_s[base,:,wave] = self.oi_vis_sc_visdata_s[base::nbase,wave] # visdata_sc_p[base,:,wave] = self.oi_vis_sc_visdata_p[base::nbase,wave] # # lambda_mean_sc = np.mean(self.wave_sc) # delta_lambda_sc = self.wave_sc[1]-self.wave_sc[0] # coh_length_sc = (lambda_mean_sc**2)/delta_lambda_sc # # self.group_delay_sc_s = np.zeros((nbase,self.nframe_sc)) # self.group_delay_sc_p = np.zeros((nbase,self.nframe_sc)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_sc): # incstep_sc_s = visdata_sc_s[base,frame,0:nwave_sc-1]*np.conjugate(visdata_sc_s[base,frame,1:nwave_sc]) # incstep_sc_p = visdata_sc_p[base,frame,0:nwave_sc-1]*np.conjugate(visdata_sc_p[base,frame,1:nwave_sc]) # self.group_delay_sc_s[base,frame] = np.angle(np.sum(incstep_sc_s))/(2*np.pi)*coh_length_sc # self.group_delay_sc_p[base,frame] = np.angle(np.sum(incstep_sc_p))/(2*np.pi)*coh_length_sc # # # FT group delay signals # visdata_ft_s = np.zeros((nbase,self.nframe_ft,self.nwave_ft),dtype=complex) # visdata_ft_p = np.zeros((nbase,self.nframe_ft,self.nwave_ft),dtype=complex) # # for base in range(0,nbase): # for wave in range(0,self.nwave_ft): # visdata_ft_s[base,:,wave] = self.oi_vis_ft_visdata_s[base::nbase,wave] # visdata_ft_p[base,:,wave] = self.oi_vis_ft_visdata_p[base::nbase,wave] # # lambda_mean_ft = np.mean(self.wave_ft) # delta_lambda_ft = self.wave_ft[1]-self.wave_ft[0] # coh_length_ft = (lambda_mean_ft**2)/delta_lambda_ft # # incstep_ft_s = np.zeros((nbase,self.nframe_ft),dtype=complex) # incstep_ft_p = np.zeros((nbase,self.nframe_ft),dtype=complex) # self.group_delay_ft_s = np.zeros((nbase,self.nframe_ft)) # self.group_delay_ft_p = np.zeros((nbase,self.nframe_ft)) # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # incstep_ft_s[base,frame] = np.sum(visdata_ft_s[base,frame,0:nwave_ft-1]*np.conjugate(visdata_ft_s[base,frame,1:nwave_ft])) # incstep_ft_p[base,frame] = np.sum(visdata_ft_p[base,frame,0:nwave_ft-1]*np.conjugate(visdata_ft_p[base,frame,1:nwave_ft])) # # for base in range(0,nbase): # smoothing of the inter-spectra incstep # incstep_ft_s[base,:] = np.convolve(incstep_ft_s[base,:], np.ones((smooth_gd_ft,))/smooth_gd_ft, mode='same') # incstep_ft_p[base,:] = np.convolve(incstep_ft_p[base,:], np.ones((smooth_gd_ft,))/smooth_gd_ft, mode='same') # # for base in range(0,nbase): # for frame in range(0,self.nframe_ft): # self.group_delay_ft_s[base,frame] = np.angle(incstep_ft_s[base,frame])/(2*np.pi)*coh_length_ft # self.group_delay_ft_p[base,frame] = np.angle(incstep_ft_p[base,frame])/(2*np.pi)*coh_length_ft # #============================================================================== # # Computation of the K factors # #============================================================================== # # k2factor = np.zeros((nbase,self.nframe_sc)) # dt_sc = self.sc_time[nbase]-self.sc_time[0] # dt_ft = self.ft_time[nbase]-self.ft_time[0] # # if self.polarsplit == True: # # FT frames before the first SC integration (5s) # skip_ft = (self.sc_time[0]-dt_sc/2) // dt_ft # visdata_ft_s = np.mean(self.oi_vis_ft_visdata_s.reshape(self.ft_time.shape[0],6,self.nwave_ft),axis=2) # time, baseline, wavelength # visdata_ft_p = np.mean(self.oi_vis_ft_visdata_p.reshape(self.ft_time.shape[0],6,self.nwave_ft),axis=2) # time, baseline, wavelength # visdata_ft = np.mean([visdata_ft_s, visdata_ft_p],axis=0) # else: # dt_sc = self.oi_vis_sc_time[nbase]-self.oi_vis_sc_time_s[0] # dt_ft = self.oi_vis_ft_time[nbase]-self.oi_vis_ft_time_s[0] # # FT frames before the first SC integration (5s) # skip_ft = (self.oi_vis_sc_time[0]-dt_sc/2) // dt_ft # visdata_ft = np.mean(self.oi_vis_ft_visdata.reshape(self.ft_time.shape[0],6,self.nwave_ft),axis=2) # time, baseline, wavelength # # # Ratio of the inter-frame time between SC and FT # ratio_sc_ft = dt_sc/dt_ft # # # for baseline in range(0,nbase): # for frame_sc in range(0,self.nframe_sc): # minindex_ft = int(ratio_sc_ft * frame_sc)+skip_ft # maxindex_ft = minindex_ft + int(ratio_sc_ft) # # bias_coh = np.var(np.abs(visdata_ft[minindex_ft:maxindex_ft,baseline])) # v2_coh = np.mean(np.abs(np.square(visdata_ft[minindex_ft:maxindex_ft,baseline]))) # v2_incoh = np.square(np.abs(np.mean(visdata_ft[minindex_ft:maxindex_ft,baseline]))) # k2factor[baseline, frame_sc] = v2_coh/v2_incoh # # self.kfactor = np.sqrt(k2factor) # # print self.kfactor gravi_p2vmreduced.close() class Rawdata: def __init__(self, filename): # open fits file self.filename = filename gravi_raw = pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_raw[0].header # Single mode of the instrument (no metrology data) # if 'SINGLE' in str(self.header['HIERARCH ESO DO CATG']): # self.single = True # else: # self.single = False # Raw data types if 'HIERARCH ESO DPR TYPE' in self.header: dtype = str(self.header['HIERARCH ESO DPR TYPE']) self.catg = dtype # default case # if 'DARK' in dtype: self.catg = 'DARK' # if 'FLAT' in dtype: self.catg = 'FLAT' # if 'WAVE' in dtype: self.catg = 'WAVE' # if 'P2VM' in dtype: self.catg = 'P2VM' else: self.catg = 'UNKNOWN' # Polarization setup if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False else: self.polarsplit = True fddl = None # in case the FDDL data are not in the file metrology = None detector_sc = None detector_ft = None self.gdelay_in = False # Number of SC frames self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] self.exptime_sc = self.header['HIERARCH ESO DET2 SEQ1 DIT'] from astropy.time import Time def mjd_date(datestring): datestring2 = datestring[:10]+' '+datestring[12:] t = Time([datestring2],format='iso',scale='utc') return t.mjd[0] # get the ACQ, ARRAY_GEOMETRY, FDDL, OPDC, imaging_detector, imaging_data and METROLOGY extentions for hdu in gravi_raw : if hdu.name == 'ARRAY_GEOMETRY' : # array geometry oi_array=hdu.data if hdu.name == 'IMAGING_DATA_ACQ' : # acquisition camera frames self.acq = hdu.data # direct assignation #if hdu.name == 'OPDC': # opdc = hdu.data if hdu.name == 'FDDL': fddl = hdu.data if hdu.name == 'IMAGING_DETECTOR_SC' : detector_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_FT' : detector_ft=hdu.data # IMAGING_DATA SC frames direct assignation if hdu.name == 'IMAGING_DATA_SC' : self.data_sc = hdu.data self.time_sc = np.array(list(range(0,self.nframe_sc))) * self.header['HIERARCH ESO INS TIM1 PERIOD'] + \ mjd_date(self.header['HIERARCH ESO INS TIM1 START'])*86400. # reference in the middle of the inter-frame if hdu.name == 'IMAGING_DATA_FT' : imaging_ft=hdu.data if hdu.name == 'METROLOGY' : metrology=hdu.data if hdu.name == 'OPDC' : opdc=hdu.data self.gdelay_in = True # OPDC parameters if self.gdelay_in == True: self.opdc_time = opdc.field('TIME')/1e6 # in seconds self.opdc_piezo_offset = opdc.field('PIEZO_DL_OFFSET') # in VOLTS self.opdc_fddl_offset = opdc.field('VLTI_DL_OFFSET')*1e3 # FDDL (NOT VLTI) offsets converted in microns self.opdc_state = opdc.field('STATE') # OI_ARRAY parameters self.array_tel_name = oi_array.field('TEL_NAME') self.array_sta_name = oi_array.field('STA_NAME') self.array_staxyz = oi_array.field('STAXYZ') # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1e6 + mjd_date(self.header['HIERARCH ESO PCR ACQ START'])*86400. # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # OPDC table # TBD # IMAGING_DETECTOR parameters if detector_sc is not None: self.regname_sc = detector_sc.field('REGNAME') else: self.regname_sc = '' if detector_ft is not None: self.regname_ft = detector_ft.field('REGNAME') else: self.regname_ft = '' # IMAGING_DATA_FT self.time_ft = imaging_ft.field('TIME')/1e6 + mjd_date(self.header['HIERARCH ESO PCR ACQ START'])*86400. # in seconds self.exptime_ft = self.header['HIERARCH ESO DET3 SEQ1 DIT'] if self.polarsplit == True: #self.data_ft = imaging_ft.field('PIX').reshape(self.time_ft.shape[0],5,48) data = imaging_ft.field('PIX')#.reshape(self.time_ft.shape[0],24,10) # print data.shape self.data_ft = np.zeros((self.time_ft.shape[0],48,5)) for i in range(0,24): for j in range(0,5): self.data_ft[:,2*i,j] = data[:,j,i] for j in range(5,10): self.data_ft[:,2*i+1,j-5] = data[:,j,i] # raw.data_ft[time,output,wave_channel] # print self.data_ft else: self.data_ft = imaging_ft.field('PIX').reshape(self.time_ft.shape[0],24,5) self.nframe_ft = self.time_ft.shape[0] # number of FT frames # print self.data_ft # METROLOGY parameters if metrology is not None: self.metrology_time = metrology.field('TIME')/1e6 + mjd_date(self.header['HIERARCH ESO PCR ACQ START'])*86400. # in seconds self.metrology_volt = metrology.field('VOLT') self.metrology_powerlaser = metrology.field('POWER_LASER') self.metrology_lambdalaser = metrology.field('LAMBDA_LASER') gravi_raw.close() class Badpix: def __init__(self, filename): # open fits file self.filename = filename gravi_badpix=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_badpix[0].header for hdu in gravi_badpix : if hdu.name=='IMAGING_DATA_SC' : h = hdu.header self.ncol_sc= h['NAXIS1'] self.nrow_sc= h['NAXIS2'] if hdu.name=='IMAGING_DATA_FT' : h = hdu.header self.nregion_ft= h['TFIELDS']-1 self.nwave_ft= 5 # Image of the badpix map of the SC detector self.badpiximage_sc = gravi_badpix['IMAGING_DATA_SC'].data.squeeze() # squeeze to remove the spurious dimension self.badpix_exptime_sc = gravi_badpix['IMAGING_DATA_SC'].header['EXPTIME'] badpixtable_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') badpix_ft = gravi_badpix['IMAGING_DATA_FT'].data for output in range(1,self.nregion_ft+1): badpixtable_ft[output-1,:] = badpix_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft) self.badpixtable_ft = badpixtable_ft gravi_badpix.close() class Profile: def __init__(self, filename): # open fits file self.filename = filename gravi_profile=pyfits.open(filename, mode='readonly') # Main header self.header = gravi_profile[0].header if "COMBINED" in str(self.header['HIERARCH ESO FT POLA MODE']): self.polarsplit = False self.nchannels = 24 else: self.polarsplit = True self.nchannels = 48 imgaging_det_sc = gravi_profile['IMAGING_DETECTOR_SC'].data self.regname_sc = imgaging_det_sc.field('REGNAME') self.center_sc = (imgaging_det_sc.field('CENTER')) imgaging_det_ft = gravi_profile['IMAGING_DETECTOR_FT'].data self.regname_ft = imgaging_det_ft.field('REGNAME') self.center_ft = (imgaging_det_ft.field('CENTER')) # get the PROFILE_DATA useful keywords h = gravi_profile['PROFILE_DATA'].header self.startx = h['STARTX'] self.nx = h['NX'] del h self.imaging_data_sc = np.squeeze(gravi_profile['IMAGING_DATA_SC'].data) self.imaging_data_ft = np.squeeze(gravi_profile['IMAGING_DATA_FT'].data.field('PIX')) # DOES NOT WORK ON MY MAC FOR UNKNOWN REASON # self.profile_data = [] #np.zeros((self.nchannels)) # prof = gravi_profile['PROFILE_DATA'].data # print ' *** Reading profiles...' # for i in range(0,self.nchannels): # self.profile_data.append(np.squeeze(prof['DATA'+str(i+1)])) # print ' *** Reading profiles: OK !' gravi_profile.close() del gravi_profile class Dark: def __init__(self, filename): # open fits file self.filename = filename gravi_dark=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_dark[0].header for hdu in gravi_dark : if hdu.name=='IMAGING_DATA_SC' : h = hdu.header self.ncol_sc= h['NAXIS1'] self.ncol_sc= h['NAXIS2'] if hdu.name=='IMAGING_DATA_FT' : h = hdu.header self.nregion_ft= h['TFIELDS']-1 self.nwave_ft= 5 # Image of the dark of the SC detector self.darkimage_sc = gravi_dark['IMAGING_DATA_SC'].data.squeeze() # squeeze to remove the spurious dimension self.dark_exptime_sc = gravi_dark['IMAGING_DATA_SC'].header['EXPTIME'] for hdu in gravi_dark : if hdu.name=='IMAGING_ERR_SC' : # now dark_error is in RMS self.darkerror_sc = (gravi_dark['IMAGING_ERR_SC'].data.squeeze()) darktable_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') dark_ft = gravi_dark['IMAGING_DATA_FT'].data for output in range(1,self.nregion_ft+1): darktable_ft[output-1,:] = dark_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft) self.darktable_ft = darktable_ft for hdu in gravi_dark : # the ERR tables contain the variance if hdu.name=='IMAGING_ERR_FT' : darkerror_ft = np.zeros([self.nregion_ft,self.nwave_ft],dtype='d') dark_ft = gravi_dark['IMAGING_ERR_FT'].data for output in range(1,self.nregion_ft+1): # now dark_error is in RMS darkerror_ft[output-1,:] = (np.abs(dark_ft.field('DATA'+str(output))[0].reshape(1,self.nwave_ft))) self.darkerror_ft = darkerror_ft gravi_dark.close() class PhaseFTSC: def __init__(self, filename): # open fits file self.filename = filename gravi_phaseftsc=pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_phaseftsc[0].header spectrum_sc = None spectrum_ft = None detector_sc = None detector_ft = None # get the FDDL, OI_ARRAY, OI_WAVELENGTH and OI_VIS_MET extentions for hdu in gravi_phaseftsc : if hdu.name == 'SPECTRUM_DATA_FT' : spectrum_ft=hdu.data if hdu.name == 'SPECTRUM_DATA_SC' : spectrum_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_SC' : detector_sc=hdu.data if hdu.name == 'IMAGING_DETECTOR_FT' : detector_ft=hdu.data if hdu.name == 'FDDL': fddl = hdu.data self.fddlpresent = True # FDDL positions if fddl is not None: self.fddl_time = fddl.field('TIME')/1e6 # in seconds self.fddl_ft_pos = fddl.field('FT_POS') self.fddl_sc_pos = fddl.field('SC_POS') self.fddl_opl_air = fddl.field('OPL_AIR') # # Pipeline FT and SC average phase per baseline # ext7 = gravi_phaseftsc[7].data # ext8 = gravi_phaseftsc[8].data # self.met_phase_sc = ext7.field('PHASE') # baseline, time_sc (expressed in radians) # self.met_phase_ft = ext8.field('PHASE') # baseline, time_ft (expressed in radians) # Metrology OPL SC-FT tables per telescope from pipeline ext9 = gravi_phaseftsc[9].data self.time_met = ext9.field('TIME')/1E6 # in seconds self.nframe_met = self.time_met.shape[0] self.met_dopl = np.zeros((4,self.nframe_met)) for tel in range(0,4): self.met_dopl[tel,:] = ext9.field('OPL'+str(tel+1))*1e6 # tel, time_met IN MICRONS # IMAGING_DETECTOR parameters if detector_sc is not None: self.regname_sc = detector_sc.field('REGNAME') self.ports_sc = detector_sc.field('PORTS') else: self.regname_sc = '' if detector_ft is not None: self.regname_ft = detector_ft.field('REGNAME') self.ports_ft = detector_ft.field('PORTS') else: self.regname_ft = '' self.nregion = self.regname_sc.shape[0] self.nwave_sc = spectrum_sc.field('DATA1').shape[1] self.nwave_ft = spectrum_ft.field('DATA1').shape[1] # Extracted non resampled spectra if spectrum_sc is not None: self.spectrum_sc = [] # initialization self.time_sc = spectrum_sc.field('TIME')/1E6 # in seconds self.spectrum_sc = np.zeros((self.nregion, self.time_sc.shape[0], self.nwave_sc)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_sc[i-1,:,:] = np.mean(spectrum_sc.field('DATA'+str(i)), axis=2) if spectrum_ft is not None: self.spectrum_ft = [] # initialization self.time_ft = spectrum_ft.field('TIME')/1E6 # in seconds self.spectrum_ft = np.zeros((self.nregion, self.time_ft.shape[0], self.nwave_ft)) # output, time, wavelength for i in range(1,self.nregion+1): self.spectrum_ft[i-1,:,:] = np.mean(spectrum_ft.field('DATA'+str(i)), axis=2) # Number of frames if 'HIERARCH ESO DET2 NDIT' in self.header: self.nframe_sc = self.header['HIERARCH ESO DET2 NDIT'] else: self.nframe_sc = 0 self.nframe_ft = self.spectrum_ft.shape[1] gravi_phaseftsc.close() class Argon: def __init__(self, filename): # open fits file self.filename = filename gravi_argon = pyfits.open(filename, mode='readonly') # Main FITS header self.header = gravi_argon[0].header self.frames = np.array(gravi_argon['IMAGING_DATA_SC'].data)