Source code for antaress.ANTARESS_process.ANTARESS_main

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
from pathos.multiprocessing import cpu_count
from copy import deepcopy
import bindensity as bind
from os import makedirs
from os.path import exists as path_exist
import glob
from astropy.io import fits
from dace_query.spectroscopy import Spectroscopy
from scipy.interpolate import CubicSpline
from ..ANTARESS_analysis.ANTARESS_model_prof import calc_macro_ker_anigauss,calc_macro_ker_rt 
from ..ANTARESS_grids.ANTARESS_star_grid import model_star
from ..ANTARESS_grids.ANTARESS_occ_grid import occ_region_grid,sub_calc_plocc_ar_prop,calc_plocc_ar_prop,retrieve_ar_prop_from_param,generate_ar_prop
from ..ANTARESS_grids.ANTARESS_prof_grid import init_custom_DI_prof,custom_DI_prof,theo_intr2loc,gen_theo_atm,var_stellar_prop
from ..ANTARESS_grids.ANTARESS_coord import calc_mean_anom_TR,calc_Kstar,calc_tr_contacts,calc_rv_star,get_timeorbit,coord_expos,coord_expos_ar
from ..ANTARESS_analysis.ANTARESS_inst_resp import return_pix_size,def_st_prof_tab,cond_conv_st_prof_tab,conv_st_prof_tab,get_FWHM_inst,resamp_st_prof_tab,return_spec_nord,flag_err_inst,return_cen_wav_ord
from ..ANTARESS_analysis.ANTARESS_ana_comm import par_formatting_inst_vis
from ..ANTARESS_general.minim_routines import par_formatting
from ..ANTARESS_corrections.ANTARESS_sp_reduc import red_sp_data_instru
from ..ANTARESS_analysis.ANTARESS_joined_star import joined_Star_ana
from ..ANTARESS_analysis.ANTARESS_joined_atm import joined_Atm_ana
from ..ANTARESS_plots.ANTARESS_plots_all import ANTARESS_plot_functions
from ..ANTARESS_corrections.ANTARESS_calib import calc_gcal
from ..ANTARESS_process.ANTARESS_plocc_spec import def_in_plocc_profiles,def_diff_profiles,eval_diff_profiles
from ..ANTARESS_conversions.ANTARESS_masks_gen import def_masks
from ..ANTARESS_conversions.ANTARESS_conv import DI_CCF_from_spec,DiffIntr_CCF_from_spec,Atm_CCF_from_spec,conv_2D_to_1D_spec
from ..ANTARESS_conversions.ANTARESS_binning import process_bin_prof
from ..ANTARESS_corrections.ANTARESS_detrend import detrend_prof,pc_analysis
from ..ANTARESS_process.ANTARESS_data_process import align_profiles,rescale_profiles,extract_diff_profiles,extract_intr_profiles,extract_pl_profiles,EvE_outputs 
from ..ANTARESS_analysis.ANTARESS_ana_comm import MAIN_single_anaprof
from ..ANTARESS_conversions.ANTARESS_sp_cont import process_spectral_cont
from ..ANTARESS_general.utils import air_index,dataload_npz,gen_specdopshift,stop,np_where1D,closest,datasave_npz,def_edge_tab,check_data,npint,path_singslash
from ..ANTARESS_general.constant_data import Rsun,Rjup,c_light,G_usi,Msun,AU_1


[docs] def ANTARESS_settings_overwrite(data_dic,mock_dic,gen_dic,theo_dic,plot_dic,glob_fit_dic,detrend_prof_dic,input_dic): r"""**ANTARESS settings overwrite.** Overwrites ANTARESS default settings with custom inputs. Args: TBD Returns: None """ #Overwriting full dictionaries if 'gen_dic' in input_dic['settings']:gen_dic.update(input_dic['settings']['gen_dic']) if 'mock_dic' in input_dic['settings']:mock_dic.update(input_dic['settings']['mock_dic']) if 'data_dic' in input_dic['settings']: for key in ['DI','Intr']: if key in input_dic['settings']['data_dic']:data_dic[key].update(input_dic['settings']['data_dic'][key]) if 'glob_fit_dic' in input_dic['settings']: for key in ['DIProp','IntrProf','IntrProp','DiffProf']: if key in input_dic['settings']['glob_fit_dic']:glob_fit_dic[key].update(input_dic['settings']['glob_fit_dic'][key]) if 'plot_dic' in input_dic['settings']:plot_dic.update(input_dic['settings']['plot_dic']) if 'detrend_prof_dic' in input_dic['settings']:detrend_prof_dic.update(input_dic['settings']['detrend_prof_dic']) #Overwriting specific fields if ('orders4ccf' in input_dic) and (len(input_dic['orders4ccf'])>0): for inst in input_dic['orders4ccf']:gen_dic['orders4ccf'][inst] = input_dic['orders4ccf'][inst] return None
[docs] def ANTARESS_main(data_dic,mock_dic,gen_dic,theo_dic,plot_dic,glob_fit_dic,detrend_prof_dic,system_param,input_dic,custom_plot_settings): r"""**Main ANTARESS function.** Runs ANTARESS workflow. The pipeline is defined as modules than can be run independently. Each module takes as input the datasets produced by earlier modules, transforms or analyzes them, and saves the outputs to disk. This approach allows the user to re-run the pipeline from any module, which is useful when several flow options are available. It is even critical with large datasets such as the ones produced by ESPRESSO, which can take several hours to process with a given module. Finally, this approach also allows a user to retrieve data at any stage of the process flow for external use. Args: TBD Returns: None """ ############################################################################## #Initializations ############################################################################## coord_dic,data_prop = init_gen(data_dic,mock_dic,gen_dic,system_param,theo_dic,plot_dic,glob_fit_dic,detrend_prof_dic) #################################################################################################################### #Processing datasets for each visit of each instrument # - binned visits are processed at the end #################################################################################################################### if gen_dic['sequence'] not in ['system_view']: for inst in data_dic['instrum_list']: print('') print('-----------------------') print('Processing instrument :',inst) print('-----------------------') #Initialize instrument tables and dictionaries init_inst(mock_dic,inst,gen_dic,data_dic,theo_dic,data_prop,coord_dic,system_param,plot_dic) #Spectral data if ('spec' in data_dic[inst]['type']): #Estimating instrumental calibration if gen_dic['gcal']: calc_gcal(gen_dic,data_dic,inst,plot_dic,coord_dic,data_prop) #Global corrections of spectral data # - performed before the loop on individual visits because some corrections exploit information from all visits and require the full range of the data red_sp_data_instru(inst,data_dic,plot_dic,gen_dic,data_prop,coord_dic,system_param) #------------------------------------------------- #Processing all visits for current instrument for vis in data_dic[inst]['visit_list']: print(' -----------------') print(' Processing visit: '+vis) print(' -----------------') #Initialization of visit properties init_vis(data_prop,data_dic,vis,coord_dic,inst,system_param,gen_dic) #------------------------------------------------- #Processing disk-integrated stellar profiles data_type_gen = 'DI' #------------------------------------------------- #Spectral detrending if gen_dic['detrend_prof'] and (detrend_prof_dic['full_spec']): detrend_prof(detrend_prof_dic,data_dic,coord_dic,inst,vis,data_dic,data_prop,gen_dic,plot_dic) #Converting stellar spectra into CCFs if gen_dic[data_type_gen+'_CCF']: DI_CCF_from_spec(inst,vis,data_dic,gen_dic) #Single line detrending if gen_dic['detrend_prof'] and (not detrend_prof_dic['full_spec']): detrend_prof(detrend_prof_dic,data_dic,coord_dic,inst,vis,data_dic,data_prop,gen_dic,plot_dic) #Calculating theoretical properties of the planet-occulted and/or active regions if gen_dic['theoPlOcc'] and ( (len(data_dic[inst][vis]['studied_pl'])>0) or (len(data_dic[inst][vis]['studied_ar'])>0) ): calc_plocc_ar_prop(system_param,gen_dic,theo_dic,coord_dic,inst,vis,data_dic,calc_pl_atm=gen_dic['calc_pl_atm'],ar_dic=theo_dic) #Analyzing original disk-integrated profiles if gen_dic['fit_'+data_type_gen]: MAIN_single_anaprof('',data_type_gen+'orig',data_dic,gen_dic,inst,vis,coord_dic,theo_dic,plot_dic,system_param['star']) #Aligning disk-integrated profiles to star rest frame if (gen_dic['align_'+data_type_gen]): align_profiles(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic) #Rescaling profiles to their correct flux level # - the module is used to calculate light curves even if no scaling is needed if gen_dic['flux_sc']: rescale_profiles(data_dic[inst],inst,vis,data_dic,coord_dic,coord_dic[inst][vis]['t_dur_d'],gen_dic,plot_dic,system_param,theo_dic,ar_dic=theo_dic) #Calculating master spectrum of the disk-integrated star used in weighted averages and continuum-normalization if gen_dic['DImast_weight']: process_bin_prof('',data_type_gen,gen_dic,inst,vis,data_dic,coord_dic,data_prop,system_param,theo_dic,plot_dic,masterDIweigh=True,ar_dic=theo_dic) #Converting 2D disk-integrated profiles into 1D, and analyzing them if gen_dic['spec_1D']: conv_2D_to_1D_gen_functions(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic,theo_dic,plot_dic,system_param) #Processing binned disk-integrated profiles if gen_dic['bin']: bin_gen_functions(data_type_gen,'',inst,gen_dic,data_dic,coord_dic,data_prop,system_param,theo_dic,plot_dic,vis=vis) #-------------------------------------------------------------------------------------------------- #Processing differential profiles data_type_gen = 'Diff' #-------------------------------------------------------------------------------------------------- #Extracting differential profiles if (gen_dic['diff_data']): extract_diff_profiles(gen_dic,data_dic,inst,vis,data_prop,coord_dic) #Converting differential spectra into CCFs if gen_dic[data_type_gen+'_CCF']: DiffIntr_CCF_from_spec(data_type_gen,inst,vis,data_dic,gen_dic) #Converting and analyzing 2D differential spectra into 1D if gen_dic['spec_1D']: conv_2D_to_1D_gen_functions(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic,theo_dic,plot_dic,system_param) #-------------------------------------------------------------------------------------------------- #Processing intrinsic stellar profiles data_type_gen = 'Intr' #-------------------------------------------------------------------------------------------------- #Extracting intrinsic stellar profiles if gen_dic['intr_data']: extract_intr_profiles(data_dic,gen_dic,inst,vis,system_param['star'],coord_dic,theo_dic,plot_dic) #Converting out-of-transit differential and intrinsic spectra into CCFs if gen_dic[data_type_gen+'_CCF']: DiffIntr_CCF_from_spec(data_type_gen,inst,vis,data_dic,gen_dic,data_type_gen) #Applying PCA to out-of transit differential profiles if (gen_dic['pca_ana']): pc_analysis(gen_dic,data_dic,inst,vis,data_prop,coord_dic) #Analyzing intrinsic stellar profiles in the star rest frame if gen_dic['fit_'+data_type_gen]: MAIN_single_anaprof('',data_type_gen+'orig',data_dic,gen_dic,inst,vis,coord_dic,theo_dic,plot_dic,system_param['star']) #Aligning intrinsic stellar profiles to their local rest frame if gen_dic['align_'+data_type_gen]: align_profiles(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic) #Converting and analyzing 2D out-of-transit differential and intrinsic spectra into 1D, and analyzing them if gen_dic['spec_1D']: conv_2D_to_1D_gen_functions(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic,theo_dic,plot_dic,system_param) #Processing binned intrinsic profiles if gen_dic['bin']: bin_gen_functions(data_type_gen,'',inst,gen_dic,data_dic,coord_dic,data_prop,system_param,theo_dic,plot_dic,vis=vis) #Building estimates for planet-occulted stellar profiles in in-transit exposures if gen_dic['loc_prof_est']: def_in_plocc_profiles(inst,vis,gen_dic,data_dic,data_prop,coord_dic,system_param,theo_dic,glob_fit_dic,plot_dic) #Building estimates for differential profiles in all exposures # - in-transit profiles include planet-occulted and active region stellar profiles if gen_dic['diff_prof_est']: def_diff_profiles(inst,vis,gen_dic,data_dic,data_prop,coord_dic,system_param,theo_dic,glob_fit_dic,plot_dic) #Correcting exposures from (occulted/non-occulted) active region contamination using the differential profile estimates if gen_dic['corr_diff']: eval_diff_profiles(inst, vis, gen_dic, data_dic,data_prop,coord_dic,system_param,theo_dic,glob_fit_dic,plot_dic,'corr') #Building best-fit differential profile time series if gen_dic['eval_bestfit']: eval_diff_profiles(inst, vis, gen_dic, data_dic,data_prop,coord_dic,system_param,theo_dic,glob_fit_dic,plot_dic,'bestfit') #-------------------------------------------------------------------------------------------------- #Processing atmospheric profiles data_type_gen = 'Atm' #-------------------------------------------------------------------------------------------------- #Extracting atmospheric profiles if gen_dic['pl_atm']: extract_pl_profiles(data_dic,inst,vis,gen_dic) #Converting atmospheric spectra into CCFs if gen_dic[data_type_gen+'_CCF'] and ('spec' in data_dic[inst][vis]['type']): Atm_CCF_from_spec(inst,vis,data_dic,gen_dic) #Fitting atmospheric profiles in the star rest frame if gen_dic['fit_'+data_type_gen]: MAIN_single_anaprof('',data_type_gen+'orig',data_dic,gen_dic,inst,vis,coord_dic,theo_dic,plot_dic,system_param['star']) #Aligning atmospheric profiles to the planet rest frame if gen_dic['align_'+data_type_gen]: align_profiles(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic) #Converted 2D atmospheric profiles into 1D, and analyzing them if gen_dic['spec_1D']: conv_2D_to_1D_gen_functions(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic,theo_dic,plot_dic,system_param) #Processing binned atmospheric profiles if gen_dic['bin']: bin_gen_functions(data_type_gen,'',inst,gen_dic,data_dic,coord_dic,data_prop,system_param,theo_dic,plot_dic,vis=vis) ### end of visits #Update instrument properties update_inst(data_dic,inst,gen_dic) #Processing binned profiles over multiple visits # - beware that data from different visits should be comparable to be binned # this is not the case, e.g, with blazed 2D spectra or if the stellar line shape changed if (data_dic[inst]['n_visits_inst']>1) and (gen_dic['binmultivis']): print(' --------------------------') print(' Processing combined visits') print(' --------------------------') for data_type_gen in ['DI','Intr','Atm']: bin_gen_functions(data_type_gen,'multivis',inst,gen_dic,data_dic,coord_dic,data_prop,system_param,theo_dic,plot_dic) ### end of instruments #Update generic properties if gen_dic['spec_1D'] or gen_dic['CCF_from_sp']: update_gen(data_dic,gen_dic) #Saving EvE outputs if gen_dic['EvE_outputs']: EvE_outputs(gen_dic,data_dic) #################################################################################################################### #Call to analysis function over combined instruments and/or visits #################################################################################################################### if gen_dic['joined_ana']: print('-------------------------------') print('Analyzing combined datasets') print('-------------------------------') #Wrap-up function to fit stellar profiles and their properties if gen_dic['fit_DIProp'] or gen_dic['fit_IntrProf'] or gen_dic['fit_IntrProp'] or gen_dic['fit_DiffProf'] : joined_Star_ana(glob_fit_dic,system_param,theo_dic,data_dic,gen_dic,plot_dic,coord_dic,data_prop) #Wrap-up function to fit atmospheric profiles and their properties if gen_dic['fit_AtmProf'] or gen_dic['fit_AtmProp']: joined_Atm_ana(gen_dic) ############################################################################## #Call to plot functions ############################################################################## if gen_dic['plots_on']: ANTARESS_plot_functions(system_param,plot_dic,data_dic,gen_dic,coord_dic,theo_dic,data_prop,glob_fit_dic,mock_dic,input_dic,custom_plot_settings) return None
[docs] def conv_2D_to_1D_gen_functions(data_type_gen,data_dic,inst,vis,gen_dic,coord_dic,theo_dic,plot_dic,system_param): """**Wrap-up function for 2D->1D datasets.** Args: TBD: Returns: None """ #Processing 2D spectra into new 1D spectra if gen_dic['spec_1D_'+data_type_gen] and (data_dic[inst][vis]['type']=='spec2D'): conv_2D_to_1D_spec(data_type_gen,inst,vis,gen_dic,data_dic,data_dic[data_type_gen],coord_dic) #Analyzing converted profiles if (data_type_gen!='Diff') and gen_dic['fit_'+data_type_gen+'_1D']: MAIN_single_anaprof('',data_type_gen+'_1D',data_dic,gen_dic,inst,vis,coord_dic,theo_dic,plot_dic,system_param['star']) return None
[docs] def bin_gen_functions(data_type_gen,mode,inst,gen_dic,data_dic,coord_dic,data_prop,system_param,theo_dic,plot_dic,mock_dic={},vis=None): """**Wrap-up function for binned datasets.** Args: TBD: Returns: None """ #Binning profiles for analysis purpose if gen_dic[data_type_gen+'bin'+mode]: process_bin_prof(mode,data_type_gen,gen_dic,inst,vis,data_dic,coord_dic,data_prop,system_param,theo_dic,plot_dic,ar_dic=mock_dic) #Analyzing binned profiles if gen_dic['fit_'+data_type_gen+'bin'+mode]: MAIN_single_anaprof(mode,data_type_gen+'bin',data_dic,gen_dic,inst,vis,coord_dic,theo_dic,plot_dic,system_param['star']) #Calculating generic stellar continuum from binned master spectrum if (data_type_gen in ['DI','Intr']) and gen_dic[data_type_gen+'_stcont']: process_spectral_cont(mode,data_type_gen,inst,data_dic,gen_dic,vis) #Defining CCF mask # - over all visits if possible, or over the single processed visit if (data_type_gen in ['DI','Intr']) and gen_dic['def_'+data_type_gen+'masks']: if (mode=='' and (data_dic[inst]['n_visits_inst']==1)) or (mode=='multivis'): def_masks(mode,gen_dic,data_type_gen,inst,vis,data_dic,plot_dic,system_param,data_prop) return None
[docs] def init_gen(data_dic,mock_dic,gen_dic,system_param,theo_dic,plot_dic,glob_fit_dic,detrend_prof_dic): r"""**Initialization: generic** Initializes generic fields for the workflow. Args: TBD Returns: None """ #Positional dictionary coord_dic={} #Additional properties data_prop={} #Data types gen_dic['type_name']={'DI':'disk-integrated','Diff':'differential','Intr':'intrinsic','Atm':'atmospheric','Absorption':'absorption','Emission':'emission'} gen_dic['type_list'] = np.array(['DI','Diff','Intr','Atm']) #Planets with known orbits and transit properties # - this does not imply that the planet is actually transiting gen_dic['def_pl']=[] for obj in system_param: if (obj!='star') and ('inclination' in system_param[obj]) and (obj in data_dic['DI']['system_prop']['achrom']):gen_dic['def_pl']+=[obj] #Deactivate all plot modules if requested if not gen_dic['plots_on']: for key in plot_dic: if plot_dic[key] in ['png','pdf','jpg']:plot_dic[key]='' #Check for active plots to activate call to global plot function if gen_dic['plots_on']: gen_dic['plots_on'] = False for key in plot_dic: if plot_dic[key] in ['png','pdf','jpg']: gen_dic['plots_on'] = True break #------------------------------------------------------------------------------ #Data is processed by the workflow if gen_dic['sequence'] not in ['system_view']: #Multi-threading print('Multi-threading: '+str(cpu_count())+' threads available') #Instruments to which are associated the different datasets if gen_dic['mock_data']: data_dic['instrum_list']=list(mock_dic['visit_def'].keys()) if len(data_dic['instrum_list'])==0:stop('ERROR : no instrument was requested for mock datasets.') print('Running with artificial data') else: data_dic['instrum_list'] = list(gen_dic['data_dir_list'].keys()) if len(data_dic['instrum_list'])==0:stop('ERROR : no instrument was requested for observed datasets in gen_dic["data_dir_list"].') print('Running with observational data') #Used visits for inst in data_dic['instrum_list']: if inst not in gen_dic['unused_visits']:gen_dic['unused_visits'][inst]=[] #Total number of processed instruments gen_dic['n_instru'] = len(data_dic['instrum_list']) #All processed types # - if not set, assumed to be CCF by default # - gen_dic['type'] traces the original types of datasets throughout the pipeline gen_dic['all_types'] = [] for inst in data_dic['instrum_list']: if (inst not in gen_dic['type']):gen_dic['type'][inst]='CCF' gen_dic['all_types']+=[gen_dic['type'][inst]] gen_dic['all_types'] = list(np.unique(gen_dic['all_types'])) gen_dic['specINtype'] = any('spec' in s for s in gen_dic['all_types']) gen_dic['ccfINtype'] = ('CCF' in gen_dic['all_types']) #Use of covariance gen_dic['flag_err_inst'] = {} if gen_dic['mock_data']:gen_dic['use_cov'] = False else: #Used by default with spectral datasets, unless CCF datasets are processed or user requests otherwise if gen_dic['ccfINtype']:gen_dic['use_cov'] = False if not gen_dic['use_cov']:print('Covariances discounted') #Automatic activation/deactivation if gen_dic['pca_ana']:gen_dic['intr_data'] = True if gen_dic['intr_data']: if not gen_dic['diff_data']: print('Automatic activation of differential profile extraction') gen_dic['diff_data'] = True else: for key in ['map_Intr_prof','all_intr_data','Intr_prof']:plot_dic[key]='' #Deactivate general conditions if (not gen_dic['specINtype']): if gen_dic['DI_CCF']: print('WARNING: all datasets in CCF mode, automatic disabling of DI profiles conversion into CCFs.') gen_dic['DI_CCF']=False if detrend_prof_dic['full_spec']: print('WARNING: all datasets in CCF mode, automatic disabling of spectral detrending.') detrend_prof_dic['full_spec'] = False if gen_dic['Diff_CCF']: if (not gen_dic['specINtype']) or (gen_dic['DI_CCF']) or (not gen_dic['diff_data']) or (gen_dic['intr_data']): print('WARNING: automatic disabling of Diff profiles conversion into CCFs.') gen_dic['Diff_CCF']=False if gen_dic['Intr_CCF']: if (not gen_dic['specINtype']) or (gen_dic['DI_CCF']) or (not gen_dic['intr_data']): print('WARNING: automatic disabling of Intr profiles conversion into CCFs.') gen_dic['Intr_CCF']=False #Deactivate all calculation options if not gen_dic['calc_all']: for key in gen_dic: if ('calc_' in key) and (gen_dic[key]):gen_dic[key] = False #Deactivate spectral corrections in CCF mode sp_corr_list = ['corr_tell','glob_mast','corr_FbalOrd','corr_Fbal','corr_Ftemp','corr_cosm','mask_permpeak','corr_wig','corr_fring','trim_spec'] if (not gen_dic['specINtype']): for key in sp_corr_list:gen_dic[key]=False gen_dic['gcal']=False gen_dic['cond_plot_gcal'] = False #Deactivate EvE outputs if gen_dic['EvE_outputs']: print('WARNING (EvE outputs): data must be in spectral mode') gen_dic['EvE_outputs'] = False else: #1D spectra if ('spec2D' not in gen_dic['all_types']): #Deactive if spectra are not 2D: # - spectral calibration, spectral balance correction over orders, wiggles gen_dic['gcal'] = False gen_dic['cond_plot_gcal'] = False gen_dic['corr_FbalOrd']=False plot_dic['wig_order_fit']='' #2D spectra else: #Automatic calculation of calibration profiles gen_dic['gcal'] = True gen_dic['cond_plot_gcal'] = (plot_dic['gcal_all']!='') or (plot_dic['gcal_ord']!='') or (plot_dic['noises_ord']!='') #Activate global master calculation if required for flux balance corrections if (gen_dic['corr_Fbal']) or (gen_dic['corr_FbalOrd']): gen_dic['glob_mast']=True #Default bin size for global flux balance (in 1e13 s-1) if (inst not in gen_dic['Fbal_bin_nu']): Fbal_bin_nu_inst = { 'SOPHIE_HE':1., 'SOPHIE_HR':1., 'CORALIE':1., 'HARPN':1., 'HARPS':1., 'ESPRESSO':1., 'ESPRESSO_MR':1., 'CARMENES_VIS':1., 'EXPRES':1., 'IGRINS2_Red':1., 'IGRINS2_Blue':1., 'MAROONX_Red':1., 'MAROONX_Blue':1., 'MIKE_Red':1., 'MIKE_Blue':1., 'NIGHT':1., 'NIRPS_HA':1., 'NIRPS_HE':1.} if inst not in Fbal_bin_nu_inst:stop('ERROR : default value for "gen_dic["Fbal_bin_nu"]" undefined for '+inst+' (source : ANTARESS_main.py > init_gen()') gen_dic['Fbal_bin_nu'][inst] = Fbal_bin_nu_inst[inst] #Default bin size for order flux balance (in A) if (inst not in gen_dic['FbalOrd_binw']): FbalOrd_binw_inst = { 'SOPHIE_HE':2., 'SOPHIE_HR':2., 'CORALIE':2., 'HARPN':2., 'HARPS':2., 'ESPRESSO':2., 'ESPRESSO_MR':2., 'CARMENES_VIS':2., 'EXPRES':2., 'NIGHT':2., 'NIRPS_HA':2., 'NIRPS_HE':2.} if inst not in FbalOrd_binw_inst:stop('ERROR : default value for "gen_dic["FbalOrd_binw"]" undefined for '+inst+' (source : ANTARESS_main.py > init_gen()') gen_dic['FbalOrd_binw'][inst] = FbalOrd_binw_inst[inst] #Deactivate arbitrary scaling if scaling not required if (not gen_dic['corr_Fbal']) and (not gen_dic['corr_FbalOrd']):gen_dic['Fbal_vis']=None #Check for arbitrary scaling elif (gen_dic['Fbal_vis']!='ext') and (gen_dic['n_instru']>1):stop('Flux balance must be set to a theoretical reference for multiple instruments') #Deactivate transit scaling if temporal flux correction is applied if gen_dic['corr_Ftemp']:data_dic['DI']['rescale_DI'] = True #Activate correction flag gen_dic['corr_data'] = False for key in sp_corr_list: gen_dic['corr_data'] |=gen_dic[key] #Deactivate modules and fields irrelevant to mock datasets if gen_dic['mock_data']: for key in ['gcal','corr_tell','tell_weight','glob_mast','corr_Fbal','corr_FbalOrd','corr_Ftemp','corr_cosm','mask_permpeak','corr_wig','EvE_outputs']:gen_dic[key] = False #Deactivate plot if module is not called if (plot_dic['glob_mast']!='') and gen_dic['glob_mast']: plot_dic['glob_mast']=='' print('Disabling "glob_mast" plot') if (plot_dic['Fbal_corr_vis']!='') and ((not gen_dic['corr_Fbal']) or (gen_dic['Fbal_vis'] is None)): plot_dic['Fbal_corr_vis']=='' print('Disabling "Fbal_corr_vis" plot') #Set general condition for conversions gen_dic['spec_1D'] = False gen_dic['calc_spec_1D'] = False gen_dic['CCF_from_sp'] = False gen_dic['calc_CCF_from_sp'] = False gen_dic['type2var'] = {'DI':'EFsc2','Diff':'EFdiff2','Intr':'EFintr2','Emission':'EFem2','Absorption':'EAbs2'} gen_dic['earliertypes4var'] = {'DI':['DI'],'Diff':['DI','Diff'],'Intr':['DI','Diff','Intr'],'Atm':['DI','Diff','Atm']} for key in gen_dic['type_list']: gen_dic['spec_1D'] |= gen_dic['spec_1D_'+key] gen_dic['calc_spec_1D'] |= gen_dic['calc_spec_1D_'+key] gen_dic['CCF_from_sp'] |= gen_dic[key+'_CCF'] gen_dic['calc_CCF_from_sp'] |= gen_dic['calc_'+key+'_CCF'] #Initialize flags data_dic[key]['type'] = {inst:deepcopy(gen_dic['type'][inst]) for inst in gen_dic['type']} data_dic[key]['spec_to_CCF'] = {inst:False for inst in gen_dic['type']} data_dic[key]['spec2D_to_spec1D'] = {inst:False for inst in gen_dic['type']} #Deactive 1D plots if no 1D conversion is applied if not gen_dic['spec_1D_'+key] and (plot_dic['sp_'+key+'_1D']!=''): print('WARNING: no 1D conversion of '+key+' profiles is applied, disabling plot_dic["sp_'+key+'_1D"]') plot_dic['sp_'+key+'_1D'] = '' if gen_dic['CCF_from_sp']:gen_dic['ccfINtype'] = True #Raise warning if gen_dic['EvE_outputs']: if not gen_dic['spec_1D_DI']:print('WARNING (EvE outputs): 1D conversion must be applied (gen_dic["spec_1D_DI"] = True)') #Set specific fields per instrument gen_dic['corr_FbalOrd_inst'] = {} for inst in gen_dic['type']: #Set specific conversion flags and final data mode for each type of profiles and instruments # - the general instrument and visit dictionaries contain a field type that represents the mode of the data at the current stage of the pipeline # - here we set the final mode for each type of profile, so that they can be retrieved in their original mode for plotting # eg, if profiles are converted into CCF at the differential stage, we keep the information that disk-integrated profiles were in spectral mode # - we also set a flag for each instrument related to the effective conversion of its datasets for ikey,key in enumerate(gen_dic['type_list']): #Conversion of current type into CCF is requested if gen_dic[key+'_CCF']: #Checking global data format if gen_dic['type'][inst]=='CCF':print('WARNING: input profiles for '+inst+' are in CCF format, no CCF conversion of '+key+' profiles will be applied.') else: #Checking format of previous data types conv_ok = True for iprev in range(0,ikey): #Stop if one of previous data types was already converted if data_dic[gen_dic['type_list'][iprev]]['spec_to_CCF'][inst]: conv_ok = False print('WARNING: '+gen_dic['type_name'][gen_dic['type_list'][iprev]]+' profiles for '+inst+' were converted in CCF format, no CCF conversion of '+gen_dic['type_name'][key]+' profiles will be applied.') break #Flag conversion as ok and affect data format to current and later types if conv_ok: #Flag is set to True only for the specific data type that was converted data_dic[key]['spec_to_CCF'][inst]=True for isub in range(ikey,len(gen_dic['type_list'])):data_dic[gen_dic['type_list'][isub]]['type'][inst]='CCF' #Conversion of Intr profiles is also applied to OT Diff profiles if key=='Intr': data_dic['Diff']['spec_to_CCF'][inst]=True data_dic['Diff']['type'][inst]='CCF' #Conversion of current type into 1D spectra is requested if gen_dic['spec_1D_'+key]: if gen_dic['type'][inst]=='CCF':print('WARNING: input profiles for '+inst+' are in CCF format, no 1D conversion of '+key+' profiles will be applied.') elif gen_dic['type'][inst]=='spec1D':print('WARNING: input profiles for '+inst+' are in 1D format, no 1D conversion of '+key+' profiles will be applied.') else: conv_ok = True for iprev in range(0,ikey): if data_dic[gen_dic['type_list'][iprev]]['spec2D_to_spec1D'][inst]: conv_ok = False print('WARNING: '+gen_dic['type_name'][gen_dic['type_list'][iprev]]+' profiles for '+inst+' were converted in 1D format, no 1D conversion of '+gen_dic['type_name'][key]+' profiles will be applied.') break if conv_ok: data_dic[key]['spec2D_to_spec1D'][inst]=True for isub in range(ikey,len(gen_dic['type_list'])):data_dic[gen_dic['type_list'][isub]]['type'][inst]='spec1D' if key=='Intr': data_dic['Diff']['spec2D_to_spec1D'][inst]=True data_dic['Diff']['type'][inst]='spec1D' #Set general condition to activate binning modules for mode in ['','multivis']: gen_dic['bin'+mode]=False gen_dic['calc_bin'+mode]=False for data_type_gen in ['DI','Intr','Atm']: gen_dic['bin'+mode] |= ( gen_dic[data_type_gen+'bin'+mode] | gen_dic['fit_'+data_type_gen+'bin'+mode]) gen_dic['calc_bin'+mode] |= ( gen_dic['calc_'+data_type_gen+'bin'+mode] | gen_dic['calc_fit_'+data_type_gen+'bin'+mode]) #Set general condition for single profiles fits for key in ['DI','Intr','Atm']: gen_dic['fit_'+key+'_gen'] = gen_dic['fit_'+key] | gen_dic['fit_'+key+'bin'] | gen_dic['fit_'+key+'binmultivis'] gen_dic['fit_Diff_gen'] = False #Automatic continuum and fit range for key in ['DI','Intr','Atm']: if gen_dic['fit_'+key+'_gen']: if ('cont_range' not in data_dic[key]):data_dic[key]['cont_range']={} if ('fit_range' not in data_dic[key]):data_dic[key]['fit_range']={} if (gen_dic['diff_data']) and ('cont_range' not in data_dic['Diff']):data_dic['Diff']['cont_range']={} #Deactivate conditions if gen_dic['DIbin']==False: data_dic['DI']['fit_MCCFout']=False if (not gen_dic['specINtype']) or (gen_dic['DI_CCF']):plot_dic['spectral_LC']='' if (not gen_dic['diff_data']): for key in ['map_Diff_prof','Diff_prof']:plot_dic[key]='' if not gen_dic['pl_atm']: for key in ['map_Atm_prof','sp_atm','CCFatm']:plot_dic[key]='' else: if gen_dic['Intr_CCF']:stop('Error: atmospheric extraction cannot be performed after Diff./Intr. CCF conversion') #------------------------------------------------------------------------------------------------------------------------ #Weighing conditions # - automatic activation to calculate or use contributors to variance profiles of disk-integrated spectra, based on the use of modules that average profiles requesting this variance # + extraction of Diff profiles, requiring DI master (if intrinsic or atmospheric extraction is required, the extraction of Diff profiles is as well) # + estimate of Intr profiles based on DI master # + conversion of profiles from 2D into 1D # + binning profiles together # - we disable calculation of DI weighing master if it is requested but no other module requiring this weight are recalculated weights_on = '' weights_off = '' cond_def = (gen_dic['diff_data'] | (gen_dic['loc_prof_est'] & (data_dic['Intr']['opt_loc_prof_est']['corr_mode'] in ['DIbin','Intrbin'])) | gen_dic['spec_1D'] | gen_dic['bin'] | gen_dic['binmultivis']) gen_dic['DImast_weight'] = cond_def if gen_dic['DImast_weight']:weights_on+='(DI master) ' else:weights_off+='(DI master) ' if gen_dic['DImast_weight'] and gen_dic['calc_DImast']:gen_dic['calc_DImast'] = gen_dic['calc_diff_data'] | (gen_dic['calc_loc_prof_est'] & (data_dic['Intr']['opt_loc_prof_est']['corr_mode'] in ['DIbin','Intrbin'])) | gen_dic['calc_spec_1D'] | gen_dic['calc_bin'] | gen_dic['calc_binmultivis'] if ('spec2D' in gen_dic['all_types']):gen_dic['cal_weight'] = cond_def | gen_dic['corr_wig'] else:gen_dic['cal_weight'] = False if gen_dic['cal_weight']:weights_on+='(Calibration) ' else:weights_off+='(Calibration) ' if gen_dic['specINtype']:gen_dic['tell_weight'] = cond_def & gen_dic['corr_tell'] else:gen_dic['tell_weight'] = False if gen_dic['tell_weight']:weights_on+='(Tellurics) ' else:weights_off+='(Tellurics) ' print('Weight contributors :') if weights_on != '':print(' On :'+weights_on) else:print(' On : None') if weights_off != '':print(' Off :'+weights_off) else:print(' Off : None') #------------------------------------------------------------------------------------------------------------------------ #Automatic activation of light curve computation # - required to calculate binned disk-integrated profiles or convert them from 2D to 1D, to convert Diff profiles into CCFs or from 2D to 1D, and to compute (and process) intrinsic and atmospheric profiles if (not gen_dic['flux_sc']) and (gen_dic['DImast_weight'] or gen_dic['spec_1D_DI'] or gen_dic['DIbin'] or gen_dic['Diff_CCF'] or gen_dic['spec_1D_Diff'] or gen_dic['intr_data'] or (gen_dic['pl_atm'] and data_dic['Atm']['pl_atm_sign']=='Absorption')): print('Automatic activation of flux scaling calculation') gen_dic['flux_sc'] = True #Set general conditions to activate joined datasets modules gen_dic['joined_ana']=False gen_dic['fit_DIProf'] = False #module undefined for now gen_dic['fit_DiffProp'] = False #module undefined for now for key in gen_dic['type_list']:gen_dic['joined_ana'] |= (gen_dic['fit_'+key+'Prof'] or gen_dic['fit_'+key+'Prop']) #Standard-deviation curves with bin size for the out-of-transit differential CCFs # - defining the maximum size of the binning window, and the binning size for the sliding window (we will average the bins in windows of width bin_size from 1 to 40 (arbitrary)) if gen_dic['scr_search']: gen_dic['scr_srch_max_binwin']=40. gen_dic['scr_srch_nperbins'] = 1.+np.arange(gen_dic['scr_srch_max_binwin']) #Default options for instrumental calibration if gen_dic['gcal']: for inst in data_dic['instrum_list'] : #Spectral bin width (A) if (inst not in gen_dic['gcal_binw']): def_gcal_binw = {'CORALIE':1.,'ESPRESSO':0.15,'HARPS':0.25,'NIGHT':0.3,'NIRPS_HE':0.3} if inst in def_gcal_binw:gen_dic['gcal_binw'][inst] = def_gcal_binw[inst] else:stop('ERROR: no default value of "gen_dic["gcal_binw"]" for '+inst+'. Define one, or run the module with a custom value for this field.') #Check for phantom bins # - bin size of calibration profiles must be small enough to sample variations at the edges dw_edge = {'HARPS':1.7,'ESPRESSO':1.} if inst in dw_edge:dw_edge_inst = dw_edge[inst] else:stop('ERROR: no default value of "dw_edge" for '+inst+'. Define one in ANTARESS_main.py > init_gen().') gen_dic['gcal_n_edge'] = int(dw_edge_inst/gen_dic['gcal_binw'][inst]) if gen_dic['gcal_n_edge']<5:stop('ERROR : too few edge points ('+str(gen_dic['gcal_n_edge'])+') over '+str(dw_edge_inst)+' A : decrease gen_dic["gcal_binw"] ('+str(gen_dic['gcal_binw'][inst])+') for '+inst) #Spline smoothing factor if (inst not in gen_dic['gcal_smooth']): def_gcal_smooth = {'HARPS':1e-3,'ESPRESSO':1e-3} if inst in def_gcal_smooth:gen_dic['gcal_smooth'][inst] = def_gcal_smooth[inst] else:stop('ERROR: no default value of "gen_dic["gcal_smooth"]" for '+inst+'. Define one, or run the module with a custom value for this field.') #Threshold if (inst not in gen_dic['gcal_thresh']):gen_dic['gcal_thresh'][inst] = {'outliers':5.,'global':1e10} #Stellar continuum for key in ['DI','Intr']: if (gen_dic[key+'_stcont'] and (not (gen_dic[key+'bin'] or gen_dic[key+'binmultivis']))): print("WARNING: binned "+key+" spectrum must be calculated to use gen_dic['"+key+"_stcont']") #Default options for stellar continuum calculation if gen_dic['mask_permpeak'] or gen_dic['DI_stcont'] or gen_dic['Intr_stcont']: for inst in data_dic['instrum_list'] : #Size of rolling window for peak exclusion # - in A if (inst not in gen_dic['contin_roll_win']):gen_dic['contin_roll_win'][inst] = 2. #Size of smoothing window # - in A if (inst not in gen_dic['contin_smooth_win']):gen_dic['contin_smooth_win'][inst] = 0.5 #Size of local maxima window # - in A if (inst not in gen_dic['contin_locmax_win']):gen_dic['contin_locmax_win'][inst] = 0.5 #Flux/wavelength stretching if (inst not in gen_dic['contin_stretch']):gen_dic['contin_stretch'][inst] = 10. #Rolling pin radius # - value corresponds to the bluest wavelength of the processed spectra if (inst not in gen_dic['contin_pinR']):gen_dic['contin_pinR'][inst] = 5. #Raise warning if gen_dic['EvE_outputs']: if not gen_dic['DIbin']:print('WARNING (EvE outputs): a master-out spectrum must be calculated (gen_dic["DIbin"] = True)') #------------------------------------------------------------------------------------------------------------------------ #Stellar mask generation is requested gen_dic['def_st_masks'] = gen_dic['def_DImasks'] | gen_dic['def_Intrmasks'] #Mask used to compute CCF on stellar lines if len(gen_dic['CCF_mask'])>0: gen_dic['CCF_mask_wav']={} gen_dic['CCF_mask_wgt']={} #Condition to exclude ranges contaminated by planetary absorption # - independent of atmospheric signal extraction if len(data_dic['Atm']['no_plrange'])>0: data_dic['Atm']['exc_plrange']=True data_dic['Atm']['plrange']=np.array(data_dic['Atm']['plrange']) else:data_dic['Atm']['exc_plrange']=False #Properties associated with planetary signal extraction if gen_dic['calc_pl_atm']: #Deactivate signal reduction for input CCFs if (not gen_dic['specINtype']):gen_dic['atm_CCF']=False #Mask used to compute atmospheric CCFs and/or exclude atmospheric planetary ranges if data_dic['Atm']['CCF_mask'] is not None: #Upload CCF mask ext = data_dic['Atm']['CCF_mask'].split('.')[-1] if (ext=='fits'): hdulist = fits.open(data_dic['Atm']['CCF_mask']) data_loc = hdulist[1].data data_dic['Atm']['CCF_mask_wav'] = data_loc['lambda'] if data_dic['Atm']['use_maskW']:data_dic['Atm']['CCF_mask_wgt'] = data_loc['contrast'] hdulist.close() elif (ext=='csv'): data_loc = np.genfromtxt(data_dic['Atm']['CCF_mask'], delimiter=',', names=True) data_dic['Atm']['CCF_mask_wav'] = data_loc['wave'] if data_dic['Atm']['use_maskW']:data_dic['Atm']['CCF_mask_wgt'] = data_loc['contrast'] elif (ext in ['txt','dat']): data_loc = np.loadtxt(data_dic['Atm']['CCF_mask']).T data_dic['Atm']['CCF_mask_wav'] = data_loc[0] if data_dic['Atm']['use_maskW']:data_dic['Atm']['CCF_mask_wgt'] = data_loc[1] else: stop('CCF mask extension TBD') #No weighing of lines if not data_dic['Atm']['use_maskW']:data_dic['Atm']['CCF_mask_wgt'] = np.repeat(1.,len(data_dic['Atm']['CCF_mask_wav'])) #------------------------------------------------------------------------------ #Star #------------------------------------------------------------------------------ if gen_dic['sequence'] not in ['st_master_tseries']: #Conversions and calculations system_param['star']['Rstar_km'] = system_param['star']['Rstar']*Rsun #Spherical star if ('f_GD' not in system_param['star']): system_param['star']['f_GD']=0. system_param['star']['RpoleReq'] = 1. #Oblate star elif system_param['star']['f_GD']>0.: print('Star is oblate') system_param['star']['RpoleReq']=1.-system_param['star']['f_GD'] #Reference time for stellar phase (bjd) if 'Tcenter' not in system_param['star']:system_param['star']['Tcenter'] = 2400000. #Stellar equatorial rotation rate (rad/s) # - om = 2*pi/P = v/R system_param['star']['om_eq'] = system_param['star']['veq']/system_param['star']['Rstar_km'] for key in ['_spots','_faculae']: if 'veq'+key in system_param['star']:system_param['star']['om_eq'+key]=system_param['star']['veq'+key]/system_param['star']['Rstar_km'] else: system_param['star']['veq'+key]=system_param['star']['veq'] system_param['star']['om_eq'+key]=system_param['star']['om_eq'] #Stellar equatorial rotation period (d) # - P = 2*pi*R/v for key in ['','_spots','_faculae']: system_param['star']['Peq'+key] = (2.*np.pi*system_param['star']['Rstar_km'])/(system_param['star']['veq'+key]*24*3600) #No GD if ('beta_GD' not in system_param['star']):system_param['star']['beta_GD']=0. if ('Tpole' not in system_param['star']):system_param['star']['Tpole']=0. #Conversions system_param['star']['istar_rad']=system_param['star']['istar']*np.pi/180. system_param['star']['cos_istar']=np.cos(system_param['star']['istar_rad']) system_param['star']['vsini']=system_param['star']['veq']*np.sin(system_param['star']['istar_rad']) #km/s #Default parameters for key in ['alpha_rot','beta_rot','alpha_rot_spots','beta_rot_spots','alpha_rot_faculae','beta_rot_faculae','c1_CB','c2_CB','c3_CB','c1_pol','c2_pol','c3_pol','c4_pol']: if key not in system_param['star']:system_param['star'][key] = 0. #Conversion factor from the LOS velocity output (/Rstar/h) to RV in km/s system_param['star']['RV_conv_fact']=-system_param['star']['Rstar_km']/3600. #Macroturbulence if theo_dic['mac_mode'] is not None: if 'rt' in theo_dic['mac_mode']: theo_dic['mac_mode_func'] = calc_macro_ker_rt if theo_dic['mac_mode']=='rt_iso': system_param['star']['A_T']=system_param['star']['A_R'] system_param['star']['ksi_T']=system_param['star']['ksi_R'] elif 'anigauss' in theo_dic['mac_mode']: theo_dic['mac_mode_func'] = calc_macro_ker_anigauss if theo_dic['mac_mode']=='anigauss_iso': system_param['star']['eta_T']=system_param['star']['eta_R'] #------------------------------------------------------------------------------ #Planets #------------------------------------------------------------------------------ #Planets defined in the system gen_dic['all_pl'] = [pl_loc for pl_loc in system_param.keys() if pl_loc!='star'] #Planets considered for analysis gen_dic['studied_pl_list'] = list(gen_dic['studied_pl'].keys()) if len(gen_dic['studied_pl_list'])>0: txt_print = 'Study of: '+gen_dic['studied_pl_list'][0] if len(gen_dic['studied_pl_list'])>1: for pl_loc in gen_dic['studied_pl_list'][1::]:txt_print+=', '+pl_loc print(txt_print) #Keplerian motion if ('kepl_pl' not in gen_dic):gen_dic['kepl_pl']=['all'] if (gen_dic['kepl_pl']==['all']): print('Accounting for Keplerian motion from all planets') gen_dic['kepl_pl']=deepcopy(gen_dic['all_pl']) #Planet properties gen_dic['tr_pl']=[] for pl_loc in gen_dic['all_pl']: print('Initializing planet '+pl_loc) #Checking if there is a "-" in a target name if '-' in pl_loc:stop('ERROR : Invalid target name: {}. Target names should not contain a hyphen.'.format(pl_loc)) PlParam_loc=system_param[pl_loc] #Converting argument of periastron from degree to radians PlParam_loc['omega_rad']=PlParam_loc['omega_deg']*np.pi/180. #Converting orbital period from days to seconds PlParam_loc['period_s'] = PlParam_loc['period']*24.*3600. #Planet with known orbit and transit properties if pl_loc in gen_dic['def_pl']: #Adding absolute semi-major axis if ('a' not in PlParam_loc):PlParam_loc['a']=PlParam_loc['aRs']*system_param['star']['Rstar_km']/AU_1 # in au #Converting orbital inclination from degree to radians PlParam_loc['inclin_rad']=PlParam_loc['inclination']*np.pi/180. #Deriving planet mass (in Mj) if ('Msini' in PlParam_loc) and ('Mp' not in PlParam_loc):PlParam_loc['Mp']=PlParam_loc['Msini']/np.sin(PlParam_loc['inclin_rad']) #Semi-amplitude of planet orbital motion around the star (approximated with Mp << Mstar) in km/s PlParam_loc['Kp_orb'] = (2.*np.pi/PlParam_loc['period_s'])*np.sin(PlParam_loc['inclin_rad'])*PlParam_loc['a']*AU_1/np.sqrt(1.-PlParam_loc['ecc']**2.) #Stellar mass derived from semi-orbital motion around the star # - assuming Mp << Ms : # P^2/a^3 = 4*pi^2/G*Ms # Ms = 4*pi^2*a^3/(G*P^2) PlParam_loc['Mstar_orb'] = 4.*np.pi**2.*(PlParam_loc['a']*AU_1*1e3)**3./(Msun*G_usi*PlParam_loc['period_s']**2.) #Converting sky-projected angle from degree to radians PlParam_loc['lambda_rad']=PlParam_loc['lambda_proj']*np.pi/180. #Impact parameter if 'b' not in PlParam_loc:PlParam_loc['b']=PlParam_loc['aRs']*np.cos(PlParam_loc['inclin_rad']) if (PlParam_loc['ecc']==0.) & (PlParam_loc['b']>1.):print('WARNING: impact parameter for ',pl_loc,' > 1') #Calculation of transit center (or time of cunjunction for a non-transiting planet) # - when T0 is not given but Tperiastron is known # - we could calculate directly the true anomaly from Tperi, but all phases in the routine # are defined relatively with the transit center # - mean anomaly is counted from pericenter as M(t) = 2pi*(t - Tperi)/P # M(Tcenter) = 2pi*(Tcenter - Tperi)/P # Tcenter = Tperi + M(Tcenter)*P/2pi if ('TCenter' not in PlParam_loc): PlParam_loc['TCenter']=PlParam_loc['Tperi'] if (PlParam_loc['ecc']>0): PlParam_loc['Mean_anom_TR'] = calc_mean_anom_TR(PlParam_loc['ecc'],PlParam_loc['omega_rad']) PlParam_loc['TCenter']+=(PlParam_loc['Mean_anom_TR']*PlParam_loc["period"]/(2.*np.pi)) #Keplerian semi-amplitude of the star from the planet (km/s) PlParam_loc['Kstar_kms']=PlParam_loc['Kstar']/1000. if 'Kstar' in PlParam_loc else calc_Kstar(PlParam_loc,system_param['star'])/1000. #Orbital frequency, in year-1 PlParam_loc['omega_p']=2.*np.pi*365.2425/PlParam_loc['period'] #Transit check # - this call is also used to check for transit status if (pl_loc in gen_dic['studied_pl_list']): if (pl_loc not in data_dic['DI']['system_prop']['achrom']): print(' No RpRs defined for '+pl_loc+' : assuming object is not transiting') elif ('inclination' not in PlParam_loc): print(' No inclination defined for '+pl_loc+' : assuming object is not transiting') else: #Contact phases contact_phases=calc_tr_contacts(data_dic['DI']['system_prop']['achrom'][pl_loc][0],PlParam_loc,plot_dic['stend_ph'],system_param['star']) #Planet is transiting if contact_phases is not None: print(' Transiting object') gen_dic['tr_pl']+=[pl_loc] #Transit duration if ('TLength' not in PlParam_loc): PlParam_loc['TLength'] = (contact_phases[3]-contact_phases[0])*PlParam_loc['period'] print(' Automatic definition of T14['+str(pl_loc)+']='+"{0:.2f}".format(PlParam_loc['TLength']*24.)+' h') else: print(' Non-transiting object') #Calculating theoretical properties of the planet-occulted regions if (len(gen_dic['tr_pl'])>0) and (gen_dic['sequence'] not in ['st_master_tseries']) and (not gen_dic['theoPlOcc']) and ((('CCF' in data_dic['Diff']['type'].values()) and (gen_dic['fit_Intr'])) or (gen_dic['align_Intr']) or (gen_dic['calc_pl_atm'])): print('Automatic activation of "gen_dic["theoPlOcc"]"') gen_dic['theoPlOcc']=True #Oversampling factor for values from planet-occulted regions # - uses the nominal planet-to-star radius ratios, which must correspond to the band from which local properties are derived theo_dic['d_oversamp_pl']={} for pl_loc in theo_dic['n_oversamp']: if (pl_loc in gen_dic['tr_pl']) and (theo_dic['n_oversamp'][pl_loc]>0.): theo_dic['d_oversamp_pl'][pl_loc] = data_dic['DI']['system_prop']['achrom'][pl_loc][0]/theo_dic['n_oversamp'][pl_loc] #Set flag for errors on estimates for local stellar profiles (depending on whether they are derived from data or models) if data_dic['Intr']['opt_loc_prof_est']['corr_mode'] in ['DIbin','Intrbin']:data_dic['Intr']['cov_loc_star']=True else:data_dic['Intr']['cov_loc_star']=False #Transit and stellar surfce chromatic properties # - must be defined for data processing, transit scaling, calculation of planet properties # - we remove the 'chrom' dictionary if no spectral data is used as input or if a single band is defined if ('LD' not in data_dic['DI']['system_prop']['achrom']): data_dic['DI']['system_prop']['achrom']['LD']=['uniform'] print('WARNING : limb-darkening law set to "uniform". Define it with "data_dic["DI"]["system_prop"]["achrom"]"') for ideg in range(2,5): if 'LD_u'+str(ideg) not in data_dic['DI']['system_prop']['achrom']:data_dic['DI']['system_prop']['achrom']['LD_u'+str(ideg)] = [0.] if ('GD_dw' in data_dic['DI']['system_prop']['achrom']): system_param['star']['GD']=True print('Star is gravity-darkened') else:system_param['star']['GD']=False data_dic['DI']['system_prop']['achrom']['w']=[None] data_dic['DI']['system_prop']['achrom']['nw']=1 data_dic['DI']['system_prop']['achrom']['cond_in_RpRs']={} data_dic['DI']['system_prop']['chrom_mode'] = 'achrom' if ('chrom' in data_dic['DI']['system_prop']): if (not gen_dic['specINtype']) or (len(data_dic['DI']['system_prop']['chrom']['w'])==1):data_dic['DI']['system_prop'].pop('chrom') else: data_dic['DI']['system_prop']['chrom_mode'] = 'chrom' data_dic['DI']['system_prop']['chrom']['w'] = np.array(data_dic['DI']['system_prop']['chrom']['w']) data_dic['DI']['system_prop']['chrom']['nw']=len(data_dic['DI']['system_prop']['chrom']['w']) data_dic['DI']['system_prop']['chrom']['cond_in_RpRs']={} #Typical scale of chromatic variations w_edge = def_edge_tab(data_dic['DI']['system_prop']['chrom']['w'][None,:][None,:])[0,0] data_dic['DI']['system_prop']['chrom']['dw'] = w_edge[1::]-w_edge[0:-1] data_dic['DI']['system_prop']['chrom']['med_dw'] = np.median(data_dic['DI']['system_prop']['chrom']['dw']) #Store properties at the stage of broadband scaling data_dic['DI']['system_prop_sc'] = deepcopy(data_dic['DI']['system_prop']) if gen_dic['DI_CCF'] and ('chrom' in data_dic['DI']['system_prop']): data_dic['DI']['system_prop_sc']['chrom_mode'] = 'achrom' data_dic['DI']['system_prop_sc'].pop('chrom') #Default transit model if 'transit_prop' not in data_dic['DI']:data_dic['DI']['transit_prop']={} #Definition of grids discretizing planets disk to calculate planet-occulted properties theo_dic['x_st_sky_grid_pl']={} theo_dic['y_st_sky_grid_pl']={} theo_dic['Ssub_Sstar_pl']={} data_dic['DI']['system_prop']['RpRs_max']={} if ('nsub_Dpl' not in theo_dic):theo_dic['nsub_Dpl']={} for pl_loc in gen_dic['tr_pl']: #Largest possible planet size if ('chrom' in data_dic['DI']['system_prop']):data_dic['DI']['system_prop']['RpRs_max'][pl_loc] = np.max(data_dic['DI']['system_prop']['achrom'][pl_loc]+data_dic['DI']['system_prop']['chrom'][pl_loc]) else:data_dic['DI']['system_prop']['RpRs_max'][pl_loc] = data_dic['DI']['system_prop']['achrom'][pl_loc][0] #Default grid scaled from a 51x51 grid for a hot Jupiter transiting a solar-size star # - dsub_ref = (2*Rjup/Rsun)*(1/51) # nsub_Dpl = int(2*RpRs/dsub_ref) = int( 51*RpRs*Rsun/Rjup ) if (pl_loc not in theo_dic['nsub_Dpl']): theo_dic['nsub_Dpl'][pl_loc] =int( 51.*data_dic['DI']['system_prop']['RpRs_max'][pl_loc]*Rsun/Rjup ) print('Default nsub_Dpl['+str(pl_loc)+']='+str(theo_dic['nsub_Dpl'][pl_loc])) #Corresponding planet grid _,theo_dic['Ssub_Sstar_pl'][pl_loc],theo_dic['x_st_sky_grid_pl'][pl_loc],theo_dic['y_st_sky_grid_pl'][pl_loc],r_sub_pl2=occ_region_grid(data_dic['DI']['system_prop']['RpRs_max'][pl_loc],theo_dic['nsub_Dpl'][pl_loc]) #Identification of cells within the nominal and chromatic planet radii data_dic['DI']['system_prop']['achrom']['cond_in_RpRs'][pl_loc] = [(r_sub_pl2<data_dic['DI']['system_prop']['achrom'][pl_loc][0]**2.)] if ('chrom' in data_dic['DI']['system_prop']): data_dic['DI']['system_prop']['chrom']['cond_in_RpRs'][pl_loc]={} for iband in range(data_dic['DI']['system_prop']['chrom']['nw']): data_dic['DI']['system_prop']['chrom']['cond_in_RpRs'][pl_loc][iband] = (r_sub_pl2<data_dic['DI']['system_prop']['chrom'][pl_loc][iband]**2.) #------------------------------------------------------------------------------ #Active regions #------------------------------------------------------------------------------ #Generating multiple active regions if (mock_dic['auto_gen_ar'] != {}):generate_ar_prop(mock_dic, data_dic, gen_dic) #Active regions considered for analysis gen_dic['studied_ar_list'] = list(gen_dic['studied_ar'].keys()) if len(gen_dic['studied_ar_list'])!=len(np.unique(gen_dic['studied_ar_list'])):stop('Active regions must have unique names in each visit') if len(gen_dic['studied_ar_list'])>0: txt_print = 'Study of: '+gen_dic['studied_ar_list'][0] if len(gen_dic['studied_ar_list'])>1: for pl_loc in gen_dic['studied_ar_list'][1::]:txt_print+=', '+pl_loc print(txt_print) #Initialize active region use theo_dic['x_st_sky_grid_ar']={} theo_dic['y_st_sky_grid_ar']={} theo_dic['Ssub_Sstar_ar'] = {} theo_dic['d_oversamp_ar']={} #Active region activation triggered gen_dic['ar_coord_par'] = ['lat_rad_exp','sin_lat_exp','cos_lat_exp','long_rad_exp','sin_long_exp','cos_long_exp','x_sky_exp','y_sky_exp','z_sky_exp'] if len(gen_dic['studied_ar_list'])>0: print('Active regions are simulated') #Oversampling factor for active region-occulted regions # - input corresponds to the half-angular opening of the active region # taking the largest projected radius of the active region if placed along the LOS: # sin(ang) = Rproj/Rstar for ar in theo_dic['n_oversamp_ar']: if (theo_dic['n_oversamp_ar'][ar]>0.): theo_dic['d_oversamp_ar'][ar] = np.sin(data_dic['DI']['ar_prop']['achrom'][ar][0])/theo_dic['n_oversamp_ar'][ar] #Active region surface chromatic properties for ideg in range(2,5): if 'LD_u'+str(ideg) not in data_dic['DI']['ar_prop']['achrom']:data_dic['DI']['ar_prop']['achrom']['LD_u'+str(ideg)] = [0.] #Need to define chromatic band properties data_dic['DI']['ar_prop']['achrom']['w']=[None] data_dic['DI']['ar_prop']['achrom']['nw']=1 data_dic['DI']['ar_prop']['chrom_mode'] = 'achrom' if ('chrom' in data_dic['DI']['ar_prop']): if (not gen_dic['specINtype']) or (len(data_dic['DI']['ar_prop']['chrom']['w'])==1):data_dic['DI']['ar_prop'].pop('chrom') else: data_dic['DI']['ar_prop']['chrom_mode'] = 'chrom' data_dic['DI']['ar_prop']['chrom']['w'] = np.array(data_dic['DI']['ar_prop']['chrom']['w']) data_dic['DI']['ar_prop']['chrom']['nw']=len(data_dic['DI']['ar_prop']['chrom']['w']) #Typical scale of chromatic variations w_edge = def_edge_tab(data_dic['DI']['ar_prop']['chrom']['w'][None,:][None,:])[0,0] data_dic['DI']['ar_prop']['chrom']['dw'] = w_edge[1::]-w_edge[0:-1] data_dic['DI']['ar_prop']['chrom']['med_dw'] = np.median(data_dic['DI']['ar_prop']['chrom']['dw']) #Definition of grids discretizing active region disk to calculate active region properties for ar in gen_dic['studied_ar_list']: #Maximum projected active region radius RarRs_max = np.sin(data_dic['DI']['ar_prop']['achrom'][ar][0]) #Default grid scaled from a 51x51 grid for an active region with projected radius equal to a hot Jupiter transiting a solar-size star # - dsub_ref = (2*Rjup/Rsun)*(1/51) # nsub_Dar = int(2*RpRs/dsub_ref) = int( 51*RpRs*Rsun/Rjup ) = int( 51*sin(ang)*Rsun/Rjup ) if (ar not in theo_dic['nsub_Dar']): theo_dic['nsub_Dar'][pl_loc] =int( 51.*RarRs_max*Rsun/Rjup ) print('Default nsub_Dar['+str(ar)+']='+str(theo_dic['nsub_Dar'][ar])) #Default oversampling if (ar not in theo_dic['n_oversamp_ar']): if (theo_dic['n_oversamp_ar'] != {}):theo_dic['n_oversamp_ar'][ar] = theo_dic['n_oversamp_ar'][list(theo_dic['n_oversamp_ar'].keys())[0]] elif (theo_dic['n_oversamp'] != {}):theo_dic['n_oversamp_ar'][ar] = theo_dic['n_oversamp'][list(theo_dic['n_oversamp'].keys())[0]] else:theo_dic['n_oversamp_ar'][ar] = 1. #Corresponding active region grid _,theo_dic['Ssub_Sstar_ar'][ar],theo_dic['x_st_sky_grid_ar'][ar], theo_dic['y_st_sky_grid_ar'][ar],_ = occ_region_grid(RarRs_max, theo_dic['nsub_Dar'][ar],planet=False) #------------------------------------------------------------------------------ #Generic path names gen_dic['main_pl_text'] = '' gen_dic['main_pl_ar_text'] = '' for pl_loc in gen_dic['studied_pl_list']:gen_dic['main_pl_text']+=pl_loc save_system = gen_dic['save_dir']+gen_dic['star_name']+gen_dic['dir_suff']+'/' if gen_dic['main_pl_text'] != '':save_system+=gen_dic['main_pl_text']+'_' gen_dic['save_data_dir'] = save_system+'Saved_data/' gen_dic['save_plot_dir'] = save_system+'Plots/' gen_dic['add_txt_path']={'DI':'','Intr':'','Diff':'','Atm':data_dic['Atm']['pl_atm_sign']+'/'} gen_dic['data_type_gen']={'DI':'DI','Diff':'Diff','Intr':'Intr','Absorption':'Atm','Emission':'Atm'} gen_dic['typegen2type']={'DI':'DI','Diff':'Diff','Intr':'Intr','Atm':data_dic['Atm']['pl_atm_sign']} #Data is processed if gen_dic['sequence'] not in ['system_view']: #------------------------------------------------------------------------------------------------------------------------ #Model star #------------------------------------------------------------------------------------------------------------------------ grid_type=[] if any('spec' in s for s in data_dic['DI']['type'].values()):grid_type+=['spec'] if ('CCF' in data_dic['DI']['type'].values()):grid_type+=['ccf'] #Calculation of total stellar flux for use in simulated light curves if gen_dic['calc_flux_sc'] and (data_dic['DI']['transit_prop']['nsub_Dstar'] is not None): model_star('Ftot',theo_dic,grid_type,data_dic['DI']['system_prop'],data_dic['DI']['transit_prop']['nsub_Dstar'],system_param['star'],True,True) #Definition of model stellar grid to calculate local or disk-integrated properties # - used throughout the pipeline, unless stellar properties are fitted if gen_dic['theoPlOcc'] or (data_dic['DI']['ar_prop'] != {}) or (gen_dic['fit_DI_gen'] and (('custom' in data_dic['DI']['model'].values()) or ('RT_ani_macro' in data_dic['DI']['model'].values()))) or gen_dic['mock_data'] \ or gen_dic['fit_DiffProf'] or gen_dic['fit_IntrProf'] or gen_dic['loc_prof_est']: #or gen_dic['correct_ar'] #Stellar grid model_star('grid',theo_dic,grid_type,data_dic['DI']['system_prop'],theo_dic['nsub_Dstar'],system_param['star'],True,True) #Theoretical atmosphere cond_st_atm = False if gen_dic['mock_data']: for inst in mock_dic['intr_prof']: if (mock_dic['intr_prof'][inst]['mode']=='theo' ):cond_st_atm = True if (gen_dic['fit_DI_gen'] and ('custom' in data_dic['DI']['model'].values())): for inst in data_dic['DI']['mod_def']: if (data_dic['DI']['mod_def'][inst]['mode'] == 'theo'):cond_st_atm = True if gen_dic['fit_IntrProf'] and (glob_fit_dic['IntrProf']['mode'] == 'theo'):cond_st_atm = True if gen_dic['loc_prof_est'] and (data_dic['Intr']['opt_loc_prof_est']['corr_mode'] in ['glob_mod','indiv_mod']) and (data_dic['Intr']['opt_loc_prof_est']['mode']=='theo'):cond_st_atm = True if cond_st_atm: if theo_dic['st_atm']['calc']: theo_dic['sme_grid'] = gen_theo_atm(theo_dic['st_atm'],system_param['star']) datasave_npz(gen_dic['save_data_dir']+'Introrig_prop/IntrProf_grid',{'sme_grid':theo_dic['sme_grid']}) else:theo_dic['sme_grid'] = dataload_npz(gen_dic['save_data_dir']+'Introrig_prop/IntrProf_grid')['sme_grid'] #------------------------------------------------------------------------------------------------------------------------ #Define conservative maximum shift of spectral table (km/s) # - this is used to extend the common spectral table (oiginally defined in the input rest frame) so that it covers all possible wavelength shifts of the profiles in : # + the input rest frame, where profiles are shifted by up to the systemic rv + Keplerian semi-amplitude into the star rest frame # w_star ~ w_solbar * (1 - (rv_kep + rv_sys)/c) # + the star rest frame, where profiles are shifted by up to the sky-projected rotational velocity into the surface rest frame # w_surf ~ w_star * (1 - rv[surf/star]/c)) # + the star rest frame, where profiles are shifted by up to the orbital semi-amplitude into the planetary rest frame # w_pl ~ w_star * (1 - rv[pl/star]/c)) # - where we approximate the Doppler shifts to define their maximum values # - a negative rv is associated with a redshift of the table, so that we take the opposite of vsini and the semi-amplitude to evaluate the maximum positive extension gen_dic['min_rv_shift'] = 0. gen_dic['max_rv_shift'] = 0. for pl_loc in gen_dic['kepl_pl']: gen_dic['min_rv_shift']+=system_param[pl_loc]['Kstar_kms'] gen_dic['max_rv_shift']-=system_param[pl_loc]['Kstar_kms'] if gen_dic['intr_data']: gen_dic['min_rv_shift'] = np.max([gen_dic['min_rv_shift'], system_param['star']['vsini']]) gen_dic['max_rv_shift'] = np.min([gen_dic['max_rv_shift'],-system_param['star']['vsini']]) if gen_dic['pl_atm']: max_rv_shift_orb = 0. for pl_loc in gen_dic['studied_pl_list']: max_rv_shift_orb=system_param[pl_loc]['Kp_orb'] gen_dic['min_rv_shift'] = np.max([gen_dic['min_rv_shift'], max_rv_shift_orb]) gen_dic['max_rv_shift'] = np.min([gen_dic['max_rv_shift'],-max_rv_shift_orb]) gen_dic['min_rv_shift'] -= system_param['star']['sysvel'] gen_dic['min_rv_shift'] = np.max([gen_dic['min_rv_shift'],0.]) gen_dic['max_rv_shift'] += system_param['star']['sysvel'] gen_dic['max_rv_shift'] = np.min([gen_dic['max_rv_shift'],0.]) #------------------------------------------------------------------------------------------------------------------------ #Create directories where data is saved/restored if (not path_exist(gen_dic['save_data_dir']+'Processed_data/Global/')):makedirs(gen_dic['save_data_dir']+'Processed_data/Global/') if gen_dic['specINtype']: dir_name_dic={'gcal':'Processed_data/Calibration','CCF_from_sp':'Processed_data/CCFfromSpec'} for key in dir_name_dic.keys(): if (gen_dic[key]) and (not path_exist(gen_dic['save_data_dir']+dir_name_dic[key]+'/')):makedirs(gen_dic['save_data_dir']+dir_name_dic[key]+'/') if (gen_dic['corr_data']): corr_path = gen_dic['save_data_dir']+'Corr_data/' if (not path_exist(corr_path)):makedirs(corr_path) dir_name_dic={'corr_tell':'Tell','glob_mast':'Global_Master','corr_Fbal':'Fbal','corr_FbalOrd':'Fbal/Orders','corr_Ftemp':'Ftemp','corr_cosm':'Cosm','mask_permpeak':'Permpeak','corr_wig':'Wiggles'} for key in dir_name_dic.keys(): if (gen_dic[key]) and (not path_exist(corr_path+dir_name_dic[key]+'/')):makedirs(corr_path+dir_name_dic[key]+'/') if (gen_dic['corr_wig']): #Condition for exposure analysis gen_dic['wig_exp_ana'] = gen_dic['wig_exp_init']['mode'] | gen_dic['wig_exp_filt']['mode'] | gen_dic['wig_exp_samp']['mode'] | gen_dic['wig_exp_nu_ana']['mode'] | gen_dic['wig_exp_fit']['mode'] | gen_dic['wig_exp_point_ana']['mode'] if gen_dic['wig_exp_ana'] and (not path_exist(corr_path+'Wiggles/Exp_fit/')):makedirs(corr_path+'Wiggles/Exp_fit/') if gen_dic['wig_vis_fit']['mode'] and (not path_exist(corr_path+'Wiggles/Vis_fit/')):makedirs(corr_path+'Wiggles/Vis_fit/') if gen_dic['wig_corr']['mode'] and (not path_exist(corr_path+'Wiggles/Data/')):makedirs(corr_path+'Wiggles/Data/') if (gen_dic['corr_fring']) and (not path_exist(corr_path+'Fring/')):makedirs(corr_path+'Fring/') if (gen_dic['trim_spec']) and (not path_exist(corr_path+'Trim/')):makedirs(corr_path+'Trim/') if (gen_dic['detrend_prof']) and (not path_exist(gen_dic['save_data_dir']+'Detrend_prof/')):makedirs(gen_dic['save_data_dir']+'Detrend_prof/') if (gen_dic['flux_sc']) and (not path_exist(gen_dic['save_data_dir']+'Scaled_data/')):makedirs(gen_dic['save_data_dir']+'Scaled_data/') if gen_dic['DImast_weight'] and (not path_exist(gen_dic['save_data_dir']+'DI_data/Master/')):makedirs(gen_dic['save_data_dir']+'DI_data/Master/') if (gen_dic['diff_data']) and (not path_exist(gen_dic['save_data_dir']+'Diff_data/')):makedirs(gen_dic['save_data_dir']+'Diff_data/') if gen_dic['pca_ana'] and (not path_exist(gen_dic['save_data_dir']+'PCA_results/')):makedirs(gen_dic['save_data_dir']+'PCA_results/') if (gen_dic['intr_data']) and (not path_exist(gen_dic['save_data_dir']+'Intr_data/')):makedirs(gen_dic['save_data_dir']+'Intr_data/') if gen_dic['loc_prof_est']: if (not path_exist(gen_dic['save_data_dir']+'Loc_estimates/')):makedirs(gen_dic['save_data_dir']+'Loc_estimates/') if (not path_exist(gen_dic['save_data_dir']+'Loc_estimates/'+data_dic['Intr']['opt_diff_prof_est']['corr_mode']+'/')):makedirs(gen_dic['save_data_dir']+'Loc_estimates/'+data_dic['Intr']['opt_loc_prof_est']['corr_mode']+'/') if gen_dic['diff_prof_est']: if (not path_exist(gen_dic['save_data_dir']+'Diff_estimates/')):makedirs(gen_dic['save_data_dir']+'Diff_estimates/') if (not path_exist(gen_dic['save_data_dir']+'Diff_estimates/'+data_dic['Diff']['opt_diff_prof_est']['corr_mode']+'/')):makedirs(gen_dic['save_data_dir']+'Diff_estimates/'+data_dic['Diff']['opt_loc_prof_est']['corr_mode']+'/') if (gen_dic['pl_atm']): if (not path_exist(gen_dic['save_data_dir']+'Atm_data/')):makedirs(gen_dic['save_data_dir']+'Atm_data/') if (not path_exist(gen_dic['save_data_dir']+'Atm_data/'+data_dic['Atm']['pl_atm_sign']+'/')):makedirs(gen_dic['save_data_dir']+'Atm_data/'+data_dic['Atm']['pl_atm_sign']+'/') for data_type in ['DI','Intr','Atm']: if gen_dic['align_'+data_type] and (not path_exist(gen_dic['save_data_dir']+'Aligned_'+data_type+'_data/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+'Aligned_'+data_type+'_data/'+gen_dic['add_txt_path'][data_type]) if gen_dic['spec_1D_'+data_type]: if (not path_exist(gen_dic['save_data_dir']+data_type+'_data/1Dfrom2D/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+data_type+'_data/1Dfrom2D/'+gen_dic['add_txt_path'][data_type]) if (data_type=='Intr') and (not path_exist(gen_dic['save_data_dir']+'Diff_data/1Dfrom2D/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+'Diff_data/1Dfrom2D/'+gen_dic['add_txt_path'][data_type]) if (gen_dic[data_type+'bin'] or gen_dic[data_type+'binmultivis']) and (not path_exist(gen_dic['save_data_dir']+data_type+'bin_data/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+data_type+'bin_data/'+gen_dic['add_txt_path'][data_type]) if (gen_dic['fit_'+data_type+'bin'] or gen_dic['fit_'+data_type+'binmultivis']) and (not path_exist(gen_dic['save_data_dir']+data_type+'bin_prop/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+data_type+'bin_prop/'+gen_dic['add_txt_path'][data_type]) if gen_dic['def_'+data_type+'masks'] and (not path_exist(gen_dic['save_data_dir']+'CCF_masks_'+data_type+'/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+'CCF_masks_'+data_type+'/'+gen_dic['add_txt_path'][data_type]) if gen_dic[data_type+'_CCF']: if (not path_exist(gen_dic['save_data_dir']+data_type+'_data/CCFfromSpec/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+data_type+'_data/CCFfromSpec/'+gen_dic['add_txt_path'][data_type]) if (data_type=='Intr') and (not path_exist(gen_dic['save_data_dir']+'Diff_data/CCFfromSpec/'+gen_dic['add_txt_path'][data_type])):makedirs(gen_dic['save_data_dir']+'Diff_data/CCFfromSpec/'+gen_dic['add_txt_path'][data_type]) if (gen_dic['fit_DI'] or gen_dic['sav_keywords']) and (not path_exist(gen_dic['save_data_dir']+'DIorig_prop/')):makedirs(gen_dic['save_data_dir']+'DIorig_prop/') if ((gen_dic['fit_Intr']) or (gen_dic['theoPlOcc'])) and (not path_exist(gen_dic['save_data_dir']+'Introrig_prop/')):makedirs(gen_dic['save_data_dir']+'Introrig_prop/') if (gen_dic['fit_Atm']) and (not path_exist(gen_dic['save_data_dir']+'Atmorig_prop/'+data_dic['Atm']['pl_atm_sign']+'/')):makedirs(gen_dic['save_data_dir']+'Atmorig_prop/'+data_dic['Atm']['pl_atm_sign']+'/') for key in ['IntrProp','IntrProf','AtmProf','AtmProp']: if (gen_dic['fit_'+key]) and (not path_exist(gen_dic['save_data_dir']+'Joined_fits/'+key+'/')):makedirs(gen_dic['save_data_dir']+'Joined_fits/'+key+'/') return coord_dic,data_prop
[docs] def init_inst(mock_dic,inst,gen_dic,data_dic,theo_dic,data_prop,coord_dic,system_param,plot_dic): r"""**Initialization: instrument** Initializes instrument-specific fields for the workflow. Args: TBD Returns: None """ if (gen_dic['type'][inst]=='CCF'):inst_data='CCFs' elif (gen_dic['type'][inst]=='spec1D'):inst_data='1D spectra' elif (gen_dic['type'][inst]=='spec2D'):inst_data='2D echelle spectra' print(' Reading and initializing '+inst_data) #Initialize dictionaries for current instrument for key_dic in [gen_dic,coord_dic,data_dic['DI'],data_dic,data_dic['Diff'],data_dic['Intr'],data_dic['Atm'],data_prop,theo_dic]:key_dic[inst]={} DI_data_inst = data_dic['DI'][inst] if (inst not in gen_dic['scr_lgth']):gen_dic['scr_lgth'][inst]={} #Facility and reduction id facil_inst_all = { 'HARPS':'ESO', 'HARPN':'TNG', 'CARMENES_VIS':'CAHA','CARMENES_VIS_CCF':'CAHA', 'CORALIE':'ESO', 'ESPRESSO':'ESO','ESPRESSO_MR':'ESO', 'EXPRES':'DCT', 'MIKE_Red':'LCO','MIKE_Blue':'LCO', 'IGRINS2_Blue':'Gemini-N','IGRINS2_Red':'Gemini-N', 'MAROONX_Blue':'Gemini-N','MAROONX_Red':'Gemini-N', 'NIGHT':'ESO', #UPDATE TO 'OHP' WHEN FINAL FORMAT AVAILABLE 'NIRPS_HE':'ESO','NIRPS_HA':'ESO', 'SOPHIE':'OHP', } if inst not in facil_inst_all:stop('ERROR : define facility ID for '+inst+' in "facil_inst_all"') facil_inst = facil_inst_all[inst] #Error definition if not gen_dic['mock_data']: gen_dic['flag_err_inst'][inst] = flag_err_inst(inst) if (not gen_dic['flag_err_inst'][inst]) and gen_dic['gcal']:stop('ERROR : error table must be available to estimate calibration') if (inst in gen_dic['force_flag_err']):gen_dic['flag_err_inst'][inst]=False if gen_dic['flag_err_inst'][inst]:print(' > Errors propagated from raw data') else:print(' > Beware: custom definition of errors') else: if (inst in mock_dic['set_err']) and (mock_dic['set_err'][inst]):gen_dic['flag_err_inst'][inst] = True else:gen_dic['flag_err_inst'][inst] = False #Mask used to compute CCF on stellar lines if (gen_dic['CCF_from_sp']) and (inst in gen_dic['CCF_mask']): #Upload CCF mask ext = gen_dic['CCF_mask'][inst].split('.')[-1] if (ext=='fits'): hdulist = fits.open(gen_dic['CCF_mask'][inst]) data_loc = hdulist[1].data gen_dic['CCF_mask_wav'][inst] = data_loc['lambda'] if gen_dic['use_maskW']:gen_dic['CCF_mask_wgt'][inst] = data_loc['contrast'] hdulist.close() elif (ext=='csv'): data_loc = np.genfromtxt(gen_dic['CCF_mask'][inst], delimiter=',', names=True) gen_dic['CCF_mask_wav'][inst] = data_loc['wave'] if gen_dic['use_maskW']:gen_dic['CCF_mask_wgt'][inst] = data_loc['contrast'] elif (ext in ['txt','dat']): data_loc = np.loadtxt(gen_dic['CCF_mask'][inst]).T gen_dic['CCF_mask_wav'][inst] = data_loc[0] if gen_dic['use_maskW']:gen_dic['CCF_mask_wgt'][inst] = data_loc[1] else: stop('ERROR : CCF mask extension TBD') #No weighing of lines if not gen_dic['use_maskW']:gen_dic['CCF_mask_wgt'][inst] = np.repeat(1.,len(gen_dic['CCF_mask_wav'][inst])) #Same convention as the ESPRESSO DRS, which takes the 'contrast' column from mask files and use its squares as weights # - masks weights should have been normalized so that mean( weights ) = 1, where weight = CCF_mask_wgt^2 gen_dic['CCF_mask_wgt'][inst]=gen_dic['CCF_mask_wgt'][inst]**2. #Resampling exposures of a given visit on a common table or not # - CCFs will be resampled on a common RV table, as we consider they should only be used for preliminary analysis if gen_dic['type'][inst]=='CCF':gen_dic['comm_sp_tab'][inst]=True elif (inst not in gen_dic['comm_sp_tab']):gen_dic['comm_sp_tab'][inst]=False if (not gen_dic['comm_sp_tab'][inst]):print(' > Data processed on individual spectral tables for each exposure') else:print(' > Data resampled on a common spectral table') #Flux balance correction gen_dic['corr_FbalOrd_inst'][inst] = deepcopy(gen_dic['corr_FbalOrd']) if gen_dic['corr_FbalOrd_inst'][inst]: if (gen_dic['type'][inst]!='spec2D'): print('WARNING: input profiles for '+inst+' are not in 2D format, order flux balance correction will not be applied.') gen_dic['corr_FbalOrd_inst'][inst] = False elif (inst=='ESPRESSO'): print('WARNING: it is advised to disable order flux balance correction with ESPRESSO 2D spectra, as it may partly absorb wiggle signal.') #Deactivate conditions if gen_dic['pl_atm'] and (data_dic['Intr']['opt_loc_prof_est']['corr_mode'] in ['Intrbin','rec_prof']) and (data_dic['Atm']['type'][inst]!=data_dic['Intr']['type'][inst]): stop('Error: intrinsic profiles for '+inst+' (in '+data_dic['Intr']['type'][inst]+') must have the same format as planetary profiles (in '+data_dic['Atm']['type'][inst]+') if also requested for their extraction)') ############################################################################################################################## #Retrieval and pre-processing of data ############################################################################################################################## #Data is calculated and not retrieved if (gen_dic['calc_proc_data']): print(' Calculating data') #Initialize dictionaries for current instrument data_inst={'visit_list':[]} data_dic[inst]=data_inst #Initialize current data type for the instrument # - data need to be processed in their input mode at the start of each new instrument or visit # - if mode is switched to CCFs/s1d at some point in the reduction of a visit, the visit type will be switched but the instrument type will remain the same, so that the next visits of the instrument are processed in their original mode # only after all visits have been processed is the instrument type switched data_inst['type']=deepcopy(gen_dic['type'][inst]) if (inst=='NIGHT') and (data_inst['type']!='spec2D'):stop('ERROR : wrong data type for NIGHT (only S2D are available)') #Total number of orders for current instrument # - we define an artificial order that contains the CCF or 1D spectrum, so that the pipeline can process in the same way as with 2D spectra # - if orders are fully removed from an instrument its default structure is updated # if orders are trimmed after the spectral reduction the default structure is kept # - mock data created with a single order if gen_dic['mock_data']: data_inst['nord'] = 1 gen_dic[inst]['norders_instru'] = 1 else: norders_instru_ref = return_spec_nord(inst) data_inst['idx_ord_ref']=deepcopy(np.arange(norders_instru_ref)) #to keep track of the original orders indexes idx_ord_kept = list(np.arange(norders_instru_ref)) gen_dic[inst]['wav_ord_inst'] = return_cen_wav_ord(inst) if (data_inst['type'] in ['spec1D','CCF']): data_inst['nord'] = 1 gen_dic[inst]['norders_instru']=norders_instru_ref elif (data_inst['type']=='spec2D'): if (inst in gen_dic['del_orders']) and (len(gen_dic['del_orders'][inst])>0): print(' Removing '+str(len(gen_dic['del_orders'][inst]))+' orders from the analysis') idx_ord_kept = list(np.delete(np.arange(norders_instru_ref),gen_dic['del_orders'][inst])) gen_dic[inst]['wav_ord_inst'] = gen_dic[inst]['wav_ord_inst'][idx_ord_kept] data_inst['idx_ord_ref'] = data_inst['idx_ord_ref'][idx_ord_kept] gen_dic[inst]['norders_instru'] = len(idx_ord_kept) data_inst['nord']=deepcopy(len(idx_ord_kept)) data_inst['nord_spec']=deepcopy(data_inst['nord']) #to keep track of the number of orders after spectra are trimmed data_inst['nord_ref']=deepcopy(gen_dic[inst]['norders_instru']) #to keep track of the original number of orders (after deletion) #Orders contributing to calculation of spectral CCFs # - set automatically in case of S1D to use the generic pipeline structure, even though S1D have no orders # - indexes are made relative to order list after initial selection if data_inst['type']=='spec1D':gen_dic[inst]['orders4ccf']=[0] elif data_inst['type']=='spec2D': if (inst in gen_dic['orders4ccf']): gen_dic[inst]['orders4ccf'] = np.intersect1d(gen_dic['orders4ccf'][inst],idx_ord_kept,return_indices=True)[2] else:gen_dic[inst]['orders4ccf'] = np.arange(gen_dic[inst]['norders_instru'],dtype=int) #Weighing calibration condition # - S2D calibration is calculated for scaling even if not needed for weights # - condition below is deactivated if conversion into CCFs or 2D->1D data_inst['gcal_blaze_vis']=[] if (data_inst['type']=='spec2D') and gen_dic['cal_weight']:data_inst['cal_weight'] = True else:data_inst['cal_weight'] = False #Initialize flag that exposures in all visits of the instrument share a common spectral table data_inst['comm_sp_tab'] = True #Initialize list of instrument visits contained within a single night data_inst['single_night'] = [] #Calibration settings if (inst not in gen_dic['gcal_nooutedge']):gen_dic['gcal_nooutedge'][inst] = [0.,0.] #Processing each visit if gen_dic['mock_data']:vis_list=list(mock_dic['visit_def'][inst].keys()) else:vis_list=list(gen_dic['data_dir_list'][inst].keys()) data_inst['dates']={} data_inst['midpoints']={} ivis=-1 for vis in vis_list: #Remove/keep visits if (vis in gen_dic['unused_visits'][inst]): print(' > Visit "'+vis+'" is removed') else: print(' Initializing visit "'+vis+'"') data_dic[inst]['visit_list']+=[vis] ivis+=1 #Artificial data if gen_dic['mock_data']: #Generate time table if 'bin_low' in mock_dic['visit_def'][inst][vis]: bjd_exp_low = np.array(mock_dic['visit_def'][inst][vis]['bin_low']) bjd_exp_high = np.array(mock_dic['visit_def'][inst][vis]['bin_high']) n_in_visit = len(bjd_exp_low) elif 'exp_range' in mock_dic['visit_def'][inst][vis]: dbjd = (mock_dic['visit_def'][inst][vis]['exp_range'][1]-mock_dic['visit_def'][inst][vis]['exp_range'][0])/mock_dic['visit_def'][inst][vis]['nexp'] n_in_visit = int((mock_dic['visit_def'][inst][vis]['exp_range'][1]-mock_dic['visit_def'][inst][vis]['exp_range'][0])/dbjd) bjd_exp_low = mock_dic['visit_def'][inst][vis]['exp_range'][0] + dbjd*np.arange(n_in_visit) bjd_exp_high = bjd_exp_low+dbjd bjd_exp_all = 0.5*(bjd_exp_low+bjd_exp_high) #Observational data else: #Adding "/" in case user forgets vis_path_root = path_singslash(gen_dic['data_dir_list'][inst][vis]+'/') #List of all exposures for current instrument if inst in ['SOPHIE']: vis_path= vis_path_root+ { 'CCF':'*ccf*', 'spec2D':'*e2ds_', }[data_inst['type']] elif inst in ['CORALIE']: vis_path= vis_path_root+ { 'CCF':'*ccf*', 'spec2D':'*e2ds_', }[data_inst['type']] elif inst=='CARMENES_VIS': vis_path= vis_path_root+ { 'spec2D':'*sci-allr-vis_', }[data_inst['type']] elif inst=='CARMENES_VIS_CCF': #current reduction is not standard, better to consider as independent instrument vis_path= vis_path_root+ { 'CCF':'*CCF*', }[data_inst['type']] elif inst in ['ESPRESSO','ESPRESSO_MR','HARPN','HARPS']: vis_path= vis_path_root+ { 'CCF':'*CCF_', 'spec1D':'*S1D_', 'spec2D':'*S2D_', }[data_inst['type']] elif inst=='NIGHT': vis_path= vis_path_root+'*S2D_' elif inst in ['NIRPS_HA','NIRPS_HE']: vis_path= vis_path_root+ { 'CCF':'*CCF_TELL_CORR_', #CCF from telluric-corrected spectra 'spec1D':'*S1D_', 'spec2D':'*S2D_', }[data_inst['type']] elif inst in ['EXPRES']: vis_path= vis_path_root+ { 'spec2D':'*', }[data_inst['type']] else:stop('ERROR : root names undefined for '+inst) #Nominal data files if (inst not in ['EXPRES','CARMENES_VIS_CCF']):nom_ext = 'A.fits' else:nom_ext = '.fits' vis_path_exp = np.array(glob.glob(vis_path+nom_ext)) n_in_visit=len(vis_path_exp) if n_in_visit==0:stop('ERROR : No data found at '+vis_path+nom_ext+'. Check path.') #Retrieved file root names # - so that we can retrieve sky-corrected and blaze files in the same order exp_rootnames = [path.split(vis_path_root)[1].split(nom_ext)[0] for path in vis_path_exp] #Retrieving sky-corrected data if (inst in gen_dic['fibB_corr']) and (vis in gen_dic['fibB_corr'][inst]) and (len(gen_dic['fibB_corr'][inst][vis])>0): if inst in ['ESPRESSO','ESPRESSO_MR','HARPS','HARPN','NIRPS_HA','NIRPS_HE']:skcorr_ext = 'SKYSUB_A.fits' else:skcorr_ext = 'C.fits' vis_path_skysub_exp = np.array([ vis_path_root+ exp_rootname+skcorr_ext for exp_rootname in exp_rootnames ]) exp_filenames = np.array([ exp_rootname+skcorr_ext for exp_rootname in exp_rootnames ]) if len(vis_path_skysub_exp)==0:stop('ERROR : No sky-sub data found. Check path.') #Orders to be replaced if gen_dic['fibB_corr'][inst][vis]=='all':idx_ord_skysub = 'all' else:idx_ord_skysub = gen_dic['fibB_corr'][inst][vis] else: vis_path_skysub_exp = None exp_filenames = np.array([ exp_rootname+nom_ext for exp_rootname in exp_rootnames ]) #Retrieving blazed data if (data_inst['type']=='spec2D') and gen_dic['gcal_blaze']: if (inst in ['ESPRESSO','ESPRESSO_MR','HARPN','HARPS','NIRPS_HA','NIRPS_HE']): vis_path_blaze_exp = np.array([ vis_path_root+ exp_rootname+'BLAZE_A.fits' for exp_rootname in exp_rootnames ]) if len(vis_path_blaze_exp)==0:print('WARNING : No blaze data found. Switching to estimate mode.') else:data_inst['gcal_blaze_vis']+=[vis] elif inst=='CORALIE': vis_path_blaze_exp = [] for file_path in vis_path_exp: hdulist =fits.open(file_path) hdr =hdulist[0].header vis_path_blaze_exp+=[vis_path_root + hdr['ESO DRS BLAZE FILE']] vis_path_blaze_exp = np.array(vis_path_blaze_exp) if len(vis_path_blaze_exp)==0:print('WARNING : No blaze data found. Switching to estimate mode.') data_inst['gcal_blaze_vis']+=[vis] else:print(' Blaze names undefined for '+inst+'. Switching to estimate mode.') if (vis in data_inst['gcal_blaze_vis']):print(' Using blaze files for calibration profiles') #Print visit date if available # - taking the day at the start of the night, rather than the day at the start of the exposure series if (not gen_dic['mock_data']) and (inst in ['CORALIE','ESPRESSO','ESPRESSO_MR','HARPS','HARPN','CARMENES_VIS','CARMENES_VIS_CCF','EXPRES','NIGHT','NIRPS_HA','NIRPS_HE']): vis_year_exp_all=[] vis_month_exp_all=[] vis_day_exp_all=[] vis_hour_exp_all=[] bjd_exp_all = [] for file_path in vis_path_exp: hdulist =fits.open(file_path) hdr =hdulist[0].header if inst=='CORALIE': vis_year_exp_all+= [int(hdr['HIERARCH ESO CORA SHUTTER START DATE'][0:4])] vis_month_exp_all+= [int(hdr['HIERARCH ESO CORA SHUTTER START DATE'][4:6])] vis_day_exp_all+= [int(hdr['HIERARCH ESO CORA SHUTTER START DATE'][6:8])] vis_hour_exp_all+=[float(hdr['HIERARCH ESO CORA SHUTTER START HOUR'])] bjd_exp_all += [ hdr['HIERARCH ESO DRS BJD'] - 2400000. ] elif inst=='EXPRES': vis_year_exp_all+= [int(hdr['DATE-OBS'].split(' ')[0].split('-')[0])] vis_month_exp_all+= [int(hdr['DATE-OBS'].split(' ')[0].split('-')[1])] vis_day_exp_all+= [int(hdr['DATE-OBS'].split(' ')[0].split('-')[2]) ] vis_hour_exp_all+=[int(hdr['DATE-OBS'].split(' ')[1].split(':')[0])] elif inst=='CARMENES_VIS_CCF':bjd_exp_all += [ hdr['BJD'] ] elif inst in ['ESPRESSO','ESPRESSO_MR','HARPS','HARPN','CARMENES_VIS','NIGHT','NIRPS_HA','NIRPS_HE']: vis_year_exp_all+= [int(hdr['DATE-OBS'].split('T')[0].split('-')[0])] vis_month_exp_all+= [int(hdr['DATE-OBS'].split('T')[0].split('-')[1])] vis_day_exp_all+= [int(hdr['DATE-OBS'].split('T')[0].split('-')[2]) ] vis_hour_exp_all+=[int(hdr['DATE-OBS'].split('T')[1].split(':')[0])] if inst in ['ESPRESSO','ESPRESSO_MR','HARPS','HARPN','NIGHT','NIRPS_HA','NIRPS_HE']:bjd_exp_all +=[ hdr['HIERARCH '+facil_inst+' QC BJD'] - 2400000. ] elif inst=='CARMENES_VIS':bjd_exp_all += [ hdr['CARACAL BJD'] ] if inst=='CORALIE': vis_yr = hdr['HIERARCH ESO CORA SHUTTER START DATE'][0:4] vis_mt = hdr['HIERARCH ESO CORA SHUTTER START DATE'][4:6] elif inst=='EXPRES': vis_yr = hdr['DATE-OBS'].split(' ')[0].split('-')[0] vis_mt = hdr['DATE-OBS'].split(' ')[0].split('-')[1] elif (inst not in ['CARMENES_VIS_CCF']): vis_yr = hdr['DATE-OBS'].split('T')[0].split('-')[0] vis_mt = hdr['DATE-OBS'].split('T')[0].split('-')[1] #BJD midpoint bjd_exp_all = np.array(bjd_exp_all) bjd_vis = np.mean(bjd_exp_all) data_inst['midpoints'][vis] = '%.5f'%(bjd_vis+2400000.)+' BJD' #Take the day before as reference if exposure is past midnight (ie, not between 12 and 23) if len(vis_day_exp_all)>0: vis_year_exp_all = np.array(vis_year_exp_all) vis_month_exp_all = np.array(vis_month_exp_all) vis_day_exp_all = np.array(vis_day_exp_all) vis_hour_exp_all = np.array(vis_hour_exp_all) idx_day_shift = np_where1D(vis_hour_exp_all<=12) vis_day_exp_all[idx_day_shift]-=1 idx_month_shift =np_where1D(vis_day_exp_all == 0) vis_day_exp_all[idx_month_shift]=32 #for the purpose of setting the exposure day at the end of the previous month vis_month_exp_all[idx_month_shift]-=1 idx_year_shift =np_where1D(vis_month_exp_all == 0) vis_month_exp_all[idx_year_shift]=12 vis_year_exp_all[idx_year_shift]-=1 #Check wether visit is contained within a single night # - either: # + all exposures on same day before midnight or after midnight # + all exposures within two consecutive days between noon and midnight min_day = np.min(vis_day_exp_all) max_day = np.max(vis_day_exp_all) min_hr = np.min(vis_hour_exp_all) max_hr = np.max(vis_hour_exp_all) if ((max_day==min_day) and ((min_hr>=12) or (max_hr<=12))) or ((max_day==min_day+1) and ((min_hr>=12) and (max_hr<=12))): data_inst['single_night']+=[vis] vis_day = np.min(vis_day_exp_all) vis_day_txt = '0'+str(vis_day) if vis_day<10 else str(vis_day) data_inst['dates'][vis] = str(vis_year_exp_all[0])+'/'+str(vis_month_exp_all[0])+'/'+str(vis_day_txt) else: vis_day_txt_all = np.array(['0'+str(vis_day) if vis_day<10 else str(vis_day) for vis_day in vis_day_exp_all]) else: print('WARNING: data could not be retrieved') if (np.max(bjd_exp_all) - np.min(bjd_exp_all))<1.:data_inst['single_night']+=[vis] else: bjd_vis = np.mean(bjd_exp_all) - 2400000. data_inst['single_night']+=[vis] #Initializing dictionaries for visit theo_dic[inst][vis]={} data_dic['Atm'][inst][vis]={} data_inst[vis] = {'n_in_visit':n_in_visit,'studied_pl':[],'studied_ar':[],'comm_sp_tab':True} coord_dic[inst][vis] = {} for pl_loc in gen_dic['studied_pl_list']: if (inst in gen_dic['studied_pl'][pl_loc]) and (vis in gen_dic['studied_pl'][pl_loc][inst]):data_inst[vis]['studied_pl']+=[pl_loc] for ar in gen_dic['studied_ar']: if (inst in gen_dic['studied_ar'][ar]) and (vis in gen_dic['studied_ar'][ar][inst]):data_inst[vis]['studied_ar']+=[ar] data_prop[inst][vis] = {} data_dic_temp={} gen_dic[inst][vis] = {} DI_data_inst[vis] = {} #Associating telluric spectrum with each exposure if not gen_dic['mock_data']:data_inst[vis]['tell_DI_data_paths']={} #Exposure-specific calibration profile for weighing # - the path is created even if weights are not needed to store temporarily blaze-derived calibration profiles if data_inst['cal_weight'] or (vis in data_inst['gcal_blaze_vis']):data_inst[vis]['sing_gcal_DI_data_paths']={} #Initialize current data type for the visit data_inst[vis]['type']=deepcopy(data_inst['type']) #Initialize exposure tables # - spatial positions in x/y in front of the stellar disk # + x position along the node line # y position along the projection of the normal to the orbital plane # +calculated also for the out-of-transit positions, as in case of exposures binning some of the # ingress/egress exposures may contribute to the binned exposures average positions # - phase of in/out exposures, per instrument and per visit # - radial velocity of planet in star rest frame for key in ['bjd','t_dur','RV_star_solCDM','RV_star_stelCDM','cen_ph_st','st_ph_st','end_ph_st']:coord_dic[inst][vis][key] = np.zeros(n_in_visit,dtype=float)*np.nan for pl_loc in list(set(gen_dic['studied_pl_list']+gen_dic['kepl_pl'])): coord_dic[inst][vis][pl_loc]={} if pl_loc in data_inst[vis]['studied_pl']: for key in ['ecl','cen_ph','st_ph','end_ph','ph_dur','rv_pl','v_pl']:coord_dic[inst][vis][pl_loc][key] = np.zeros(n_in_visit,dtype=float)*np.nan for key in ['cen_pos','st_pos','end_pos']:coord_dic[inst][vis][pl_loc][key] = np.zeros([3,n_in_visit],dtype=float)*np.nan #Definition of mid-transit times for each planet associated with the visit if (pl_loc in gen_dic['Tcenter_visits']) and (inst in gen_dic['Tcenter_visits'][pl_loc]) and (vis in gen_dic['Tcenter_visits'][pl_loc][inst]): if vis not in data_inst['single_night']:stop('ERROR : gen_dic["Tcenter_visits"] not available for multi-epochs visits') coord_dic[inst][vis][pl_loc]['Tcenter'] = gen_dic['Tcenter_visits'][pl_loc][inst][vis] else: norb = round((bjd_vis+2400000.-system_param[pl_loc]['TCenter'])/system_param[pl_loc]["period"]) coord_dic[inst][vis][pl_loc]['Tcenter'] = system_param[pl_loc]['TCenter']+norb*system_param[pl_loc]["period"] #Observation properties for key in ['AM','IWV_AM','TEMP','PRESS','seeing','colcorrmin','colcorrmax','mean_SNR','alt','az','BERV']:data_prop[inst][vis][key] = np.zeros(n_in_visit,dtype=float)*np.nan data_prop[inst][vis]['exp_filenames']=np.empty([n_in_visit],dtype='U50') if inst=='ESPRESSO_MR': for key in ['AM_UT','seeing_UT']:data_prop[inst][vis][key] = np.zeros(n_in_visit,dtype=object)*np.nan data_prop[inst][vis]['PSF_prop'] = np.zeros([n_in_visit,4],dtype=float) data_prop[inst][vis]['SNRs']=np.zeros([n_in_visit,data_inst['nord_ref']],dtype=float)*np.nan if inst in ['ESPRESSO','HARPN','HARPS','NIRPS_HA','NIRPS_HE']: data_prop[inst][vis]['colcorr_ord']=np.zeros([n_in_visit,data_inst['nord_ref']],dtype=float)*np.nan for key in ['BLAZE_A','BLAZE_B','WAVE_MATRIX_THAR_FP_A','WAVE_MATRIX_THAR_FP_B']:data_prop[inst][vis][key]=np.empty([n_in_visit],dtype='U35') data_prop[inst][vis]['satur_check']=np.zeros([n_in_visit],dtype=int)*np.nan data_prop[inst][vis]['adc_prop']=np.zeros([n_in_visit,6],dtype=float)*np.nan data_prop[inst][vis]['piezo_prop']=np.zeros([n_in_visit,4],dtype=float)*np.nan if inst=='ESPRESSO': data_prop[inst][vis]['guid_coord']=np.zeros([n_in_visit,2],dtype=float)*np.nan #Set error flag per visit # - same for a given instrument, but this allows dealing with binned visits gen_dic[inst][vis]['flag_err']=deepcopy(gen_dic['flag_err_inst'][inst]) #Raw CCF properties if inst in ['HARPN','HARPS','CORALIE','SOPHIE','ESPRESSO','ESPRESSO_MR','NIRPS_HA','NIRPS_HE']: for key in ['RVdrift','rv_pip','FWHM_pip','ctrst_pip']:DI_data_inst[vis][key] = np.zeros(n_in_visit,dtype=float)*np.nan for key in ['erv_pip','eFWHM_pip','ectrst_pip']:DI_data_inst[vis][key] = np.zeros(n_in_visit,dtype=float) #Activity indexes data_inst[vis]['act_idx'] = [] if inst in ['HARPN','HARPS','CORALIE','SOPHIE','ESPRESSO','ESPRESSO_MR','NIRPS_HA','NIRPS_HE'] and gen_dic['DACE_sp']: DACE_sp = True #Corresponding DACE star name, instrument and visit dace_name = hdr['HIERARCH '+facil_inst+' OBS TARG NAME'] DACE_sp = Spectroscopy.get_timeseries(dace_name, sorted_by_instrument=False) #Condition to retrieve current visit dace_red_inst = hdr['INSTRUME'] if vis in data_inst['single_night']: if dace_red_inst=='ESPRESSO': if bjd_vis<=58649.5:dace_red_inst+='18' else:dace_red_inst+='19' elif dace_red_inst=='HARPS': if bjd_vis<=57167.5:dace_red_inst+='03' else:dace_red_inst+='15' idx_dace_vis = np_where1D((np.asarray(DACE_sp['ins_name']) == dace_red_inst) & (np.asarray(DACE_sp['date_night']) == vis_yr+'-'+vis_mt+'-'+vis_day_txt)) else: if dace_red_inst=='ESPRESSO':inst_off_all = ['18' if bjd_loc<58649.5 else '19' for bjd_loc in bjd_exp_all] elif dace_red_inst=='HARPS':inst_off_all = ['03' if bjd_loc<57167.5 else '15' for bjd_loc in bjd_exp_all] dace_red_inst_all = np.array([dace_red_inst+inst_off for inst_off in inst_off_all]) idx_dace_vis = np.zeros(0,dtype=int) vis_day_txt_un,idx_vis_day_txt_un = np.unique(vis_day_txt_all,return_index=True) for vis_day_txt,dace_red_inst in zip(vis_day_txt_un,dace_red_inst_all[idx_vis_day_txt_un]): idx_dace_vis = np.append(idx_dace_vis, np_where1D((np.asarray(DACE_sp['ins_name']) == dace_red_inst) & (np.asarray(DACE_sp['date_night']) == vis_yr+'-'+vis_mt+'-'+vis_day_txt))) if len(idx_dace_vis)==0: print(' Visit not found in DACE') DACE_sp = False else: #Retrieve activity indexes and bjd DACE_idx = {} DACE_idx['bjd'] = np.array(DACE_sp['rjd'],dtype=float)[idx_dace_vis] data_inst[vis]['act_idx'] = ['ha','na','ca','s','rhk'] for key in data_inst[vis]['act_idx']: data_prop[inst][vis][key] = np.zeros([n_in_visit,2],dtype=float)*np.nan if key=='rhk':dace_key = 'rhk' else:dace_key = key+'index' DACE_idx[key] = np.array(DACE_sp[dace_key],dtype=float)[idx_dace_vis] DACE_idx[key+'_err'] = np.array(DACE_sp[dace_key+'_err'],dtype=float)[idx_dace_vis] else: DACE_sp = False if inst=='EXPRES': data_inst[vis]['act_idx'] = ['ha','s'] for key in data_inst[vis]['act_idx']:data_prop[inst][vis][key] = np.zeros([n_in_visit,2],dtype=float)*np.nan #------------------------------------------------------------------------------------------------------------ #Simulated active regions if (len(data_inst[vis]['studied_ar'])>0): if gen_dic['mock_data']: if (inst in mock_dic['ar_prop']) and (vis in mock_dic['ar_prop'][inst]):ar_prop_nom = retrieve_ar_prop_from_param(mock_dic['ar_prop'][inst][vis], inst, vis) else:stop('ERROR: active regions are required in visit '+vis+' but their mock properties are not defined') else: if (inst in theo_dic['ar_prop']) and (vis in theo_dic['ar_prop'][inst]):ar_prop_nom = retrieve_ar_prop_from_param(theo_dic['ar_prop'][inst][vis], inst, vis) else:stop('ERROR: active regions are required in visit '+vis+' but their theoretical properties are not defined') ar_prop_nom['cos_istar']=system_param['star']['cos_istar'] for ar in data_inst[vis]['studied_ar']: coord_dic[inst][vis][ar]={} for key in gen_dic['ar_coord_par']:coord_dic[inst][vis][ar][key] = np.zeros([3,n_in_visit],dtype=float)*np.nan coord_dic[inst][vis][ar]['is_visible'] = np.zeros([3,n_in_visit],dtype=bool) for key in ['Tc_ar', 'ang_rad','lat_rad', 'fctrst']:coord_dic[inst][vis][ar][key] = ar_prop_nom[ar][key] #Pre-processing all exposures in visit for isub_exp,iexp in enumerate(range(n_in_visit)): #Artificial data if gen_dic['mock_data']: coord_dic[inst][vis]['bjd'][iexp] = bjd_exp_all[iexp] - 2400000. coord_dic[inst][vis]['t_dur'][iexp] = (bjd_exp_high[iexp]-bjd_exp_low[iexp])*24.*3600. #Observational data else: #File name data_prop[inst][vis]['exp_filenames'][iexp] = exp_filenames[iexp] #Data of the fits file hdulist =fits.open(vis_path_exp[iexp]) #Header hdr =hdulist[0].header #Retrieve exposure bjd in instrument/visit table if inst in ['HARPS','HARPN','ESPRESSO','ESPRESSO_MR','NIGHT','NIRPS_HA','NIRPS_HE']: bjd_exp = hdr['HIERARCH '+facil_inst+' QC BJD'] - 2400000. elif inst in ['CORALIE']:bjd_exp = hdr['HIERARCH ESO DRS BJD'] - 2400000. elif inst in ['SOPHIE']:bjd_exp = hdr['HIERARCH OHP DRS BJD'] - 2400000. elif inst in ['CARMENES_VIS']:bjd_exp = hdr['CARACAL BJD'] elif inst=='EXPRES':bjd_exp = hdulist[1].header['BARYMJD'] +0.5 else:stop('ERROR : bjd keyword for '+inst+' is undefined') coord_dic[inst][vis]['bjd'][iexp] = bjd_exp #Exposure time (s) if inst=='SOPHIE': #start time st_hmns=(hdr['HIERARCH OHP OBS DATE START'].split('T')[1]).split(':') h_start=float(st_hmns[0]) mn_start=float(st_hmns[1]) s_start=float(st_hmns[2]) #end time end_hmns=(hdr['HIERARCH OHP OBS DATE END'].split('T')[1]).split(':') h_end=float(end_hmns[0]) if h_end<h_start:h_end+=24. mn_end=float(end_hmns[1]) s_end=float(end_hmns[2]) #exposure time (s) coord_dic[inst][vis]['t_dur'][iexp]=( (h_end-h_start)*3600. + (mn_end-mn_start)*60. + (s_end-s_start) ) elif inst=='CORALIE':coord_dic[inst][vis]['t_dur'][iexp] = hdr['HIERARCH ESO OBS TEXP'] elif inst=='EXPRES':coord_dic[inst][vis]['t_dur'][iexp] = hdr['AEXPTIME'] else:coord_dic[inst][vis]['t_dur'][iexp] = hdr['EXPTIME'] #RV drift # - drift is measured from the calibration lamps or FP # - drift is corrected automatically in new reductions for ESPRESSO and kin spectrographs, direcly included in the wavelength solution if inst in ['CORALIE']: DI_data_inst[vis]['RVdrift'][iexp]=hdr['HIERARCH ESO DRS DRIFT SPE RV'] #in m/s elif inst=='SOPHIE': DI_data_inst[vis]['RVdrift'][iexp]=hdr['HIERARCH OHP DRS DRIFT RV'] #in m/s elif inst in ['HARPS','HARPN','ESPRESSO','ESPRESSO_MR','NIRPS_HA','NIRPS_HE']: if 'HIERARCH '+facil_inst+' QC DRIFT DET0 MEAN' in hdr: DI_data_inst[vis]['RVdrift'][iexp]=hdr['HIERARCH '+facil_inst+' QC DRIFT DET0 MEAN']*return_pix_size(inst)*1e3 #in pix -> m/s #Stellar phase if gen_dic['sequence'] not in ['st_master_tseries']: coord_dic[inst][vis]['st_ph_st'][iexp],coord_dic[inst][vis]['cen_ph_st'][iexp],coord_dic[inst][vis]['end_ph_st'][iexp] = get_timeorbit(system_param['star']['Tcenter'],coord_dic[inst][vis]['bjd'][iexp], {'period':system_param['star']['Peq']}, coord_dic[inst][vis]['t_dur'][iexp])[0:3] #Orbital coordinates for each studied planet for pl_loc in data_inst[vis]['studied_pl']: coord_dic[inst][vis][pl_loc]['cen_pos'][:,iexp],coord_dic[inst][vis][pl_loc]['st_pos'][:,iexp],coord_dic[inst][vis][pl_loc]['end_pos'][:,iexp],coord_dic[inst][vis][pl_loc]['ecl'][iexp],coord_dic[inst][vis][pl_loc]['rv_pl'][iexp],coord_dic[inst][vis][pl_loc]['v_pl'][iexp],\ coord_dic[inst][vis][pl_loc]['st_ph'][iexp],coord_dic[inst][vis][pl_loc]['cen_ph'][iexp],coord_dic[inst][vis][pl_loc]['end_ph'][iexp],coord_dic[inst][vis][pl_loc]['ph_dur'][iexp]=coord_expos(pl_loc,coord_dic,inst,vis,system_param['star'], system_param[pl_loc],coord_dic[inst][vis]['bjd'][iexp],coord_dic[inst][vis]['t_dur'][iexp],data_dic,data_dic['DI']['system_prop']['achrom'][pl_loc][0]) #Surface coordinates for each studied active region for ar in data_inst[vis]['studied_ar']: ar_prop_exp = coord_expos_ar(ar,coord_dic[inst][vis]['bjd'][iexp],ar_prop_nom,system_param['star'],coord_dic[inst][vis]['t_dur'][iexp],gen_dic['ar_coord_par']) for key in ar_prop_exp:coord_dic[inst][vis][ar][key][:, iexp] = [ar_prop_exp[key][0],ar_prop_exp[key][1],ar_prop_exp[key][2]] #-------------------------------------------------------------------------------------------------- #Processing all exposures in visit for isub_exp,iexp in enumerate(range(n_in_visit)): #Artificial data if gen_dic['mock_data']: #Initialize data at first exposure if isub_exp==0: print(' Building exposures ... ') data_inst[vis]['mock'] = True fixed_args = {} if inst not in mock_dic['intr_prof']: print(' Automatic definition of mock line profile') mock_dic['intr_prof'][inst] = {'mode':'ana','coord_line':'mu','model': 'gauss', 'line_trans':None,'mod_prop':{'ctrst_ord0__IS__VS_' : 0.5,'FWHM_ord0__IS__VS_' : 5. },'pol_mode' : 'modul'} fixed_args.update(mock_dic['intr_prof'][inst]) #Activation of spectral conversion and resampling cond_conv_st_prof_tab(theo_dic['rv_osamp_line_mod'],fixed_args,data_inst[vis]['type']) #Mock data spectral table # - fixed for all exposures, defined in the star rest frame # - the model table is defined in the star rest frame, so that models calculated / defined on this table remain centered independently of the systemic and keplerian RV # models for each exposure are then resampled on the original spectral table if required, and the table is only then shifted / scaled by the systemic and keplerian RV # - theoretical models are defined in wavelength space # analytical models are defined in RV space, and their table converted into spectral space afterwards if relevant data_inst[vis]['nspec'] = int(np.ceil((mock_dic['DI_table']['x_end']-mock_dic['DI_table']['x_start'])/mock_dic['DI_table']['dx'])) delta_cen_bins = mock_dic['DI_table']['x_end'] - mock_dic['DI_table']['x_start'] dcen_bins = delta_cen_bins/data_inst[vis]['nspec'] fixed_args['cen_bins']=mock_dic['DI_table']['x_start']+dcen_bins*(0.5+np.arange(data_inst[vis]['nspec'])) fixed_args['edge_bins']=def_edge_tab(fixed_args['cen_bins'][None,:][None,:])[0,0] fixed_args['ncen_bins']=data_inst[vis]['nspec'] fixed_args['dcen_bins']=fixed_args['edge_bins'][1::] - fixed_args['edge_bins'][0:-1] fixed_args['dim_exp']=[data_inst['nord'],fixed_args['ncen_bins']] #Resampled spectral table for model line profile if fixed_args['resamp']:resamp_st_prof_tab(None,None,None,fixed_args,gen_dic,1,theo_dic['rv_osamp_line_mod']) #Effective instrumental convolution fixed_args['FWHM_inst'] = get_FWHM_inst(inst,fixed_args,fixed_args['cen_bins']) #Initialize intrinsic profile properties # - removing Peq, as it is always defined in the nominal stellar properties, but used as condition to switch from veq to Peq within var_stellar_prop() when Peq is defined as model parameter instead of veq params_mock = deepcopy(system_param['star']) params_mock.pop('Peq') if inst not in mock_dic['flux_cont']:mock_dic['flux_cont'][inst]={} if vis not in mock_dic['flux_cont'][inst]:mock_dic['flux_cont'][inst][vis] = 1. params_mock.update({'rv':0.,'cont':mock_dic['flux_cont'][inst][vis]}) params_mock = par_formatting(params_mock,fixed_args['mod_prop'],None,None,fixed_args,inst,vis) params_mock = par_formatting_inst_vis(params_mock,fixed_args,inst,vis,mock_dic['intr_prof'][inst]['mode']) #Generic properties required for model calculation if inst not in mock_dic['sysvel']:mock_dic['sysvel'][inst]={} if vis not in mock_dic['sysvel'][inst]:mock_dic['sysvel'][inst][vis] = 0. fixed_args.update({ 'mac_mode':theo_dic['mac_mode'], 'type':data_inst[vis]['type'], 'nord':data_inst['nord'], 'nthreads':mock_dic['nthreads'], 'unthreaded_op':mock_dic['unthreaded_op'], 'resamp_mode' : gen_dic['resamp_mode'], 'conv2intr':False, 'inst':inst, 'vis':vis, 'fit':False, 'system_param':system_param, 'ar_coord_par':gen_dic['ar_coord_par'], 'rout_mode':'Intr_prop', }) #Active region properties if (inst in mock_dic['ar_prop']) and (data_dic['DI']['ar_prop'] != {}): params_mock['use_ar']=True ar_prop = data_dic['DI']['ar_prop'] else: params_mock['use_ar']=False ar_prop = {} #Initializing stellar properties fixed_args = var_stellar_prop(fixed_args,theo_dic,data_dic['DI']['system_prop'],ar_prop,system_param['star'],params_mock) #Observational data else: #Observational data for current exposure hdulist =fits.open(vis_path_exp[iexp]) if vis_path_skysub_exp is not None:hdulist_skysub =fits.open(vis_path_skysub_exp[iexp]) hdr =hdulist[0].header #VLT UT used by ESPRESSO if (inst=='ESPRESSO') and ((isub_exp==0) or (vis not in data_inst['single_night'])): tel_inst = { 'ESO-VLT-U1':'1', 'ESO-VLT-U2':'2', 'ESO-VLT-U3':'3', 'ESO-VLT-U4':'4' }[hdr['TELESCOP']] #Initialize data at first exposure if isub_exp==0: data_inst[vis]['mock'] = False #Instrumental mode if inst in ['ESPRESSO','HARPS','NIRPS_HA','NIRPS_HE']: data_prop[inst][vis]['ins_mod'] = "".join(hdr['HIERARCH ESO INS MODE'].split()) else:stop('ERROR : instrumental mode undefined for '+inst) #Binning factor along X if inst in ['ESPRESSO']: data_prop[inst][vis]['det_binx'] = hdr['HIERARCH ESO DET BINX'] elif inst in ['HARPS','NIRPS_HA','NIRPS_HE']: data_prop[inst][vis]['det_binx'] = hdr['HIERARCH ESO DET WIN1 BINX'] #Radial velocity tables for input CCFs # - assumed to be common for all CCFs of the visit dataset, and thus calculated only once if (data_inst['type']=='CCF'): if inst in ['HARPN','HARPS','ESPRESSO','ESPRESSO_MR','NIRPS_HA','NIRPS_HE']: start_rv = hdr['HIERARCH '+facil_inst+' RV START'] # first velocity delta_rv = hdr['HIERARCH '+facil_inst+' RV STEP'] # delta vel step n_vel_vis = (hdulist[1].header)['NAXIS1'] else: start_rv = hdr['CRVAL1'] # first wavelength/velocity delta_rv = hdr['CDELT1'] # delta lambda/vel step n_vel_vis = hdr['NAXIS1'] #Velocity table and resolution velccf = start_rv+delta_rv*np.arange(n_vel_vis) #Screening # - check for oversampling, screening the profiles so that they are defined over uncorrelated points with a resolution similar to the instrumental pixel width. # - CCFs are screened (ie, keeping one point every specified length) with a length defined as input, or the instrumental bin length # the goal is to keep only uncorrelated points, on which we can measure the dispersion as a white noise error # - we remove scr_lgth-1 points between two kept points, ie we keep one point in scr_lgth+1 # if the original table is x0,x1,x2,x3,x4.. and scr_lgth=3 we keep x0,x3.. # - we define in 'scr_lgth' the number of CCF bins to be removed depending on the chosen length, and the CCF resolution # - we assume the same value can be used for a given visit: # + if the bin size is chosen as typical length, then it will not depend on the visits but the CCFs may have been computed with a # different resolution for each visit. In that case the number of screened bins is calculated once for the visit as: # scr_lgth = bin size/CCF_step # with CCF typically oversampled, CCF_step is lower than the bin size and we keep about one CCF point every bin size # + if the correlation length is chosen, then it may depend on the visit conditions but should be constant for a given visit. In that case # the number of screened bins has already been set in 'scr_lgth' for the visit, using the empirical estimation on the differential CCFs if vis not in gen_dic['scr_lgth'][inst]: gen_dic[inst][vis]['scr_lgth']=int(round(return_pix_size(inst)/delta_rv)) if gen_dic[inst][vis]['scr_lgth']<=1:gen_dic[inst][vis]['scr_lgth']=1 if gen_dic[inst][vis]['scr_lgth']==1:print(' No CCF screening required') else:print(' Screening by '+str(gen_dic[inst][vis]['scr_lgth'])+' pixels') if gen_dic[inst][vis]['scr_lgth']>1: idx_scr_bins=np.arange(n_vel_vis,dtype=int)[gen_dic['ist_scr']::gen_dic[inst][vis]['scr_lgth']] velccf = velccf[idx_scr_bins] #Table dimensions data_inst[vis]['nspec'] = len(velccf) #---------------------------------------- #Number of wavelength bins for spectra elif ('spec' in data_inst['type']): #1D spectra if (data_inst['type']=='spec1D'): if inst in ['ESPRESSO','ESPRESSO_MR','HARPN']:data_inst[vis]['nspec'] = (hdulist[1].header)['NAXIS2'] else:stop('TBD') #2D spectra elif (data_inst['type']=='spec2D'): if inst in ['ESPRESSO','ESPRESSO_MR','HARPN','HARPS','CARMENES_VIS','NIGHT','NIRPS_HA','NIRPS_HE']:data_inst[vis]['nspec'] = (hdulist[1].header)['NAXIS1'] elif inst=='CORALIE':data_inst[vis]['nspec'] = hdr['NAXIS1'] elif inst=='EXPRES':data_inst[vis]['nspec'] = 7920 else:stop('ERROR : number of pixels for '+inst+' undefined (set header key)') #Initialize data at first exposure if isub_exp==0: #Table dimensions data_inst[vis]['dim_all'] = [n_in_visit,data_inst['nord'],data_inst[vis]['nspec']] data_inst[vis]['dim_exp'] = [data_inst['nord'],data_inst[vis]['nspec']] data_inst[vis]['dim_sp'] = [n_in_visit,data_inst['nord']] data_inst[vis]['dim_ord'] = [n_in_visit,data_inst[vis]['nspec']] #Initialize spectral table, flux table, banded covariance matrix data_dic_temp['cen_bins'] = np.zeros(data_inst[vis]['dim_all'], dtype=float) data_dic_temp['flux'] = np.zeros(data_inst[vis]['dim_all'], dtype=float)*np.nan data_dic_temp['cov'] = np.zeros(data_inst[vis]['dim_sp'], dtype=object) #Initialize blaze-derived and detector noise profiles # - detector noise is only required for weighing, but is retrieved by default to avoid recomputing the calibration module if it was ran without weighing activated if (vis in data_inst['gcal_blaze_vis']): data_dic_temp['gcal'] = np.zeros(data_inst[vis]['dim_all'], dtype=float)*np.nan data_dic_temp['sdet2'] = np.zeros(data_inst[vis]['dim_all'], dtype=float)*np.nan #Telluric spectrum if ('spec' in data_inst[vis]['type']) and gen_dic['tell_weight'] and (gen_dic['calc_tell_mode']=='input'): data_dic_temp['tell'] = np.ones(data_inst[vis]['dim_all'], dtype=float) ### End of first exposure #------------------------------------------------------------------------------------ #Retrieve data for current exposure #------------------------------------------------------------------------------------ #----------------------------------- #Artificial data #----------------------------------- if gen_dic['mock_data']: #print(' building exposure '+str(iexp)+'/'+str(n_in_visit - 1)) param_exp = deepcopy(params_mock) #Table for model calculation args_exp = def_st_prof_tab(None,None,None,deepcopy(fixed_args)) #Initializing stellar profiles args_exp = init_custom_DI_prof(args_exp,gen_dic,param_exp) #Initializing broadband scaling of intrinsic profiles into local profiles # - defined in forward mode at initialization, or defined in fit mode only if the stellar grid is not updated through the fit # - there are no default pipeline tables for this scaling because they depend on the local spectral tables of the line profiles args_exp['Fsurf_grid_spec'] = theo_intr2loc(args_exp['grid_dic'],args_exp['system_prop'],args_exp,args_exp['ncen_bins'],theo_dic['nsub_star']) #Add jitter to the intrinsic profile properties (simulating stellar activity) if (fixed_args['mode']=='ana') and (inst in mock_dic['drift_intr']) and (vis in mock_dic['drift_intr'][vis]) and (len(mock_dic['drift_intr'][inst][vis]>0)): for par_drift in mock_dic['drift_intr'][inst][vis] : if par_drift in param_exp: if (par_drift=='rv'):param_exp[par_drift] += mock_dic['drift_intr'][inst][vis][par_drift][iexp] else:param_exp[par_drift] *= mock_dic['drift_intr'][inst][vis][par_drift][iexp] #Disk-integrated stellar line base_DI_prof = custom_DI_prof(param_exp,None,args=args_exp)[0] #Deviation from nominal stellar profile surf_prop_dic, surf_prop_dic_ar,_ = sub_calc_plocc_ar_prop([data_dic['DI']['system_prop']['chrom_mode']],args_exp,['line_prof'],data_dic[inst][vis]['studied_pl'],data_dic[inst][vis]['studied_ar'],deepcopy(system_param),theo_dic,args_exp['system_prop'],param_exp,coord_dic[inst][vis],[iexp], system_ar_prop_in=args_exp['system_ar_prop']) #Adding to base profile the deviations from planet and active region contributions DI_prof_exp = base_DI_prof - surf_prop_dic[data_dic['DI']['system_prop']['chrom_mode']]['line_prof'][:,0] - surf_prop_dic_ar[data_dic['DI']['system_prop']['chrom_mode']]['line_prof'][:,0] #Instrumental response # - in RV space for analytical model, in wavelength space for theoretical profiles # - resolution can be modified to model systematic variations from the instrument or atmosphere # - disabled if measured profiles as used as proxy for the intrinsic profiles if (fixed_args['mode']!='Intrbin') and (inst in mock_dic['drift_post']) and (vis in mock_dic['drift_post'][vis]) and ('resol' in mock_dic['drift_post'][inst][vis]): fixed_args['FWHM_inst'] = fixed_args['ref_conv']/mock_dic['drift_post'][inst][vis]['resol'][iexp] #Convolution, conversion and resampling DI_prof_exp = conv_st_prof_tab(None,None,None,fixed_args,args_exp,DI_prof_exp,fixed_args['FWHM_inst']) #Set negative flux values to null DI_prof_exp[DI_prof_exp<0.] = 0. #Define number of photoelectrons extracted during the exposure # - the model is a density of photoelectrons per unit of time, with continuum set to the input mean flux density if (inst in mock_dic['gcal']):mock_gcal = mock_dic['gcal'][inst] else:mock_gcal = 1. DI_prof_exp_Ftrue = mock_gcal*DI_prof_exp*coord_dic[inst][vis]['t_dur'][iexp] #Keplerian motion and systemic shift of the disk-integrated profile # - we shift profiles from the star rest frame (source) to the solar barycentric rest frame (receiver) # see gen_specdopshift() : # w_receiver = w_source * (1+ rv[s/r]/c)) # w_solbar = w_star * (1+ rv[star/starbar]/c))* (1+ rv[starbar/solbar]/c)) # we include variations in the systemic rv if requested RV_star_stelCDM_mock = calc_rv_star(coord_dic,inst,vis,system_param,gen_dic, coord_dic[inst][vis]['bjd'][iexp],coord_dic[inst][vis]['t_dur'][iexp],mock_dic['sysvel'][inst][vis])[0] sysvel_mock = mock_dic['sysvel'][inst][vis] if (inst in mock_dic['drift_post']) and (vis in mock_dic['drift_post'][vis]) and ('rv' in mock_dic['drift_post'][inst][vis]):sysvel_mock+=mock_dic['drift_post'][inst][vis]['rv'][iexp] if ('spec' in data_inst[vis]['type']): data_dic_temp['cen_bins'][iexp,0]=fixed_args['cen_bins']*gen_specdopshift(RV_star_stelCDM_mock)*gen_specdopshift(sysvel_mock) else:data_dic_temp['cen_bins'][iexp,0] = fixed_args['cen_bins'] + (RV_star_stelCDM_mock+sysvel_mock) #Defining flux and error table # - see weights_bin_prof(), the measured (total, not density) flux can be defined as: # F_meas(t,w) = gcal(band) Nmeas(t,w) # where Nmeas(t,w), drawn from a Poisson distribution with number of events Ntrue(t,w), is the number of photo-electrons measured durint texp # the estimate of the error is # EF_meas(t,w)^2 = gcal(band) F_meas(t,w) + Ern^2 ~ gcal(band) F_meas(t,w) # which is a biased estimate of the true error but corresponds to what is returned by DRS # we neglect read noise, also because only S2D tables are available in ANTARESS # - the S/N of the mock profiles is F_meas(t,w)/EF_meas(t,w) = sqrt(F_meas(t,w)/gcal(band)) = sqrt(Nmeas(t,w)) proportional to sqrt(texp*cont) # errors in the continuum of the normalized profiles are proportional to 1/sqrt(texp*cont) # beware that this error does not necessarily match the flux dispersion if (inst in mock_dic['set_err']) and mock_dic['set_err'][inst]: DI_prof_exp_Fmeas = np.array(list(map(np.random.poisson, DI_prof_exp_Ftrue, data_inst[vis]['nspec']*[1]))).flatten() else: DI_prof_exp_Fmeas = DI_prof_exp_Ftrue DI_err_exp_Emeas2 = mock_gcal*DI_prof_exp_Fmeas data_dic_temp['flux'][iexp,0] = DI_prof_exp_Fmeas data_dic_temp['cov'][iexp,0] = DI_err_exp_Emeas2[None,:] if mock_dic['verbose_flux_cont']:print(' Exposure %s: Average SNR = %.0f, Average photon count = %.0f' % (str(isub_exp), np.mean(DI_prof_exp_Fmeas/np.sqrt(DI_err_exp_Emeas2)), np.mean(DI_prof_exp_Fmeas))) #----------------------------------- #Observational data #----------------------------------- else: cond_bad_exp = np.zeros(data_inst[vis]['dim_exp'],dtype=bool) #Retrieve BERV # - in km/s if inst=='EXPRES': data_prop[inst][vis]['BERV'][iexp] = hdulist[2].header['HIERARCH wtd_mdpt_bc']*c_light else: if inst=='CARMENES_VIS': reduc_txt = 'CARACAL' elif inst in ['HARPN','HARPS','ESPRESSO','ESPRESSO_MR','NIGHT','NIRPS_HA','NIRPS_HE']:reduc_txt = facil_inst+' QC' elif inst=='CORALIE':reduc_txt = facil_inst+' DRS' else:stop('ERROR : define BERV retrieval') data_prop[inst][vis]['BERV'][iexp] = hdr['HIERARCH '+reduc_txt+' BERV'] #Retrieve CCFs # - HARPS 'hdu' has 73 x n_velocity size # + 72 CCFs for each of the spectrograph order # + last element is the total CCF over all orders (n_velocity size) # - HARPN 'hdu' has 70 x n_velocity size # + 69 CCFs for each of the spectrograph order # + last element is the total CCF over all orders (n_velocity size) # - SOPHIE 'hdu' has 40 x n_velocity size # + 39 CCFs for each of the spectrograph order # + last element is the total CCF over all orders (n_velocity size) if data_inst['type']=='CCF': #Bin centers # - we associate the velocity table to all exposures to keep the same structure as spectra data_dic_temp['cen_bins'][iexp,0] = velccf #CCF tables per order if (vis_path_skysub_exp is not None) and (gen_dic['fibB_corr'][inst][vis]=='all'):hdulist_dat= hdulist_skysub else:hdulist_dat = hdulist if inst in ['HARPS','ESPRESSO','ESPRESSO_MR','HARPN','NIRPS_HA','NIRPS_HE']: all_ccf = hdulist_dat[1].data all_eccf = hdulist_dat[2].data if (vis_path_skysub_exp is not None) and (gen_dic['fibB_corr'][inst][vis]!='all'): all_ccf[idx_ord_skysub] = (hdulist_skysub[1].data)[idx_ord_skysub] all_eccf[idx_ord_skysub] = (hdulist_skysub[2].data)[idx_ord_skysub] else: all_ccf =hdulist_dat[0].data if (vis_path_skysub_exp is not None) and (gen_dic['fibB_corr'][inst][vis]!='all'): all_ccf[idx_ord_skysub] = (hdulist_skysub[0].data)[idx_ord_skysub] #Construction of the total CCF # - we re-calculate the CCF rather than taking the total CCF stored in the last column of the matrix to be able to account for sky-correction in specific orders # - we use the input "gen_dic['orders4ccf'][inst]" rather than "gen_dic[inst]['orders4ccf']", which is used for CCF computed from spectral data if inst not in gen_dic['orders4ccf']:flux_raw = np.nansum(all_ccf[ range(return_spec_nord(inst))],axis=0) else:flux_raw = np.nansum(all_ccf[gen_dic['orders4ccf'][inst]],axis=0) #Screening CCF if gen_dic[inst][vis]['scr_lgth']>1:flux_raw = flux_raw[idx_scr_bins] data_dic_temp['flux'][iexp,0] = flux_raw #Error table if gen_dic[inst][vis]['flag_err']: if inst not in gen_dic['orders4ccf']:err_raw = np.sqrt(np.nansum(all_eccf[range(return_spec_nord(inst))]**2.,axis=0)) else:err_raw = np.sqrt(np.nansum(all_eccf[gen_dic['orders4ccf'][inst]]**2.,axis=0)) if gen_dic[inst][vis]['scr_lgth']>1:err_raw = err_raw[idx_scr_bins] err_raw = np.tile(err_raw,[data_inst['nord'],1]) #------------------------------------------------------------------------------------------------------------------------------------------- #Retrieve spectral data elif ('spec' in data_inst['type']): #1D spectra if data_inst['type']=='spec1D': if inst in ['HARPS','ESPRESSO','ESPRESSO_MR','HARPN','NIRPS_HA','NIRPS_HE']: #Replacing single order with sky-corrected data if (vis_path_skysub_exp is not None):data_loc = hdulist_skysub[1].data else:data_loc = hdulist[1].data #Bin centers if gen_dic['sp_frame']=='air':data_dic_temp['cen_bins'][iexp] = data_loc['wavelength_air'] elif gen_dic['sp_frame']=='vacuum':data_dic_temp['cen_bins'][iexp] = data_loc['wavelength'] #Spectra # - flux and errors are per unit of wavelength but not per unit of time, ie that they correspond to the number of photoelectrons measured during the exposure data_dic_temp['flux'][iexp] = data_loc['flux'] err_raw = data_loc['error'] #Pixels quality qualdata_exp = data_loc['quality'] else: stop('ERROR : spectra upload undefined for '+inst) #------------------------------ #2D spectra elif data_inst['type']=='spec2D': #Sky-corrected data if (vis_path_skysub_exp is not None): if (gen_dic['fibB_corr'][inst][vis]=='all'):hdulist_dat= hdulist_skysub else: cond_skysub = np.repeat(False,return_spec_nord(inst)) #Conditions on original order indexes cond_skysub[idx_ord_skysub] = True idxsub_ord_skysub = np_where1D(cond_skysub[idx_ord_kept]) #Indexes in reduced tables else:hdulist_dat = hdulist if inst in ['HARPS','ESPRESSO','ESPRESSO_MR','HARPN','NIGHT','NIRPS_HA','NIRPS_HE']: #Bin centers # - dimension norder x nbins if gen_dic['sp_frame']=='air':ikey_wav = 5 elif gen_dic['sp_frame']=='vacuum':ikey_wav = 4 data_dic_temp['cen_bins'][iexp] = (hdulist_dat[ikey_wav].data)[idx_ord_kept] dll = (hdulist_dat[ikey_wav+2].data)[idx_ord_kept] if (vis_path_skysub_exp is not None) and (gen_dic['fibB_corr'][inst][vis]!='all'): data_dic_temp['cen_bins'][iexp][idxsub_ord_skysub] = (hdulist_skysub[ikey_wav].data)[idx_ord_kept][idxsub_ord_skysub] dll[idxsub_ord_skysub] = (hdulist_skysub[ikey_wav+2].data)[idx_ord_kept][idxsub_ord_skysub] #Spectra # - the overall flux level is changed in the sky-corrected data, messing up with the flux balance corrections if only some orders were corrected # we thus roughly rescale the flux in the corrected orders to their original level (assuming constant resolution over the order to calculate the mean flux) # - the flux is per unit of pixel, and must be divided by the pixel spectral size stored in dll to be passed in units of wavelength # - the flux is not per unit of time, but corresponds to the number of photoelectrons measured during the exposure corrected for the blaze function data_dic_temp['flux'][iexp] = (hdulist_dat[1].data)[idx_ord_kept]/dll err_raw = (hdulist_dat[2].data)[idx_ord_kept]/dll #Replacing with sky-corrected data if (vis_path_skysub_exp is not None) and (gen_dic['fibB_corr'][inst][vis]!='all'): dll_ord_skysub = dll[idxsub_ord_skysub] mean_Fraw_ord_skysub = np.nansum(data_dic_temp['flux'][iexp][idxsub_ord_skysub]*dll_ord_skysub,axis=1)/np.nansum(dll_ord_skysub,axis=1) data_dic_temp['flux'][iexp][idxsub_ord_skysub] = (hdulist_skysub[1].data)[idx_ord_kept][idxsub_ord_skysub]/dll_ord_skysub mean_Fcorr_ord_skysub = np.nansum(data_dic_temp['flux'][iexp][idxsub_ord_skysub]*dll_ord_skysub,axis=1)/np.nansum(dll_ord_skysub,axis=1) skysub_corrfact = (mean_Fraw_ord_skysub/mean_Fcorr_ord_skysub)[:,None] data_dic_temp['flux'][iexp][idxsub_ord_skysub]*=skysub_corrfact err_raw[idxsub_ord_skysub] = ((hdulist_skysub[2].data)[idx_ord_kept][idxsub_ord_skysub]/dll_ord_skysub)*skysub_corrfact #Pixels quality # - unextracted pixels (on each side of the shorter blue orders and in between the blue/red chips in ESPRESSO) are considered as undefined qualdata_exp = (hdulist_dat[3].data)[idx_ord_kept] if (vis_path_skysub_exp is not None) and (gen_dic['fibB_corr'][inst][vis]!='all'): qualdata_exp[idxsub_ord_skysub] = (hdulist_skysub[3].data)[idx_ord_kept][idxsub_ord_skysub] cond_bad_exp |= (qualdata_exp == 16384) #Specific checks if inst=='ESPRESSO': #Removing pixels at the edge of the red detector # - usually too poorly defined even for bright stars for iord in [90,91]: isub_ord = np_where1D(np.array(idx_ord_kept)==iord) if len(isub_ord)>0: data_dic_temp['flux'][iexp][isub_ord,0:30] = np.nan qualdata_exp[isub_ord,0:30] = 1 #Check for bad wavelength grid # - some exposures have their wavelength grid badly defined (going in reverse order or too close) # this pixels are flagged as bad with QC=16384, but we still need to change their grids or the pipeline crashes (at least in pixels with pixels going in reverse order) dcen_bins = data_dic_temp['cen_bins'][iexp,:,1::]-data_dic_temp['cen_bins'][iexp,:,0:-1] cond_bad_grid = (dcen_bins<=0.) | (np.abs(dcen_bins)<1e-4) iord_bad = np_where1D( np.sum(cond_bad_grid,axis=1)>0 ) for iord in iord_bad: cond_bad_grid_ord = np.append(cond_bad_grid[iord],True) ipix_good = np_where1D( (~cond_bad_grid_ord) & (qualdata_exp[iord] == 0) ) if ipix_good[0]>0: npix_bad = ipix_good[0]-1 dcen_pix = dcen_bins[iord,ipix_good[0]] data_dic_temp['cen_bins'][iexp,iord,0:npix_bad] = (data_dic_temp['cen_bins'][iexp,iord,ipix_good[0]] - dcen_pix*(1.+np.arange(npix_bad)))[::-1] cond_bad_exp[iord,0:npix_bad-1] = True if ipix_good[-1]<data_inst[vis]['nspec']-1: npix_bad = data_inst[vis]['nspec']-1-ipix_good[-1] dcen_pix = dcen_bins[iord,ipix_good[-1]] data_dic_temp['cen_bins'][iexp,iord,ipix_good[-1]+1::] = (data_dic_temp['cen_bins'][iexp,iord,ipix_good[-1]] + dcen_pix*(1.+np.arange(npix_bad))) cond_bad_exp[iord,ipix_good[-1]+1::] = True #Calibration and detector noise profiles # - see weights_bin_prof() for details # - non sky-corrected data must be used if (vis in data_inst['gcal_blaze_vis']): #Retrieving counts per pixels # - blazed and unblazed pixels are undefined and null in the same pixels hdulist_dat_blaze =fits.open(vis_path_blaze_exp[iexp]) count_exp = (hdulist[1].data)[idx_ord_kept] count_blaze_exp = (hdulist_dat_blaze[1].data)[idx_ord_kept] err_count_blaze_exp = (hdulist_dat_blaze[2].data)[idx_ord_kept] cond_def = (~np.isnan(count_exp)) & (~cond_bad_exp) cond_gcal = cond_def & (count_exp!=0.) cond_pos_ct = count_blaze_exp>0. cond_def_pos = cond_def & cond_pos_ct cond_def_neg = cond_def & ~cond_pos_ct #Processing each order for iord in range(data_inst['nord']): #Defining calibration profile # - gcal(w,t,v) = 1/(dw(w)*bl(w,t,v))) # with bl(w,t,v) = N_meas[bl](w,t,v)/N_meas(w,t,v) # - negative values are set to undefined, so that they are replaced by the calibration model in the calc_gcal() routine # - blaze are undefined in pixels where counts are null idx_gcal_ord = np_where1D(cond_gcal[iord]) gcal_temp = 1./(dll[iord,idx_gcal_ord]*count_blaze_exp[iord,idx_gcal_ord]/count_exp[iord,idx_gcal_ord]) cond_gcal_pos = gcal_temp > 0. data_dic_temp['gcal'][iexp,iord,idx_gcal_ord[cond_gcal_pos]] = gcal_temp[cond_gcal_pos] #Defining detector noise # - error on raw pixel counts (ie, without deblazing and dll normalization) # - Edet_meas(w,t,v)^2 = EN_meas[bl](w,t,v)^2 - N_meas[bl](w,t,v) # - set to 0 where negative or undefined within definition range, after duplicating the latest defined value at the edges # - at pixels where counts are negative the DRS consider them null to define error tables, so that the detector noise is then retrieved as # Edet_meas(w,t,v)^2 = EN_meas[bl](w,t,v)^2 # - some pixels have large errors despite not having large fluxes; this may come from hot pixels subtracted from the extracted counts idx_def_pos = np_where1D(cond_def_pos[iord]) data_dic_temp['sdet2'][iexp,iord,idx_def_pos] = err_count_blaze_exp[iord,idx_def_pos]**2. - count_blaze_exp[iord,idx_def_pos] if idx_def_pos[0]>0:data_dic_temp['sdet2'][iexp,iord,0:idx_def_pos[0]] = data_dic_temp['sdet2'][iexp,iord,idx_def_pos[0]] if idx_def_pos[-1]<data_inst[vis]['nspec']-1:data_dic_temp['sdet2'][iexp,iord,idx_def_pos[-1]+1::] = data_dic_temp['sdet2'][iexp,iord,idx_def_pos[-1]] data_dic_temp['sdet2'][iexp,iord][np.isnan(data_dic_temp['sdet2'][iexp,iord])] = 0. data_dic_temp['sdet2'][iexp,iord][data_dic_temp['sdet2'][iexp,iord]<0.] = 0. data_dic_temp['sdet2'][iexp,iord,cond_def_neg[iord]] = err_count_blaze_exp[iord,cond_def_neg[iord]]**2. elif inst=='CORALIE': if iexp==0:print('TEMPORARY UNTIL UNIFIED FORMAT') #Bin centers PolyDeg = hdr['ESO DRS CAL TH DEG LL'] coeffLL = np.zeros(data_inst['nord'] * (PolyDeg + 1)) for i in range(0, data_inst['nord'] * (PolyDeg + 1)): coeffLL[i] = hdr['ESO DRS CAL TH COEFF LL' + str(i)] wave = np.zeros((data_inst['nord'], data_inst[vis]['nspec'])) xx = np.arange(data_inst[vis]['nspec']) for j in range(0, data_inst['nord']): for i in range(0, PolyDeg + 1): wave[j, :] = wave[j, :] + coeffLL[i + j * (PolyDeg + 1)] * xx ** i #Correction BERV wave*=gen_specdopshift(data_prop[inst][vis]['BERV'][iexp])*(1.+1.55e-8) data_dic_temp['cen_bins'][iexp] = wave[idx_ord_kept] dll = np.diff(wave, axis=1) dll = np.append(dll,dll[:,-1][:,None],axis=1) #Blaze hdulist_dat_blaze = fits.open(vis_path_blaze_exp[iexp]) blaze_exp = (hdulist_dat_blaze[0].data)[idx_ord_kept]/dll cond_pos_ct =(~np.isnan(blaze_exp)) & (blaze_exp>0.) #Flux spectrum # - S2D contains the blazed counts count_blazed_exp = (hdulist[0].data)[idx_ord_kept] cond_def = (~np.isnan(count_blazed_exp)) & cond_pos_ct data_dic_temp['flux'][iexp,cond_def] = count_blazed_exp[cond_def]/(dll[cond_def]*blaze_exp[cond_def]) #Noise spectrum on blazed counts ron = hdr['ESO CORA CCD RON'] ecount_blaze_exp2 = np.ones(count_blazed_exp.shape)*ron**2. cond_def_loc = ~np.isnan(count_blazed_exp) & (count_blazed_exp>0.) ecount_blaze_exp2[cond_def_loc]+=count_blazed_exp[cond_def_loc] #noise spectrum on flux err_raw = np.zeros(count_blazed_exp.shape)*np.nan err_raw[cond_def] = np.sqrt(ecount_blaze_exp2[cond_def])/(dll[cond_def]*blaze_exp[cond_def]) #Calibration and detector noise profiles if (vis in data_inst['gcal_blaze_vis']): for iord in range(data_inst['nord']): idx_gcal_ord = np_where1D(cond_pos_ct[iord]) data_dic_temp['gcal'][iexp,iord,idx_gcal_ord] = 1./(dll[iord,idx_gcal_ord]*blaze_exp[iord,idx_gcal_ord]) data_dic_temp['sdet2'][iexp] = np.ones(err_raw.shape)*ron elif inst=='CARMENES_VIS': if (vis_path_skysub_exp is not None):stop('No sky-corrected data available') #Bin centers # - dimension norder x nbins # - in vacuum data_dic_temp['cen_bins'][iexp] = (hdulist_dat['WAVE'].data)[idx_ord_kept] if gen_dic['sp_frame']=='air':data_dic_temp['cen_bins'][iexp]/=air_index(data_dic_temp['cen_bins'][iexp], t=15., p=760.) #Spectra data_dic_temp['flux'][iexp] = (hdulist_dat['SPEC'].data)[idx_ord_kept] err_raw = (hdulist_dat['SIG'].data)[idx_ord_kept] #BERV correction to shift exposures from the Earth laboratory to the solar System barycenter # - see gen_specdopshift(): # w_receiver = w_source * (1+ rv[source/receiver]/c)) # we correct for the BERV, the radial velocity of the Earth with respect to the solar System barycenter (source is the Earth, where signals are measured, and the receiver the barycenter) # w_solbar = w_Earth * (1+ rv[Earth/sun bary]/c)) # w_solbar = w_Earth * (1+ (BERV/c)) # we include a general+special relativity correction for the eccentricity and change of potential of Earth orbit around the Sun, averaged over a revolution (from C. Lovis) data_dic_temp['cen_bins'][iexp]*=gen_specdopshift(data_prop[inst][vis]['BERV'][iexp])*(1.+1.55e-8) elif inst=='EXPRES': if (vis_path_skysub_exp is not None):stop('No sky-corrected data available') #Bin centers # - dimension norder x nbins # - we use the chromatic-barycentric-corrected Excalibur wavelengths # Excalibur is a non-parametric, hierarchical method for wavelength calibration. These wavelengths suffer less from systematic or instrumental errors, but cover a smaller spectral range # at undefined Excalibur wavelengths we nonetheless set the regular barycentric-corrected wavelength, as ANTARESS requires fully defined spectral tables # since flux values are set to nan at undefined Excalibur wavelengths this has no impact on the results data_dic_temp['cen_bins'][iexp] = ((hdulist_dat[1].data)['bary_wavelength'])[idx_ord_kept] bary_excalibur = ((hdulist_dat[1].data)['bary_excalibur'])[idx_ord_kept] cond_def_wav = (hdulist_dat[1].data)['excalibur_mask'][idx_ord_kept] data_dic_temp['cen_bins'][iexp][cond_def_wav] = bary_excalibur[cond_def_wav] if gen_dic['sp_frame']=='air':data_dic_temp['cen_bins'][iexp]/=air_index(data_dic_temp['cen_bins'][iexp], t=15., p=760.) #Spectra # - the flux is not per unit of time, but corresponds to the number of photoelectrons measured during the exposure corrected for the blaze function data_dic_temp['flux'][iexp] = (hdulist_dat[1].data)['spectrum'][idx_ord_kept] err_raw = (hdulist_dat[1].data)['uncertainty'][idx_ord_kept] #Set flux to nan at undefined Excalibur wavelengths # - for consistency with the way pixels are usually defined # - excalibur_mask: gives a mask that excludes pixels without Excalibur wavelengths, which have been set to NaNs. These are pixels that do not fall between calibration lines in the order, and so are only on the edges of each order, or outside the range of LFC lines, as described above. data_dic_temp['flux'][iexp][~cond_def_wav]=np.nan #Pixels quality # - bad pixels are set to False originally # - pixel_mask: pixels with too low signal to return a proper extracted value. This encompasses largely pixels on the edges of orders qualdata_exp = np.zeros(data_dic_temp['flux'][iexp].shape) qualdata_exp[~(hdulist_dat[1].data)['pixel_mask'][idx_ord_kept]] = 1. #Telluric spectrum if ('tell' in data_dic_temp): data_dic_temp['tell'][iexp] = (hdulist_dat[1].data)['tellurics'][idx_ord_kept] else: stop('ERROR : spectra upload undefined for '+inst) #Set bad quality pixels to nan # - bad pixels have flag > 0 # - unless requested by the user those pixels are not set to undefined, because we cannot propagate detailed quality information when a profile is resampled / modified and leaving them undefined # will 'bleed' to adjacent pixels when resampling and will limit some calculations if gen_dic['bad2nan']:cond_bad_exp |= (qualdata_exp > 0.) if True in cond_bad_exp:data_dic_temp['flux'][iexp][cond_bad_exp] = np.nan #Set undefined errors to square-root of flux as first-order approximation # - errors will be redefined for local CCF using their continuum dispersion # - pixels with negative values are set to undefined, as there is no simple way to attribute errors to them if not gen_dic[inst][vis]['flag_err']: cond_pos_exp = data_dic_temp['flux'][iexp]>=0. data_dic_temp['flux'][iexp][~cond_pos_exp] = np.nan err_raw = np.zeros(data_inst[vis]['dim_exp'],dtype=float) err_raw[cond_pos_exp]=np.sqrt(data_dic_temp['flux'][iexp][cond_pos_exp]) #------------------------------ #Processing orders for iord in range(data_inst['nord']): #Set pixels with null flux on order sides to undefined # - those pixels contain no information and can be removed to speed calculations # - we do not set to undefined all null pixels in case some are within the spectra if ('spec' in data_inst['type']): cumsum = np.nancumsum(data_dic_temp['flux'][iexp][iord]) if (True in (cumsum>0.)): idx_null_left = np_where1D(cumsum>0.)[0]-1 if (idx_null_left>=0) and (np.nansum(data_dic_temp['flux'][iexp][iord][0:idx_null_left+1])==0.):data_dic_temp['flux'][iexp][iord][0:idx_null_left+1] = np.nan if (True in (cumsum<cumsum[-1])): idx_null_right = np_where1D(cumsum<cumsum[-1])[-1]+2 if (idx_null_right<data_inst[vis]['nspec']) and (np.nansum(data_dic_temp['flux'][iexp][iord][idx_null_right::])==0.):data_dic_temp['flux'][iexp][iord][idx_null_right::] = np.nan #Banded covariance matrix # - dimensions (nd+1,nsp) for a given (exposure,order) # - here we have to assume there are no correlations in the uploaded profile, and store its error table in the matrix diagonal data_dic_temp['cov'][iexp,iord] = (err_raw[iord]**2.)[None,:] #----------------------------------------------------- #Observational properties associated with each exposure #----------------------------------------------------- if not gen_dic['mock_data']: #Airmass (mean between start and end of exposure) # - AM = sec(z) = 1/cos(z) with z the zenith angle if inst in ['HARPN','CARMENES_VIS']:data_prop[inst][vis]['AM'][iexp]=hdr['AIRMASS'] elif inst=='CORALIE':data_prop[inst][vis]['AM'][iexp]=hdr['HIERARCH ESO OBS TARG AIRMASS'] elif inst=='EXPRES':data_prop[inst][vis]['AM'][iexp]=hdr['AIRMASS'] else: if inst in ['HARPS','SOPHIE','ESPRESSO','NIRPS_HA','NIRPS_HE']: if inst in ['SOPHIE','HARPS','NIRPS_HA','NIRPS_HE']:pre_txt='' elif inst=='ESPRESSO':pre_txt = tel_inst data_prop[inst][vis]['AM'][iexp]=0.5*(hdr['HIERARCH '+facil_inst+' TEL'+pre_txt+' AIRM START']+hdr['HIERARCH '+facil_inst+' TEL'+pre_txt+' AIRM START'] ) elif inst=='ESPRESSO_MR': data_prop[inst][vis]['AM_UT'][iexp]=[] for tel_inst_loc in ['TEL1','TEL2','TEL3','TEL4']: pre_txt='HIERARCH ESO '+tel_inst_loc+' AIRM' data_prop[inst][vis]['AM_UT'][iexp]+=[0.5*(hdr[pre_txt+' START']+hdr[pre_txt+' END'] )] data_prop[inst][vis]['AM'][iexp]=np.mean(data_prop[inst][vis]['AM_UT'][iexp]) #Integrated water vapor at zenith (mm) if inst=='ESPRESSO': data_prop[inst][vis]['IWV_AM'][iexp]=0.5*(hdr['HIERARCH '+facil_inst+' TEL'+tel_inst+' AMBI IWV START']+hdr['HIERARCH '+facil_inst+' TEL'+tel_inst+' AMBI IWV START'] ) #No radiometer available at La Silla, TNG and CAHA elif inst in ['CARMENES_VIS','HARPN','HARPS','SOPHIE','NIRPS_HA','NIRPS_HE']: data_prop[inst][vis]['IWV_AM'][iexp] = 1. #Temperature (K) if inst in ['HARPS','SOPHIE','ESPRESSO','NIRPS_HA','NIRPS_HE']: if inst in ['SOPHIE','HARPS','NIRPS_HA','NIRPS_HE']:pre_txt='' elif inst=='ESPRESSO':pre_txt = tel_inst data_prop[inst][vis]['TEMP'][iexp] = hdr['HIERARCH '+facil_inst+' TEL'+pre_txt+' AMBI TEMP'] + 273.15 elif inst=='CORALIE': data_prop[inst][vis]['TEMP'][iexp] = hdr['HIERARCH '+facil_inst+' OBS AMBI TEMP'] + 273.15 elif inst=='HARPN': #Validated by Rosario Cosentino (HARPSN instrument scientist) pre_txt='' data_prop[inst][vis]['TEMP'][iexp] = hdr['HIERARCH '+facil_inst+' METEO TEMP10M'] + 273.15 elif inst=='CARMENES_VIS': pre_txt='' data_prop[inst][vis]['TEMP'][iexp] = hdr['HIERARCH '+facil_inst+' GEN AMBI TEMPERATURE'] + 273.15 #Pressure (atm) # - pressure is in hPa (=mbar) and must be converted into atm unit (1 atm = 101325 Pa = 1013.25 hPa) if inst in ['HARPS','SOPHIE','ESPRESSO','NIRPS_HA','NIRPS_HE']: if inst in ['SOPHIE','HARPS','NIRPS_HA','NIRPS_HE']:pre_txt='' elif inst=='ESPRESSO':pre_txt = tel_inst data_prop[inst][vis]['PRESS'][iexp]=0.5*(hdr['HIERARCH '+facil_inst+' TEL'+pre_txt+' AMBI PRES START']+hdr['HIERARCH '+facil_inst+' TEL'+pre_txt+' AMBI PRES END'] )/1013.2501 elif inst=='CORALIE': data_prop[inst][vis]['PRESS'][iexp] = hdr['HIERARCH '+facil_inst+' OBS AMBI PRESS'] / 1013.2501 elif inst=='HARPN': #Validated by Rosario Cosentino (HARPSN instrument scientist) pre_txt='' data_prop[inst][vis]['PRESS'][iexp] = hdr['HIERARCH '+facil_inst+' METEO PRESSURE'] / 1013.2501 elif inst=='CARMENES_VIS': pre_txt='' data_prop[inst][vis]['PRESS'][iexp] = hdr['HIERARCH '+facil_inst+' GEN AMBI PRESSURE'] / 1013.2501 #Seeing (mean between start and end of exposure) # - not saved for HARPSN data if inst=='CARMENES_VIS': data_prop[inst][vis]['seeing'][iexp]=hdr['HIERARCH CAHA GEN AMBI SEEING'] elif inst=='HARPS': pre_txt='HIERARCH ESO TEL AMBI FWHM' if (pre_txt+' START' in hdr):data_prop[inst][vis]['seeing'][iexp]=0.5*(hdr[pre_txt+' START']+hdr[pre_txt+' END'] ) elif inst=='SOPHIE': data_prop[inst][vis]['seeing'][iexp]=hdr['HIERARCH OHP GUID SEEING'] elif inst=='ESPRESSO_MR': data_prop[inst][vis]['seeing_UT'][iexp]=[] for tel_inst_loc in ['TEL1','TEL2','TEL3','TEL4']: data_prop[inst][vis]['seeing_UT'][iexp]+=[hdr['HIERARCH ESO '+tel_inst_loc+' IA FWHM']] data_prop[inst][vis]['seeing'][iexp]=np.mean(data_prop[inst][vis]['seeing_UT'][iexp]) elif inst=='ESPRESSO': data_prop[inst][vis]['seeing'][iexp]=hdr['HIERARCH ESO TEL'+tel_inst+' IA FWHM'] #seeing corrected for airmass elif inst in ['NIRPS_HA','NIRPS_HE']: data_prop[inst][vis]['seeing'][iexp]=hdr['HIERARCH ESO INS2 AOS ATM SEEING'] #Shape of PSF # first and second values: size in X and Y # third value : average (quadratic) size # fourth value : angle if inst in ['CORALIE','SOPHIE','HARPS','ESPRESSO','ESPRESSO_MR','HARPN']: if inst=='HARPN': pre_txt='HIERARCH '+facil_inst+' AG ACQ' data_prop[inst][vis]['PSF_prop'][iexp,0]=hdr[pre_txt+' FWHMX'] data_prop[inst][vis]['PSF_prop'][iexp,1]=hdr[pre_txt+' FWHMY'] elif inst in ['ESPRESSO','ESPRESSO_MR']: pre_txt='HIERARCH '+facil_inst+' OCS ACQ3' if pre_txt+' FWHMX ARCSEC' in hdr: data_prop[inst][vis]['PSF_prop'][iexp,0]=hdr[pre_txt+' FWHMX ARCSEC'] data_prop[inst][vis]['PSF_prop'][iexp,1]=hdr[pre_txt+' FWHMY ARCSEC'] data_prop[inst][vis]['PSF_prop'][iexp,2]=np.sqrt(data_prop[inst][vis]['PSF_prop'][iexp,0]**2.+data_prop[inst][vis]['PSF_prop'][iexp,1]**2.) if data_prop[inst][vis]['PSF_prop'][iexp,0]==0.: data_prop[inst][vis]['PSF_prop'][iexp,3]=np.nan else: data_prop[inst][vis]['PSF_prop'][iexp,3]=np.arctan2(data_prop[inst][vis]['PSF_prop'][iexp,1],data_prop[inst][vis]['PSF_prop'][iexp,0])*180./np.pi #Coefficient of color correction for Earth atmosphere # - a polynomial is fitted to the ratio between the spectrum averaged over each order, and a stellar template # - the coefficient below give the minimum and maximum values of this ratio (they may thus correspond to any order) # a low value corresponds to a strong flux decrease if inst in ['SOPHIE','HARPS','ESPRESSO','ESPRESSO_MR','HARPN']: if inst in ['CORALIE']:reduc_txt='ESO DRS ' elif inst=='SOPHIE':reduc_txt='OHP DRS ' elif inst in ['HARPS','ESPRESSO','ESPRESSO_MR']:reduc_txt='ESO QC ' elif inst=='HARPN':reduc_txt='ESO QC ' if ('HIERARCH ESO QC FLUX CORR') in hdr else 'TNG QC ' data_prop[inst][vis]['colcorrmin'][iexp]=hdr['HIERARCH '+reduc_txt+'FLUX CORR MIN'] data_prop[inst][vis]['colcorrmax'][iexp]=hdr['HIERARCH '+reduc_txt+'FLUX CORR MAX'] #Polynomial for color correction # - retrieved for each order if inst in ['HARPS','ESPRESSO','HARPN']: for iord in range(data_inst['nord']): data_prop[inst][vis]['colcorr_ord'][iexp,iord]=hdr['HIERARCH '+reduc_txt+'ORDER'+str(iord+1)+' FLUX CORR'] #Saturation check if inst in ['HARPS','ESPRESSO','HARPN']: if 'HIERARCH '+reduc_txt+'SATURATION CHECK' in hdr:data_prop[inst][vis]['satur_check'][iexp]=hdr['HIERARCH '+reduc_txt+'SATURATION CHECK'] #ADC and PIEZO properties # - properties are sometimes not correctly saved in headers, and are not critical for this sequence if (inst in ['ESPRESSO']) and (gen_dic['sequence'] not in ['st_master_tseries']): ref_adc = {'ESO-VLT-U1':'','ESO-VLT-U2':'2','ESO-VLT-U3':'3','ESO-VLT-U4':'4'}[hdr['TELESCOP']] data_prop[inst][vis]['adc_prop'][iexp,0]=hdr['HIERARCH ESO INS'+ref_adc+' ADC1 POSANG'] data_prop[inst][vis]['adc_prop'][iexp,3]=hdr['HIERARCH ESO INS'+ref_adc+' ADC2 POSANG'] #Conversion from hhmmss.ssss to deg data_prop[inst][vis]['adc_prop'][iexp,1]=hdr['HIERARCH ESO INS'+ref_adc+' ADC1 RA'] data_prop[inst][vis]['adc_prop'][iexp,2]=hdr['HIERARCH ESO INS'+ref_adc+' ADC1 DEC'] data_prop[inst][vis]['adc_prop'][iexp,4]=hdr['HIERARCH ESO INS'+ref_adc+' ADC2 RA'] data_prop[inst][vis]['adc_prop'][iexp,5]=hdr['HIERARCH ESO INS'+ref_adc+' ADC2 DEC'] for ikey in [1,2,4,5]: str_split = str(data_prop[inst][vis]['adc_prop'][iexp,ikey]).split('.') h_key = float(str_split[0][-6:-4]) min_key = float(str_split[0][-4:-2]) sec_key = float(str_split[0][-2:]+'.'+str_split[1]) data_prop[inst][vis]['adc_prop'][iexp,ikey] = 15.*h_key+15.*min_key/60.+15.*sec_key/3600. #PIEZO properties data_prop[inst][vis]['piezo_prop'][iexp,0]=hdr['HIERARCH ESO INS'+ref_adc+' TILT1 VAL1'] data_prop[inst][vis]['piezo_prop'][iexp,1]=hdr['HIERARCH ESO INS'+ref_adc+' TILT1 VAL2'] data_prop[inst][vis]['piezo_prop'][iexp,2]=hdr['HIERARCH ESO INS'+ref_adc+' TILT2 VAL1'] data_prop[inst][vis]['piezo_prop'][iexp,3]=hdr['HIERARCH ESO INS'+ref_adc+' TILT2 VAL2'] #Guide star coordinates if (inst in ['ESPRESSO']) and (data_inst['type']=='spec2D') : data_prop[inst][vis]['guid_coord'][iexp,0]=hdr['HIERARCH ESO ADA'+tel_inst+' GUID RA'] data_prop[inst][vis]['guid_coord'][iexp,1]=hdr['HIERARCH ESO ADA'+tel_inst+' GUID DEC'] #Properties derived by the pipeline if ('rv_pip' in DI_data_inst[vis]): if inst in ['CORALIE']: pre_txt='HIERARCH ESO DRS CCF' if data_inst['type']=='CCF':hdr_loc = hdr else:hdr_loc =fits.open(glob.glob(vis_path_root+data_prop[inst][vis]['exp_filenames'][iexp].split('e2ds')[0]+'ccf*')[0])[0].header DI_data_inst[vis]['rv_pip'][iexp]=hdr_loc[pre_txt+' RVC'] #Baryc RV (drift corrected) (km/s) DI_data_inst[vis]['erv_pip'][iexp]=hdr_loc['HIERARCH ESO DRS DVRMS']/1e3 elif inst=='SOPHIE': pre_txt='HIERARCH OHP DRS CCF' DI_data_inst[vis]['rv_pip'][iexp]=hdr[pre_txt+' RV'] #Baryc RV (drift corrected ?? ) (km/s) DI_data_inst[vis]['erv_pip'][iexp]=hdr[pre_txt+' ERR'] elif inst in ['HARPS','ESPRESSO','ESPRESSO_MR','HARPN','NIRPS_HA','NIRPS_HE']: reduc_txt=facil_inst+' QC ' pre_txt='HIERARCH '+reduc_txt+'CCF' DI_data_inst[vis]['rv_pip'][iexp]=hdr[pre_txt+' RV'] #Baryc RV (drift corrected) (km/s) DI_data_inst[vis]['erv_pip'][iexp]=hdr[pre_txt+' RV ERROR'] if inst in ['CORALIE']: DI_data_inst[vis]['FWHM_pip'][iexp]=hdr_loc[pre_txt+' FWHM'] DI_data_inst[vis]['ctrst_pip'][iexp]=hdr_loc[pre_txt+' CONTRAST'] else: DI_data_inst[vis]['FWHM_pip'][iexp]=hdr[pre_txt+' FWHM'] DI_data_inst[vis]['ctrst_pip'][iexp]=hdr[pre_txt+' CONTRAST'] if inst in ['HARPS','SOPHIE','ESPRESSO','ESPRESSO_MR','HARPN','NIRPS_HA','NIRPS_HE']: DI_data_inst[vis]['eFWHM_pip'][iexp]=hdr[pre_txt+' FWHM ERROR'] DI_data_inst[vis]['ectrst_pip'][iexp]=hdr[pre_txt+' CONTRAST ERROR'] #SNRs in all orders # - SNRs tables are kept associated with original orders if inst in ['SOPHIE','CORALIE']: if inst in ['CORALIE']:pre_txt='SPE' elif inst=='SOPHIE':pre_txt='CAL' for iorder in range(data_inst['nord_ref']): data_prop[inst][vis]['SNRs'][iexp,iorder]=hdr['HIERARCH '+facil_inst+' DRS '+pre_txt+' EXT SN'+str(iorder)] elif inst=='CARMENES_VIS': for iorder in range(data_inst['nord_ref']): data_prop[inst][vis]['SNRs'][iexp,iorder]=hdr['HIERARCH CARACAL FOX SNR '+str(iorder)] elif inst in ['HARPS','ESPRESSO','ESPRESSO_MR','HARPN','NIRPS_HA','NIRPS_HE']: for iorder in range(data_inst['nord_ref']): data_prop[inst][vis]['SNRs'][iexp,iorder]=hdr['HIERARCH '+facil_inst+' QC ORDER'+str(iorder+1)+' SNR'] elif inst=='EXPRES': #Correspondance channel : original orders # 1 : 28 ; 2 : 37 ; 3 : 45 ; 4 : 51 ; 5 : 58 ; 6 : 63 ; 7 : 68 ; 8 : 72 SNR_exp = np.zeros(86,dtype=float)*np.nan for ichannel,iord in zip(range(1,9),[28,37,45,51,58,63,68,72]): SNR_exp[iord] = (hdulist[2].header)['HIERARCH channel'+str(ichannel)+'_SNR'] data_prop[inst][vis]['SNRs'][iexp] = SNR_exp #Blaze and wavelength calibration files if inst=='ESPRESSO': for ii in range(1,34): if hdr['HIERARCH ESO PRO REC1 CAL%i CATG'%(ii)] == 'BLAZE_A': data_prop[inst][vis]['BLAZE_A'][iexp]=hdr['HIERARCH ESO PRO REC1 CAL%i DATAMD5'%(ii)] elif hdr['HIERARCH ESO PRO REC1 CAL%i CATG'%(ii)] == 'BLAZE_B': data_prop[inst][vis]['BLAZE_B'][iexp]=hdr['HIERARCH ESO PRO REC1 CAL%i DATAMD5'%(ii)] elif hdr['HIERARCH ESO PRO REC1 CAL%i CATG'%(ii)] == 'WAVE_MATRIX_THAR_FP_A': data_prop[inst][vis]['WAVE_MATRIX_THAR_FP_A'][iexp]=hdr['HIERARCH ESO PRO REC1 CAL%i DATAMD5'%(ii)] elif hdr['HIERARCH ESO PRO REC1 CAL%i CATG'%(ii)] == 'WAVE_MATRIX_FP_THAR_B': data_prop[inst][vis]['WAVE_MATRIX_THAR_FP_B'][iexp]=hdr['HIERARCH ESO PRO REC1 CAL%i DATAMD5'%(ii)] #Telescope altitude angle (degres) if inst=='ESPRESSO': data_prop[inst][vis]['alt'][iexp]=hdr['HIERARCH ESO TEL'+tel_inst+' ALT'] #Telescope azimuth angle (degrees) if inst=='ESPRESSO': data_prop[inst][vis]['az'][iexp]=hdr['HIERARCH ESO TEL'+tel_inst+' AZ'] #Activity indexes if DACE_sp: dace_idexp = closest( DACE_idx['bjd'] ,coord_dic[inst][vis]['bjd'][iexp] ) for key in data_inst[vis]['act_idx']: data_prop[inst][vis][key][iexp,:] = [DACE_idx[key][dace_idexp],DACE_idx[key+'_err'][dace_idexp]] elif inst=='EXPRES': data_prop[inst][vis]['ha'][iexp,:] = [(hdulist[1].header)['HALPHA'],1.879E-2] data_prop[inst][vis]['s'][iexp,:] = [(hdulist[1].header)['S-VALUE'], 6.611E-3] ### End of current exposure ### End of exposures in visit #------------------------------------------------------------------------------------ #Scaling of errors if not gen_dic[inst][vis]['flag_err']:print(' Uncertainties are set to sqrt(F)') if (inst not in gen_dic['g_err']):gen_dic['g_err'][inst]=1. elif gen_dic['g_err'][inst]!=1.: print(' Uncertainties are scaled by sqrt('+str(gen_dic['g_err'][inst])+')') data_dic_temp['cov']*=gen_dic['g_err'][inst] #Sort exposures by increasing time # - remove exposures if requested, using their index after time-sorting if (inst in gen_dic['used_exp']) and (vis in gen_dic['used_exp'][inst]) and (len(gen_dic['used_exp'][inst][vis])>0): remove_exp=True print(' Removing '+str(data_inst[vis]['n_in_visit']-len(gen_dic['used_exp'][inst][vis]))+' exposures from the analysis') else:remove_exp=False w_sorted=coord_dic[inst][vis]['bjd'].argsort() for pl_loc in data_inst[vis]['studied_pl']: for key in ['ecl','cen_ph','st_ph','end_ph','ph_dur','rv_pl','v_pl']: coord_dic[inst][vis][pl_loc][key]=coord_dic[inst][vis][pl_loc][key][w_sorted] if remove_exp:coord_dic[inst][vis][pl_loc][key] = coord_dic[inst][vis][pl_loc][key][gen_dic['used_exp'][inst][vis]] for key in ['cen_pos','st_pos','end_pos']: coord_dic[inst][vis][pl_loc][key]=coord_dic[inst][vis][pl_loc][key][:,w_sorted] if remove_exp:coord_dic[inst][vis][pl_loc][key] = coord_dic[inst][vis][pl_loc][key][:,gen_dic['used_exp'][inst][vis]] for ar in data_inst[vis]['studied_ar']: for key in list(gen_dic['ar_coord_par'])+['is_visible']: coord_dic[inst][vis][ar][key]=coord_dic[inst][vis][ar][key][:,w_sorted] if remove_exp:coord_dic[inst][vis][ar][key] = coord_dic[inst][vis][ar][key][:,gen_dic['used_exp'][inst][vis]] for key in ['bjd','t_dur','RV_star_solCDM','RV_star_stelCDM','cen_ph_st','st_ph_st','end_ph_st']: coord_dic[inst][vis][key]=coord_dic[inst][vis][key][w_sorted] if remove_exp:coord_dic[inst][vis][key] = coord_dic[inst][vis][key][gen_dic['used_exp'][inst][vis]] for key in data_prop[inst][vis]: if type(data_prop[inst][vis][key]) not in [str,int]: if len(data_prop[inst][vis][key].shape)==1: data_prop[inst][vis][key]=data_prop[inst][vis][key][w_sorted] if remove_exp:data_prop[inst][vis][key] = data_prop[inst][vis][key][gen_dic['used_exp'][inst][vis]] else: data_prop[inst][vis][key]=data_prop[inst][vis][key][w_sorted,:] if remove_exp:data_prop[inst][vis][key] = data_prop[inst][vis][key][gen_dic['used_exp'][inst][vis],:] for key in data_dic_temp: data_dic_temp[key]= data_dic_temp[key][w_sorted] if remove_exp:data_dic_temp[key] = data_dic_temp[key][gen_dic['used_exp'][inst][vis]] for key in DI_data_inst[vis]: DI_data_inst[vis][key]=DI_data_inst[vis][key][w_sorted] if remove_exp:DI_data_inst[vis][key] = DI_data_inst[vis][key][gen_dic['used_exp'][inst][vis]] if remove_exp: n_in_visit = len(gen_dic['used_exp'][inst][vis]) data_inst[vis]['n_in_visit'] = n_in_visit data_inst[vis]['dim_all'][0] = n_in_visit data_inst[vis]['dim_sp'][0] = n_in_visit data_inst[vis]['dim_ord'][0] = n_in_visit #Check if exposures have same duration coord_dic[inst][vis]['t_dur_d'] = coord_dic[inst][vis]['t_dur']/(3600.*24.) if (coord_dic[inst][vis]['t_dur']==coord_dic[inst][vis]['t_dur'][0]).all():coord_dic[inst][vis]['cst_tdur']=True else:coord_dic[inst][vis]['cst_tdur']=False #Mean SNR # - used to weigh CCFs (over contributing orders) # - used to weigh spectra (over all orders) if (not gen_dic['mock_data']): if ('CCF' in data_inst['type']): if inst not in gen_dic['orders4ccf']:data_prop[inst][vis]['mean_SNR'] = np.nanmean(data_prop[inst][vis]['SNRs'],axis=1) else:data_prop[inst][vis]['mean_SNR'] = np.nanmean(data_prop[inst][vis]['SNRs'][:,gen_dic['orders4ccf'][inst]],axis=1) elif ('spec' in data_inst['type']): data_prop[inst][vis]['mean_SNR'] = np.nanmean(data_prop[inst][vis]['SNRs'],axis=1) #Identify exposures west/east of meridian # - celestial meridian is at azimuth = 180, with azimuth counted clockwise from the celestial north # exposures at az > 180 are to the west of the meridian, exposures at az < 180 to the east if (inst=='ESPRESSO') and (not gen_dic['mock_data']): cond_eastmer = data_prop[inst][vis]['az'] < 180. data_inst[vis]['idx_eastmer'] = np_where1D(cond_eastmer) data_inst[vis]['idx_westmer'] = np_where1D(~cond_eastmer) data_inst[vis]['idx_mer'] = closest(data_prop[inst][vis]['az'] , 180.) #Estimate telescope altitude (deg) at meridian min_bjd = coord_dic[inst][vis]['bjd'][0] max_bjd = coord_dic[inst][vis]['bjd'][-1] dbjd_HR = 1./(3600.*24.) nbjd_HR = round((max_bjd-min_bjd)/dbjd_HR) bjd_HR=min_bjd+dbjd_HR*np.arange(nbjd_HR) az_HR = CubicSpline(coord_dic[inst][vis]['bjd'],data_prop[inst][vis]['az'])(bjd_HR) idx_mer_HR = closest(az_HR,180.) alt_HR = CubicSpline(coord_dic[inst][vis]['bjd'],data_prop[inst][vis]['alt'])(bjd_HR) data_inst[vis]['alt_mer'] = alt_HR[idx_mer_HR] data_inst[vis]['z_mer'] = np.sin(alt_HR[idx_mer_HR]*np.pi/180.) #------------------------------------------------------------------------------------ #Undefined bins data_dic_temp['cond_def'] = ~np.isnan(data_dic_temp['flux']) #Definition of edge spectral table # - spectral bins must be continuous for the resampling routine of the pipeline to work # - bin width are included in ESPRESSO files, but to have exactly the same lower/upper boundary between successive bins we redefine them manually data_dic_temp['edge_bins'] = def_edge_tab(data_dic_temp['cen_bins']) #Lower/upper spectral edges # - over all exposures in each order low_bin_ord_def = np.min(data_dic_temp['edge_bins'][:,:,0],axis=0) high_bin_ord_def = np.max(data_dic_temp['edge_bins'][:,:,-1],axis=0) #Nominal common spectral grid nom_nspec = data_inst[vis]['nspec'] cen_bins_com = data_dic_temp['cen_bins'][0] edge_bins_com = data_dic_temp['edge_bins'][0] low_edge_bins_com = edge_bins_com[:,0:-1] high_edge_bins_com = edge_bins_com[:,1::] #Bin size of spectral edges in each order of the common grid dbin_low_ord = edge_bins_com[:,1]-edge_bins_com[:,0] dbin_high_ord = edge_bins_com[:,-1]-edge_bins_com[:,-2] #Common spectral tables # - the common table is set to the table of the first exposure, defined in the input rest frame # it is then extended on both side of each order, from the min/max over all exposures, to the maximum rv shift in all frames of the workflow # - we give the table the same structure as for spectra, with dimensions nord x nspec # to keep a constant number of bins per order we use the maximum number of extended bins over all orders if ('spec' in data_inst['type']): #Upper and lower extension # - from the lower / upper spectral edges in each order # - we use the approximation w_shifted = w*(1 - (rv_shift/c) ) to define the boundaries low_bin_ord_exten = low_bin_ord_def*(1. - (gen_dic['min_rv_shift']/c_light) ) high_bin_ord_exten = high_bin_ord_def*(1. - (gen_dic['max_rv_shift']/c_light) ) #Number of extended bins # - we extend with a regular bin size in ln(w) space # see details in "Defining 1D table for conversion from 2D spectra" dlnbin_low_ord = dbin_low_ord/low_edge_bins_com[:,0] dlnbin_high_ord = dbin_high_ord/high_edge_bins_com[:,-1] nlow_bin_exten = np.max(npint( np.floor( np.log(low_edge_bins_com[:,0]/low_bin_ord_exten)/np.log( dlnbin_low_ord + 1. ) ) ) ) nhigh_bin_exten = np.max(npint( np.floor( np.log(high_bin_ord_exten/high_edge_bins_com[:,-1])/np.log( dlnbin_high_ord + 1. ) ) ) ) #Extended edges low_edge_bins_exten = ((low_edge_bins_com[:,0]*( dlnbin_low_ord + 1. )**(-np.arange(nlow_bin_exten+1))[:,None]).T)[:,::-1] high_edge_bins_exten = (high_edge_bins_com[:,-1]*( dlnbin_high_ord + 1. )**np.arange(nhigh_bin_exten+1)[:,None]).T else: #Upper and lower extension low_bin_ord_exten = low_bin_ord_def - gen_dic['min_rv_shift'] high_bin_ord_exten = high_bin_ord_def - gen_dic['max_rv_shift'] #Maximum number of extended bins # - taking the maximum over all order, between the edges of the common table and the min/max extensions nlow_bin_exten = np.max(npint((low_edge_bins_com[:,0]-low_bin_ord_exten)/dbin_low_ord)+1) nhigh_bin_exten = np.max(npint((high_bin_ord_exten-high_edge_bins_com[:,-1])/dbin_high_ord)+1) #Extended edges low_edge_bins_exten = ((low_edge_bins_com[:,0] - dbin_low_ord*np.arange(nlow_bin_exten+1)[:,None]).T)[:,::-1] high_edge_bins_exten = (high_edge_bins_com[:,-1] + dbin_high_ord*np.arange(nhigh_bin_exten+1)[:,None]).T #Extended centers low_cen_bins_exten = 0.5*(low_edge_bins_exten[:,0:-1]+low_edge_bins_exten[:,1::]) high_cen_bins_exten = 0.5*(high_edge_bins_exten[:,0:-1]+high_edge_bins_exten[:,1::]) #Number of bins in extended grid exten_nspec = nlow_bin_exten + nom_nspec + nhigh_bin_exten #Final grid cen_bins_com_exten = np.zeros([data_inst['nord'],exten_nspec],float)*np.nan edge_bins_com_exten = np.zeros([data_inst['nord'],exten_nspec+1],float)*np.nan for iord in range(data_inst['nord']): cen_bins_com_exten[iord] = np.concatenate((low_cen_bins_exten[iord], cen_bins_com[iord], high_cen_bins_exten[iord])) edge_bins_com_exten[iord] = np.concatenate((low_edge_bins_exten[iord,0:-1], edge_bins_com[iord], high_edge_bins_exten[iord,1::])) #Storing data_inst[vis]['proc_com_data_paths'] = gen_dic['save_data_dir']+'Processed_data/'+inst+'_'+vis+'_com' data_com = {'cen_bins':cen_bins_com_exten,'edge_bins':edge_bins_com_exten} data_com['nspec'] = exten_nspec data_com['dim_ord'] = [n_in_visit,exten_nspec] data_com['dim_all'] = [n_in_visit,data_inst['nord'],exten_nspec] data_com['dim_exp'] = [data_inst['nord'],exten_nspec] data_com['dim_sp'] = [n_in_visit,data_inst['nord']] #Initializes the common table in the star rest frame # - it will be updated if disk-integrated profiles are aligned data_inst[vis]['proc_com_star_data_paths'] = deepcopy(data_inst[vis]['proc_com_data_paths']) #Exposures in current visit are defined on different tables if ('spec' in data_inst['type']) and (np.nanmax(np.abs(data_dic_temp['edge_bins']-data_dic_temp['edge_bins'][0]))>1e-6): data_inst[vis]['comm_sp_tab'] = False #Resampling data on wavelength table common to the visit if requested if gen_dic['comm_sp_tab'][inst]: print(' Resampling on common table') #Processing each exposure flux_resamp = np.zeros(data_com['dim_all'], dtype=float)*np.nan cov_resamp = np.zeros(data_com['dim_sp'], dtype=object) if ('tell' in data_dic_temp):tell_resamp = np.ones(data_com['dim_all'], dtype=float) if (vis in data_inst['gcal_blaze_vis']): gcal_blaze_resamp = np.zeros(data_com['dim_all'], dtype=float)*np.nan sdet2_resamp = np.zeros(data_com['dim_all'], dtype=float) for iexp in range(n_in_visit): for iord in range(data_inst['nord']): flux_resamp[iexp,iord],cov_resamp[iexp][iord] = bind.resampling(data_com['edge_bins'][iord], data_dic_temp['edge_bins'][iexp,iord], data_dic_temp['flux'][iexp,iord] , cov = data_dic_temp['cov'][iexp][iord], kind=gen_dic['resamp_mode']) #Resampling input telluric spectrum if ('tell' in data_dic_temp): tell_resamp[iexp,iord] = bind.resampling(data_com['edge_bins'][iord], data_dic_temp['edge_bins'][iexp,iord], data_dic_temp['tell'][iexp,iord], kind=gen_dic['resamp_mode']) #Resampling blazed-derived calibration and detector noise profiles if (vis in data_inst['gcal_blaze_vis']): gcal_blaze_resamp[iexp,iord] = bind.resampling(data_com['edge_bins'][iord], data_dic_temp['edge_bins'][iexp,iord], data_dic_temp['gcal'][iexp,iord], kind=gen_dic['resamp_mode']) sdet2_resamp[iexp,iord] = bind.resampling(data_com['edge_bins'][iord], data_dic_temp['edge_bins'][iexp,iord], data_dic_temp['sdet2'][iexp,iord], kind=gen_dic['resamp_mode']) #Overwrite exposure tables data_dic_temp['cen_bins'][:] = deepcopy(data_com['cen_bins']) data_dic_temp['edge_bins'][:] = deepcopy(data_com['edge_bins']) data_dic_temp['flux'] = flux_resamp data_dic_temp['cov'] = cov_resamp data_dic_temp['cond_def'] = ~np.isnan(data_dic_temp['flux']) if ('tell' in data_dic_temp):data_dic_temp['tell'] = tell_resamp if (vis in data_inst['gcal_blaze_vis']): data_dic_temp['gcal'] = gcal_blaze_resamp data_dic_temp['sdet2'] = sdet2_resamp #Set flag that all exposures in the visit are now defined on a common table data_inst[vis]['comm_sp_tab'] = True #Save instrument table common to all visits # - set to the common table of the first instrument visit if ivis==0: data_dic[inst]['proc_com_data_path'] = gen_dic['save_data_dir']+'Processed_data/'+inst+'_com' data_dic[inst]['proc_com_star_data_path'] = deepcopy(data_dic[inst]['proc_com_data_path']) data_dic[inst]['dim_exp'] = deepcopy(data_inst[vis]['dim_exp']) data_dic[inst]['nspec'] = deepcopy(data_inst[vis]['nspec']) data_dic[inst]['com_vis'] = vis data_com_inst = deepcopy(data_com) datasave_npz(data_dic[inst]['proc_com_data_path'],data_com_inst) #Check wether all visits of current instrument share a common table # - flag remains True if the common table of current visit is the same as the common table for the instrument # - flag is set to False if : # + exposures of current visit are still defined on individual table (= no common table even for the visit) # + exposures of current visit are defined on a common table, but it is different from the instrument table (only possible if current visit is not the first, used as reference) if (not data_inst[vis]['comm_sp_tab']) or ((ivis>0) and ((data_com_inst['nspec']!=data_com['nspec']) or (np.nanmax(np.abs(data_com_inst['edge_bins']-data_com['edge_bins']))>1e-6))): data_dic[inst]['comm_sp_tab'] = False #------------------------------------------------------------------------------------ #Masking condition if ('spec' in data_inst['type']) and (inst in gen_dic['masked_pix']) and (vis in gen_dic['masked_pix'][inst]): cond_mask_gen = True print(' Masking spectral ranges from the analysis') else:cond_mask_gen = False #Final processing of exposures data_inst[vis]['proc_DI_data_paths']=gen_dic['save_data_dir']+'Processed_data/'+inst+'_'+vis+'_' for iexp in range(n_in_visit): #Set masked pixels to nan # - masked orders and positions are assumed to be the same for all exposures, but positions have to be specific to a given order as orders can overlap if cond_mask_gen: mask_exp = gen_dic['masked_pix'][inst][vis]['exp_list'] if (len(mask_exp)==0) or (iexp in mask_exp): for iord in gen_dic['masked_pix'][inst][vis]['ord_list']: mask_bd = gen_dic['masked_pix'][inst][vis]['ord_list'][iord] cond_mask = np.zeros(data_inst[vis]['nspec'],dtype=bool) for bd_int in mask_bd: cond_mask |= (data_dic_temp['edge_bins'][iexp,iord,0:-1]>=bd_int[0]) & (data_dic_temp['edge_bins'][iexp,iord,1:]<=bd_int[1]) data_dic_temp['cond_def'][iexp ,iord , cond_mask ] = False data_dic_temp['flux'][iexp ,iord , cond_mask ] = np.nan #Telluric spectrum # - if defined with input data # - path is made specific to the exposure to be able to point from in-transit to global profiles without copying them to disk if ('tell' in data_dic_temp): data_inst[vis]['tell_DI_data_paths'][iexp] = data_inst[vis]['proc_DI_data_paths']+'tell_'+str(iexp) datasave_npz(data_inst[vis]['tell_DI_data_paths'][iexp],{'tell':data_dic_temp['tell'][iexp]}) #Blazed-derived calibration and detector noise profiles # - if defined with input data # - path is made specific to the exposure to be able to point from in-transit to global profiles without copying them to disk if (vis in data_inst['gcal_blaze_vis']): data_inst[vis]['sing_gcal_DI_data_paths'][iexp] = data_inst[vis]['proc_DI_data_paths']+'sing_gcal_'+str(iexp) data_gcal_save = { 'gcal':data_dic_temp['gcal'][iexp], 'sdet2':data_dic_temp['sdet2'][iexp]} datasave_npz(data_inst[vis]['sing_gcal_DI_data_paths'][iexp],data_gcal_save) #Saving path to initialized raw data # - saving data per exposure to prevent size issue with npz files data_exp = {key:data_dic_temp[key][iexp] for key in ['cen_bins','edge_bins','flux','cov','cond_def']} datasave_npz(data_inst[vis]['proc_DI_data_paths']+str(iexp),data_exp) #Check for empty orders idx_ord_empty = np_where1D(np.sum(data_dic_temp['cond_def'],axis=(0,2)) == 0.) if len(idx_ord_empty)>0: txt_str = '[' for iord in idx_ord_empty[0:-1]:txt_str+=str(iord)+',' txt_str+=str(idx_ord_empty[-1])+']' print(' Empty orders : '+txt_str) #Defined edges of common spectral table low_bins_def = deepcopy(data_dic_temp['edge_bins'][:,:,0:-1]) high_bins_def = deepcopy(data_dic_temp['edge_bins'][:,:,1::]) low_bins_def[~data_dic_temp['cond_def']] = np.nan high_bins_def[~data_dic_temp['cond_def']] = np.nan data_com['min_edge_ord'] = np.nanmin(low_bins_def,axis=(0,2)) data_com['max_edge_ord'] = np.nanmax(high_bins_def,axis=(0,2)) datasave_npz(data_inst[vis]['proc_com_data_paths'],data_com) #Saving useful keyword file if gen_dic['sav_keywords']: tup_save=(np.append('index',[str(iexp) for iexp in range(n_in_visit)]),) form_save=('%-6s',) for key,form in zip(['BLAZE_A','BLAZE_B','WAVE_MATRIX_THAR_FP_A','WAVE_MATRIX_THAR_FP_B'], ['%-35s' ,'%-35s' ,'%-35s' ,'%-35s' ]): if key in data_prop[inst][vis]: form_save+=(form,) col_loc = np.append(key,[str(data_prop[inst][vis][key][iexp]) for iexp in range(n_in_visit)]) tup_save+=(col_loc,) np.savetxt(gen_dic['save_data_dir']+'DIorig_prop/'+inst+'_'+vis+'_keywords.dat', np.column_stack(tup_save),fmt=form_save) #Saving dictionary elements defined within the routine np.savez_compressed(gen_dic['save_data_dir']+'Processed_data/Global/'+inst+'_'+vis,data_add=data_inst[vis],coord_add=coord_dic[inst][vis],data_prop_add=data_prop[inst][vis],gen_add=gen_dic[inst][vis],DI_data_add=DI_data_inst[vis],allow_pickle=True) ### End of visits #Saving general information np.savez_compressed(gen_dic['save_data_dir']+'Processed_data/Global/'+inst,gen_inst_add=gen_dic[inst],data_inst_add=data_dic[inst],allow_pickle=True) #Retrieving data else: check_data({'path':gen_dic['save_data_dir']+'Processed_data/Global/'+inst}) data_load = np.load(gen_dic['save_data_dir']+'Processed_data/Global/'+inst+'.npz',allow_pickle=True) data_dic[inst]=data_load['data_inst_add'].item() gen_dic[inst]= data_load['gen_inst_add'].item() #Remove visits # - this allows removing visits after processing all of them once via gen_dic['calc_proc_data'] for vis in gen_dic['unused_visits'][inst]: if vis in data_dic[inst]['visit_list']:data_dic[inst]['visit_list'].remove(vis) #Retrieve/initialize visit dictionaries for vis in data_dic[inst]['visit_list']: data_load = np.load(gen_dic['save_data_dir']+'Processed_data/Global/'+inst+'_'+vis+'.npz',allow_pickle=True) data_dic[inst][vis]=data_load['data_add'].item() data_dic[inst][vis]['proc_DI_data_paths'] = gen_dic['save_data_dir']+'Processed_data/'+inst+'_'+vis+'_' #temporary redef; retrieved in data_load['data_add'].item() no ? data_dic[inst][vis]['proc_com_data_paths'] = gen_dic['save_data_dir']+'Processed_data/'+inst+'_'+vis+'_com' #temporary redef; retrieved in data_load['data_add'].item() no ? coord_dic[inst][vis]=data_load['coord_add'].item() data_prop[inst][vis]=data_load['data_prop_add'].item() gen_dic[inst][vis]=data_load['gen_add'].item() data_dic['DI'][inst][vis]=data_load['DI_data_add'].item() theo_dic[inst][vis]={} data_dic['Atm'][inst][vis]={} #Default transit model if (inst not in data_dic['DI']['transit_prop']):data_dic['DI']['transit_prop'][inst] = {} #Duplicate chromatic properties so that they are not overwritten by conversions data_dic[inst]['system_prop'] = deepcopy(data_dic['DI']['system_prop']) data_dic[inst]['ar_prop'] = deepcopy(data_dic['DI']['ar_prop']) #Final processing if len(data_dic[inst]['visit_list'])>1: if (not data_dic[inst]['comm_sp_tab']):print(' Visits do not share a common spectral table') else:print(' All visits share a common spectral table') ar_check = False for vis in data_dic[inst]['visit_list']: if (vis in data_dic[inst]['single_night']): print(' Processing visit '+vis) if not gen_dic['mock_data']: if vis in data_dic[inst]['dates']:print(' Date (night start) : '+data_dic[inst]['dates'][vis]) if vis in data_dic[inst]['midpoints']:print(' Visit midpoint: '+data_dic[inst]['midpoints'][vis]) else: print(' Processing multi-epoch visit '+vis) data_vis = data_dic[inst][vis] coord_vis = coord_dic[inst][vis] gen_vis = gen_dic[inst][vis] if (not data_vis['comm_sp_tab']):print(' Exposures do not share a common spectral table') else:print(' All exposures share a common spectral table') #Keeping track of original dimensions data_vis['nspec_ref'] = deepcopy(data_vis['nspec']) #Initialize all exposures as being defined data_dic['DI'][inst][vis]['idx_def'] = np.arange(data_vis['n_in_visit'],dtype=int) #------------------------------------------------------------------------------------ #Default systemic velocity if (inst not in data_dic['DI']['sysvel']):data_dic['DI']['sysvel'][inst] = {} if (vis not in data_dic['DI']['sysvel'][inst]): data_dic['DI']['sysvel'][inst][vis] = system_param['star']['sysvel'] print(' WARNING : sysvel set to default '+str(system_param['star']['sysvel'])+' km/s') #Min/max BERV if ('spec' in data_dic[inst]['type']): data_vis['min_BERV'] = np.nanmin(data_prop[inst][vis]['BERV']) data_vis['max_BERV'] = np.nanmax(data_prop[inst][vis]['BERV']) #Define CCF resolution, range, and sysvel here so that spectral data needs not be uploaded again after processing if gen_dic['CCF_from_sp'] or gen_dic['Intr_CCF']: #Velocity table in input frame if gen_dic['dRV'] is None:dvelccf= return_pix_size(inst) else:dvelccf= gen_dic['dRV'] n_vel = int((gen_dic['end_RV']-gen_dic['start_RV'])/dvelccf)+1 data_vis['nvel'] = n_vel data_vis['velccf'] = gen_dic['start_RV']+dvelccf*np.arange(n_vel) data_vis['edge_velccf'] = np.append(data_vis['velccf']-0.5*dvelccf,data_vis['velccf'][-1]+0.5*dvelccf) #Common velocity table in original frame nlow_bin_exten = int( gen_dic['min_rv_shift']/dvelccf)+1 nhigh_bin_exten = int(-gen_dic['max_rv_shift']/dvelccf)+1 low_edge_bins_exten = (data_vis['edge_velccf'][0] - dvelccf*np.arange(nlow_bin_exten+1))[::-1] high_edge_bins_exten = (data_vis['edge_velccf'][-1] + dvelccf*np.arange(nhigh_bin_exten+1)) low_cen_bins_exten = 0.5*(low_edge_bins_exten[0:-1]+low_edge_bins_exten[1::]) high_cen_bins_exten = 0.5*(high_edge_bins_exten[0:-1]+high_edge_bins_exten[1::]) data_vis['velccf_com'] = np.concatenate((low_cen_bins_exten, data_vis['velccf'], high_cen_bins_exten)) data_vis['edge_velccf_com'] = np.concatenate((low_edge_bins_exten[0:-1], data_vis['edge_velccf'], high_edge_bins_exten[1::])) #Shifting the CCFs velocity table by RV(CDM_star/CDM_sun) for key in ['velccf','edge_velccf','velccf_com','edge_velccf_com']:data_vis[key+'_star']=data_vis[key]-data_dic['DI']['sysvel'][inst][vis] #Automatic continuum and fit range # - done here rather than in the 'calc_proc_data' condition so that ranges can be defined for already-processed observed or mock datasets, even if the analysis modules were not activated # at the time of processing # - called if fit of single or joined profiles are requested, and if the continuum of differential profiles is not defined in CCF mode for key in gen_dic['type_list']: if ((key=='Diff') and gen_dic['diff_data'] and (data_dic[key]['type'][inst]=='CCF')) or gen_dic['fit_'+key+'_gen'] or gen_dic['fit_'+key+'Prof']: autom_cont = True if (inst not in data_dic[key]['cont_range']) else False autom_fit = True if (key!='Diff') and ((inst not in data_dic[key]['fit_range']) or (vis not in data_dic[key]['fit_range'][inst])) else False if autom_cont or autom_fit: #Path to data used for automatic definition # - if spectra are converted into CCFs, we cannot use the original spectral data if (data_dic[key]['type'][inst]=='CCF') and (data_dic[inst]['type']=='spec2D'): data_path = gen_dic['save_data_dir']+key+'_data/CCFfromSpec/'+gen_dic['add_txt_path'][key]+'/'+inst+'_'+vis+'_' nspec = data_vis['nvel'] else: data_path = data_vis['proc_DI_data_paths'] nspec = data_vis['nspec'] #Processing CCF order or selected spectral order if (data_dic[key]['type'][inst]=='CCF'):iord_sel = 0 else:iord_sel = data_dic[key]['fit_order'][inst] cen_bins_all = np.zeros([0,nspec],dtype=float) edge_bins_all = np.zeros([0,nspec+1],dtype=float) flux_all = np.zeros([0,nspec],dtype=float) for iexp in range(data_vis['n_in_visit']): data_exp = dataload_npz(data_path+str(iexp)) cen_bins_all=np.append(cen_bins_all,[data_exp['cen_bins'][iord_sel]],axis=0) #dimension nexp x nspec edge_bins_all=np.append(edge_bins_all,[data_exp['edge_bins'][iord_sel]],axis=0) #dimension nexp x nspec+1 flux_all=np.append(flux_all,[data_exp['flux'][iord_sel]],axis=0) #dimension nexp x nspec #Range space # - continuum is defined in RV space if data is in CCF mode, and in wavelength space if data is in spectral mode # - fit range is defined in RV space if data is in CCF mode or in spectral mode and fitted with an analytical model on a single line, and in wavelength space otherwise if autom_fit: if ('spec' in data_dic[key]['type'][inst]) and (data_dic[key]['fit_prof']['mode']=='ana'): if data_dic[key]['fit_prof']['line_trans'] is None:stop('Define "line_trans" to fit spectral data with "mode = ana"') spec2rv = True else:spec2rv = False else:spec2rv = False #Condition to define line-specific boundaries if autom_fit or (autom_cont and ((data_dic[key]['type'][inst]=='CCF') or spec2rv)): if (data_dic[key]['type'][inst]=='CCF'): RVcen_bins = cen_bins_all RVedge_bins = edge_bins_all else: RVcen_bins = c_light*((cen_bins_all/data_dic[key]['line_trans']) - 1.) RVedge_bins = c_light*((edge_bins_all/data_dic[key]['line_trans']) - 1.) #Estimate of systemic velocity for disk-integrated profiles # - as median of the RV corresponding to the CCF minima over the visit # - differential, intrinsic and atmospheric profiles are analyzed in the star rest frame cond_check = (RVcen_bins[0]>-500.) & (RVcen_bins[0]<500.) idx_min_all = np.argmin(flux_all[:,cond_check],axis=1) #dimension nexp (index of minimum over nspec_check for each exposure) sys_RV = np.nanmedian(np.take_along_axis(RVcen_bins[:,cond_check],idx_min_all[:,None],axis=1)) if key=='DI':cen_RV = sys_RV else:cen_RV = 0. #Minimum/maximum velocity of the CCF range in the input rest frame min_bin = np.max(np.nanmin(RVedge_bins)) max_bin = np.min(np.nanmax(RVedge_bins)) if key in ['Diff','Intr','Atm']: min_bin-=sys_RV max_bin-=sys_RV #Excluded range for disk-integrated, differential, and intrinsic profiles # - we assume +-3 vsini accounts for both rotational and thermal broadening of DI profiles, and for rotational shift + thermal broadening if key in ['DI','Diff','Intr']: min_exc = np.max([cen_RV - 3.*system_param['star']['vsini'] -5.,min_bin+5.]) max_exc = np.min([cen_RV + 3.*system_param['star']['vsini'] +5.,max_bin-5.]) #Excluded planet range # - we account for a width of 20km/s over the range of studied orbital velocities, for all planets considered for atmospheric signals if (key=='Atm') or (len(data_dic['Atm']['no_plrange'])>0): min_pl_RV = 1e100 max_pl_RV = -1e100 if data_dic['Atm']['pl_atm_sign']=='Absorption':idx_atm = gen_dic[inst][vis]['idx_in'] elif data_dic['Atm']['pl_atm_sign']=='Emission': idx_atm = data_dic[inst][vis]['n_in_tr'] for pl_atm in data_dic[inst][vis]['pl_with_atm']: min_pl_RV = np.min([min_pl_RV,np.min(coord_dic[inst][vis][pl_atm]['rv_pl'][idx_atm])]) max_pl_RV = np.max([max_pl_RV,np.max(coord_dic[inst][vis][pl_atm]['rv_pl'][idx_atm])]) #Exclusion in atmospheric profiles if key=='Atm': min_exc = np.min(min_pl_RV -20.,min_bin+5.) max_exc = np.max(max_pl_RV +20.,max_bin-5.) else: #Exclusion in other profiles if (key=='DI') and ('DI' in key_loc for key_loc in data_dic['Atm']['no_plrange']): min_exc = np.min((min_exc,min_pl_RV -20.+sys_RV,min_bin+5.)) max_exc = np.max((max_exc,max_pl_RV +20.+sys_RV,max_bin-5.)) elif ((key=='Diff') and ('Diff' in key_loc for key_loc in data_dic['Atm']['no_plrange'])) or ((key=='Intr') and ('Intr' in key_loc for key_loc in data_dic['Atm']['no_plrange'])): min_exc = np.min((min_exc,min_pl_RV -20.,min_bin+5.)) max_exc = np.max((max_exc,max_pl_RV +20.,max_bin-5.)) #Outer boundaries of continuum and fitted ranges in RV space min_contfit = np.max([min_bin, min_exc - 50.]) max_contfit = np.min([max_bin, max_exc + 50.]) #Continuum range # - in spectral mode, the continuum is set to the full line range if a fit is requested in spectral mode with an analytical model if autom_cont: if (data_dic[key]['type'][inst]=='CCF'): data_dic[key]['cont_range'][inst] = {iord_sel:[[min_contfit,min_exc],[max_exc,max_contfit]]} else: if spec2rv: min_contrange = data_dic[key]['line_trans']*gen_specdopshift(min_contfit) min_exrange = data_dic[key]['line_trans']*gen_specdopshift(min_exc) max_contrange = data_dic[key]['line_trans']*gen_specdopshift(max_contfit) max_exrange = data_dic[key]['line_trans']*gen_specdopshift(max_exc) data_dic[key]['cont_range'][inst] = {iord_sel:[[min_contrange,min_exrange],[max_exrange,max_contrange]]} else: data_dic[key]['cont_range'][inst] = {iord_sel:[[edge_bins_all[0],edge_bins_all[-1]]]} #Fitting range if autom_fit: if inst not in data_dic[key]['fit_range']:data_dic[key]['fit_range'][inst]={} if spec2rv: min_fitrange = data_dic[key]['line_trans']*gen_specdopshift(min_contfit) max_fitrange = data_dic[key]['line_trans']*gen_specdopshift(max_contfit) else: min_fitrange = min_contfit max_fitrange = max_contfit data_dic[key]['fit_range'][inst][vis] = [[min_fitrange,max_fitrange]] #------------------------------------------------------------------------------------ #Keplerian motion relative to the stellar CDM and the Sun (km/s) for iexp in range(data_dic[inst][vis]['n_in_visit']): coord_vis['RV_star_stelCDM'][iexp],coord_vis['RV_star_solCDM'][iexp] = calc_rv_star(coord_dic,inst,vis,system_param,gen_dic,coord_dic[inst][vis]['bjd'][iexp],coord_dic[inst][vis]['t_dur'][iexp],data_dic['DI']['sysvel'][inst][vis]) #Using sky-corrected data for current visit if (inst in gen_dic['fibB_corr']) and (vis in gen_dic['fibB_corr'][inst]):print(' ',vis,'is sky-corrected') #Indices of in/out exposures in global tables print(' '+str(data_vis['n_in_visit'])+' exposures') if ('in' in data_dic['DI']['idx_ecl']) and (inst in data_dic['DI']['idx_ecl']['in']) and (vis in data_dic['DI']['idx_ecl']['in'][inst]): for iexp in data_dic['DI']['idx_ecl']['in'][inst][vis]: for pl_loc in data_vis['studied_pl']:coord_vis[pl_loc]['ecl'][iexp] = 3*np.sign(coord_vis[pl_loc]['ecl'][iexp]) if ('out' in data_dic['DI']['idx_ecl']) and (inst in data_dic['DI']['idx_ecl']['out']) and (vis in data_dic['DI']['idx_ecl']['out'][inst]): for iexp in data_dic['DI']['idx_ecl']['out'][inst][vis]: for pl_loc in data_vis['studied_pl']:coord_vis[pl_loc]['ecl'][iexp] = 1*np.sign(coord_vis[pl_loc]['ecl'][iexp]) cond_in = np.zeros(data_vis['n_in_visit'],dtype=bool) cond_pre = np.ones(data_vis['n_in_visit'],dtype=bool) cond_post = np.ones(data_vis['n_in_visit'],dtype=bool) for pl_loc in data_vis['studied_pl']: cond_in|=(np.abs(coord_vis[pl_loc]['ecl'])>1) cond_pre&=(coord_vis[pl_loc]['ecl']==-1.) cond_post&=(coord_vis[pl_loc]['ecl']==1.) gen_vis['idx_in']=np_where1D(cond_in) gen_vis['idx_out']=np_where1D(~cond_in) gen_vis['idx_pretr']=np_where1D(cond_pre) gen_vis['idx_posttr']=np_where1D(cond_post) data_vis['n_in_tr'] = len(gen_vis['idx_in']) data_vis['n_out_tr'] = len(gen_vis['idx_out']) gen_vis['idx_exp2in'] = np.zeros(data_vis['n_in_visit'],dtype=int)-1 gen_vis['idx_exp2in'][gen_vis['idx_in']]=np.arange(data_vis['n_in_tr']) gen_vis['idx_in2exp'] = np.arange(data_vis['n_in_visit'],dtype=int)[gen_vis['idx_in']] #Duplicate chromatic properties so that they are not overwritten by conversions data_vis['system_prop'] = deepcopy(data_dic['DI']['system_prop']) data_vis['ar_prop'] = deepcopy(data_dic['DI']['ar_prop']) if len(data_vis['studied_ar'])>0:ar_check = True if (len(gen_dic['studied_ar'])>0) and (not ar_check):stop('ERROR: no active regions are associated with any visit but are requested in simulation - reset workflow with "gen_dic["calc_proc_data"]"') if (len(gen_dic['studied_ar'])==0) and (ar_check):stop('ERROR: active regions are associated with visits but not requested in simulation - reset workflow with "gen_dic["calc_proc_data"]"') #Total number of visits for current instrument gen_dic[inst]['n_visits'] = len(data_dic[inst]['visit_list']) data_dic[inst]['n_visits_inst']=len(data_dic[inst]['visit_list']) #Defining 1D table for conversion from 2D spectra # - uniformly spaced in ln(w) # - we impose d[ln(w)] = dw/w = ( w(i+1)-w(i) ) / w(i) # thus the table can be built with: # d[ln(w)]*w(i) = ( w(i+1) - w(i) ) # w(i+1) = ( d[ln(w)] + 1 )*w(i) # w(i+1) = ( d[ln(w)] + 1 )^2*w(i-1) # w(i) = ( d[ln(w)] + 1 )^i*w(0) # the table ends at the largest n so that : # w(n) <= wmax < w(n+1) # ( d[ln(w)] + 1 )^n <= wmax/w(0) < ( d[ln(w)] + 1 )^(n+1) # n*ln( d[ln(w)] + 1 ) <= ln(wmax/w(0)) < (n+1)*ln( d[ln(w)] + 1 ) # n <= ln(wmax/w(0))/ln( d[ln(w)] + 1 ) < (n+1) # n = E( ln(wmax/w(0))/ln( d[ln(w)] + 1 ) ) # - we do not set an automatic definition because tables for intrinsic and atmospheric profiles will be shifted compared to the original tables if gen_dic['spec_1D']: def def_1D_bins(prop_dic): spec_1D_prop = prop_dic['spec_1D_prop'][inst] spec_1D_prop['nspec'] = int( np.floor( np.log(spec_1D_prop['w_end']/spec_1D_prop['w_st'])/np.log( spec_1D_prop['dlnw'] + 1. ) ) ) spec_1D_prop['edge_bins'] = spec_1D_prop['w_st']*( spec_1D_prop['dlnw'] + 1. )**np.arange(spec_1D_prop['nspec']+1) spec_1D_prop['cen_bins'] = 0.5*(spec_1D_prop['edge_bins'][0:-1]+spec_1D_prop['edge_bins'][1::]) return None for key in ['DI','Intr','Atm']: if gen_dic['spec_1D_'+key]: #Default 1D table # - tables are uniformely spaced in ln(w), with d[ln(w)] = dw/w = drv/c # we use the instrument pixel size defined in rv space to set a default value if (inst not in data_dic[key]['spec_1D_prop']): dlnw_inst = return_pix_size(inst)/c_light w_st_inst_all = {'CORALIE':3500.,'ESPRESSO':3000.,'HARPN':3000.,'HARPS':3000.,'CARMENES_VIS':4500., 'NIRPS_HA':9500.} if (inst not in w_st_inst_all):print('ERROR : define w_st for "'+inst+'" in ANTARESS_main.py > init_inst()') w_end_inst_all = {'CORALIE':7500.,'ESPRESSO':8500.,'HARPN':7500.,'HARPS':7500.,'CARMENES_VIS':11500.,'NIRPS_HA':21000.} if (inst not in w_end_inst_all):print('ERROR : define w_end for "'+inst+'" in ANTARESS_main.py > init_inst()') data_dic[key]['spec_1D_prop'][inst] = {'dlnw':dlnw_inst,'w_st':w_st_inst_all[inst],'w_end':w_end_inst_all[inst]} print(' Automatic definition of 1D '+gen_dic['type_name'][key]+' table : dlnw = '+"{0:.3e}".format(dlnw_inst)+ ' ; w_st = '+"{0:.3f}".format(w_st_inst_all[inst])+' A ; w_end = '+"{0:.3f}".format(w_end_inst_all[inst])+' A.') #Defining table def_1D_bins(data_dic[key]) return None
[docs] def init_vis(data_prop,data_dic,vis,coord_dic,inst,system_param,gen_dic): r"""**Initialization: visit** Initializes visit-specific fields for the workflow. Args: TBD Returns: None """ #Reset data mode to input data_dic[inst]['nord'] = deepcopy(data_dic[inst]['nord_spec']) #Dictionaries initialization data_vis=data_dic[inst][vis] coord_vis = coord_dic[inst][vis] data_dic['Diff'][inst][vis]={} gen_vis = gen_dic[inst][vis] data_dic['Intr'][inst][vis]={} #Current rest frame data_dic['DI'][inst][vis]['rest_frame'] = 'input' #Current data data_vis['raw_exp_data_paths'] = deepcopy(data_vis['proc_DI_data_paths']) #Indices of in/out exposures in global tables print(' > '+str(data_vis['n_in_visit'])+' exposures') print(' ',data_vis['n_in_tr'],'in-transit') txt_loc = ' '+str(data_vis['n_out_tr'])+' out-of-transit' if data_vis['n_in_tr']>0:txt_loc+=' ('+str(len(gen_vis['idx_pretr']))+' pre / '+str(len(gen_vis['idx_posttr']))+' post)' print(txt_loc) data_vis['dim_in'] = [data_vis['n_in_tr']]+data_vis['dim_exp'] #Definition of CCF continuum, and ranges contributing to master out # - the continuum of a given CCFs account for the ranges selected as input, the planetary exclusion ranges, and the bins defined in the CCFs # - the final continuum is common to all CCFs # if at least one bin is undefined in one exposure it is undefined in the continuum plAtm_vis = data_dic['Atm'][inst][vis] #Generic definition of planetary range exclusion plAtm_vis['exclu_range_input']={'CCF':{},'spec':{}} plAtm_vis['exclu_range_star']={'CCF':{},'spec':{}} if (data_dic['Atm']['exc_plrange']): #Exposures to be processed # - ranges are excluded from in-transit data only if left undefined as input, since most of the time only in-transit absorption signals are analyzed # beware if an emission or absorption signal is present outside of the transit, in this case out-of-transit should be included manually if (inst in data_dic['Atm']['iexp_no_plrange']) and (vis in data_dic['Atm']['iexp_no_plrange'][inst]):plAtm_vis['iexp_no_plrange'] = data_dic['Atm']['iexp_no_plrange'][inst][vis] else:plAtm_vis['iexp_no_plrange'] = gen_vis['idx_in'] iexp_no_plrange = plAtm_vis['iexp_no_plrange'] #Planetary range in star and input rest frame # - we define here the lower/upper wavelength boundaries of excluded ranges for each planetary mask line in each requested exposure, based on a common rv range excluded in the planet rest frame # - the excluded range 'plrange' is defined in the planet rest frame, ie as RV(plrange/pl) # - masked profiles are either the input disk-integrated profiles, defined in their original sun barycentric rest frame (over RV(M/star) + RV(star/CDM_sun) ) or # processed profiles defined in the star rest frame (over RV(M/star)) # - tables have dimension (nline,2,nexp) in spectral mode, (2,nexp) in rv mode # - if input data is in spectral mode, ranges are defined in both spectral and rv space so that they are available in case of CCF conversion # if input data is in CCF mode, ranges are only defined in rv plrange_star={} for pl_loc in data_vis['studied_pl']: #Planetary range shifted from planet (source) to star (receiver) rest frame # - we define the range excluded in the planet rest frame as rv(atom/pl), where we consider an absorbing or emitting atom in the planet as the source # see gen_specdopshift(): # w_receiver = w_source * (1+ rv[s/r]/c)) # w_star = w_pl * (1+ (rv[atom/pl]/c)) * (1+ (rv[pl/star]/c)) plrange_star[pl_loc] = np.vstack((np.repeat(1e10,data_vis['n_in_visit']),np.repeat(-1e10,data_vis['n_in_visit']))) plrange_star[pl_loc][:,iexp_no_plrange] = data_dic['Atm']['plrange'][:,None] + coord_vis[pl_loc]['rv_pl'][iexp_no_plrange] if (True in np.isnan(coord_vis[pl_loc]['rv_pl'])):stop(' Run gen_dic["calc_proc_data"] again to calculate "rv_pl"') plAtm_vis['exclu_range_star']['CCF'][pl_loc] = plrange_star[pl_loc] #Planetary range shifted from planet (source) to solar (receiver) barycentric rest frame # w_solbar = w_star * (1+ rv[star/starbar]/c)) * (1+ rv[starbar/solbar]/c)) # = w_pl * (1+ (rv[atom/pl]/c)) * (1+ (rv[pl/star]/c)) * (1+ (rv[star/starbar]/c)) * (1+ (rv[starbar/solbar]/c)) plAtm_vis['exclu_range_input']['CCF'][pl_loc] = plrange_star[pl_loc] + coord_vis['RV_star_solCDM'][None,:] #Shifted spectral ranges if ('spec' in data_dic[inst]['type']): plAtm_vis['exclu_range_star']['spec'][pl_loc] = np.zeros([len(data_dic['Atm']['CCF_mask_wav']),2,data_vis['n_in_visit']])*np.nan plAtm_vis['exclu_range_star']['spec'][pl_loc][:,:,iexp_no_plrange] = data_dic['Atm']['CCF_mask_wav'][:,None,None]*gen_specdopshift(data_dic['Atm']['plrange'])[:,None]*gen_specdopshift(coord_vis[pl_loc]['rv_pl'][iexp_no_plrange],v_s = coord_vis[pl_loc]['v_pl'][iexp_no_plrange]) plAtm_vis['exclu_range_input']['spec'][pl_loc] = np.zeros([len(data_dic['Atm']['CCF_mask_wav']),2,data_vis['n_in_visit']])*np.nan plAtm_vis['exclu_range_input']['spec'][pl_loc][:,:,iexp_no_plrange] = plAtm_vis['exclu_range_star']['spec'][pl_loc][:,:,iexp_no_plrange]*gen_specdopshift(coord_vis['RV_star_stelCDM'][iexp_no_plrange])*gen_specdopshift(data_dic['DI']['sysvel'][inst][vis]) return None
[docs] def update_inst(data_dic,inst,gen_dic): r"""**Update: instrument** Updates instrument-specific fields for the workflow, once all visits have been processed. Args: TBD Returns: None """ data_inst = data_dic[inst] #Common table for the reference visit # - after this stage, all modules call the common table in the star rest frame, so we update this one while keeping the original one for plotting purposes # if profiles were not aligned the common visit table below points toward the input one data_com_ref = dataload_npz(data_inst[data_inst['com_vis']]['proc_com_star_data_paths']) #Check alignment status if gen_dic['align_DI']:com_star='_star' else:com_star='' #Defining instrument table # - we define a specific path to prevent overwriting the original common grid if data_inst['comm_sp_tab'] and (gen_dic['align_DI']) and (data_dic['DI']['type'][inst]=='CCF'): data_inst['proc_com_star_data_path'] = gen_dic['save_data_dir']+'Processed_data/'+inst+'_com_up'+com_star elif gen_dic['spec_1D'] or gen_dic['CCF_from_sp']: if gen_dic['spec_1D']: data_inst['type']='spec1D' data_inst['proc_com_star_data_path'] = gen_dic['save_data_dir']+'Processed_data/'+inst+'_com_up'+com_star elif gen_dic['CCF_from_sp']: data_inst['type']='CCF' data_inst['proc_com_star_data_path'] = gen_dic['save_data_dir']+'Processed_data/CCFfromSpec/'+inst+'_com_up'+com_star if ('chrom' in data_inst['system_prop']): data_inst['system_prop']['chrom_mode'] = 'achrom' data_inst['system_prop'].pop('chrom') else: data_inst['proc_com_star_data_path'] = gen_dic['save_data_dir']+'Processed_data/'+inst+'_com_up'+com_star #Common instrument table is set to the table of the visit taken as reference datasave_npz(data_inst['proc_com_star_data_path'],data_com_ref) #Updating dimensions if gen_dic['spec_1D'] or gen_dic['CCF_from_sp']: data_inst['comm_sp_tab']=True data_inst['dim_exp'] = deepcopy(data_com_ref['dim_exp']) data_inst['nspec'] = deepcopy(data_com_ref['nspec']) data_inst['nord'] = 1 return None
[docs] def update_gen(data_dic,gen_dic): r"""**Update: generic** Updates generic fields for the workflow, once all instruments have been processed. Args: TBD Returns: None """ if gen_dic['CCF_from_sp']: if ('chrom' in data_dic['DI']['system_prop']): data_dic['DI']['system_prop']['chrom_mode'] = 'achrom' data_dic['DI']['system_prop'].pop('chrom') return None