# -*- coding: utf-8 -*- """ Created on Tue Aug 25 16:40:21 2015 @author: kervella """ import numpy as np #============================================================================== # BEGIN - USER PARAMETERS #============================================================================== # To print the output of the program in a text file, use python gravi_wavelength.py >> text_file_name.txt # Fiber coupler metrology ONLY is used if True, otherwise BOTH the FC and telescope metrology signals are used fiber_coupler_only = True boxsmooth_telmet = 150 # 150 = smoothing to the effective SC DIT of 300ms integrate_met = True # integrate the metrology of the SC telescopes over the DIT of the SC (= MET model) # Additional delay on FT for synchronization delay_ft = +0.00145 # seconds, additional delay of the center of the FT exposure with respect to the center between two time stamps ## Refraction index law in the FDDL fibers, in relative terms with respect to the laser wavelength ## n_glass(Lambda)/n_glass(Lambda_met) = 1.0 + slope_k*(Lambda/Lambda_met) 0.021 after Sept. 2 slope_k = 0.00000 # unitless slope fit_ft_delay = False nsteps, step_ft = 500, 0.000001 # +/- nsteps around delay_ft # Best fit FT time offsets: (Sept. 11 data) # > Baseline 34-S: +1.340 ms 34-P: +1.340 ms # > Baseline 32-S: +1.350 ms 32-P: +1.341 ms # > Baseline 24-S: +1.441 ms 24-P: +1.340 ms # > Baseline 31-S: +1.640 ms 31-P: +1.551 ms # > Baseline 41-S: +1.450 ms 41-P: +1.450 ms # > Baseline 21-S: +1.342 ms 21-P: +1.245 ms # #*** IMPLEMENT THE CORRECTION FOR THE SPECTRAL SHAPE OF THE SPECTRUM OF THE SOURCE ? # #-------------- ## for the data up to August 18, delay to be applied to each baseline data #slope_k = 0.0184 # unitless slope ### number of steps around values of delay_sc and step size in seconds for SC delay optimization #fit_sc_delay = False #nsteps, step_sc = 20, 0.0001 #pixmin, pixmax = 0, 233 # min and max pixel in wavelength scale for delay fitting (max 0:233) #-------------- # for the data after Sept. 3 # ["42","41","12","43","23","13"] for Sept. 3 # delay_sc = np.array([-0.056,-0.032,-0.008,+0.010,+0.040,+0.062]) # in DETECTOR order in seconds POST-2 Sept. FOR AVERAGE OF ALL 0:233 #pixmin, pixmax = 0, 233 # min and max pixel in wavelength scale for average wavelength fitting (max 0:nwave_sc) #-------------- # AD HOC SMOOTHING of resampled FT phases and MET to REDUCE THE OSCILLATIONS boxsmooth_ft = 1 # Moving average window for boxcar smoothing of the FT signals (usually 15, in unit of MET steps = 2ms) boxsmooth_met = 1 # Moving average window for boxcar smoothing of the MET signals #============================================================================== # END - USER PARAMETERS #============================================================================== import matplotlib.pyplot as plt #import os import datetime import pyfits as fits #import gravi_visual from scipy.interpolate import interp1d from scipy.ndimage.filters import gaussian_filter1d # gravi_visual_class.py import in other directory import sys sys.path.insert(0, '../../gravi_visual') import gravi_visual_class from gravi_argon import gravi_fit_argon_mr #============================================================================== # Function to clean the already unwrapped phase signals from "soft jumps" #============================================================================== def clean_unwrap_phase(phase_in,clean_step,threshold): phase = np.copy(phase_in) nframe = phase.shape[0] spring_phase = np.zeros(nframe) bump = np.zeros(nframe) time=1 while (time < (nframe-clean_step)): bump[time] = phase[time+clean_step]-phase[time] if np.abs(bump[time])>np.abs(threshold): for i in range(0,clean_step+1): # linear interpolation over the phase jump spring_phase[time+i] = spring_phase[time-1] + np.sign(bump[time])*2*np.pi*((1.0*i)/(1.0*clean_step)) time = time + clean_step + 1 # skip the bump region else: spring_phase[time]=spring_phase[time-1] time = time + 1 phase = phase - spring_phase return phase #============================================================================== # WAVELENGTH function #============================================================================== def gravi_wavelength(filename, delay_ft, boxsmooth_ft, boxsmooth_met, fiber_coupler_only, window_sc, frame_shift, binwidth): # ntel = 4 # number of telescopes nbase = 6 # number of baselines noutputs = 48 data_dir = '/Users/kervella/Pipelines/python_tools/gravi_misc/gravi_wavelength/data/' results_dir = '/Users/kervella/Pipelines/python_tools/gravi_misc/gravi_wavelength/results/' delay_sc = np.array(list(range(0,noutputs)))*0.002801 + 0.011 # delay (s) for each of the 48 individual outputs wrt FT # if window_sc == 1: # windowed readout in MR # delay_sc = np.array(range(0,noutputs))*0.002801 + 0.011 # delay (s) for each of the 48 individual outputs wrt FT # else: # non-windowed readout in MR >>> is not correct !! the windowed mode should be adapted for older epochs # delay_sc = np.array(range(0,noutputs))*0.004400 + 0.011 # delay (s) for each of the 48 individual outputs wrt FT # Create the output PDF file for the figures and the text log file import matplotlib.backends.backend_pdf pdf = matplotlib.backends.backend_pdf.PdfPages(results_dir+filename+'-FIGURES.pdf') logfile = open(results_dir+filename+'-LOG.txt','w') print("-"*80) print("*** GRAVITY wavelength calibration testbed\n") # print " Preproc file : %s " % (filename+'.fits') print(" phase_sc_ft file: %s \n" % (filename+'.fits')) logfile.write("*** GRAVITY wavelength calibration testbed\n") # logfile.write(" Preproc file : %s \n" % (filename+'.fits')) logfile.write(" phase_sc_ft file: %s \n\n" % (filename+'.fits')) #============================================================================== # Load the data files: preproc and phase_ft_sc #============================================================================== # preproc = gravi_visual_class.Preproc(data_dir+filename+'.fits') phaseftsc = gravi_visual_class.PhaseFTSC(data_dir+filename+'.fits') #============================================================================== # Reordering of the baseline outputs of SC and FT and signs #============================================================================== # List of baseline names base_list_sc=[] base_list_ft=[] for baseline in range(0,nbase): base_list_sc.append(phaseftsc.regname_sc[baseline*phaseftsc.nregion/nbase][:2]) # ["42","41","12","43","23","13"] for Sept. 3 base_list_ft.append(phaseftsc.regname_ft[baseline*phaseftsc.nregion/nbase][:2]) # ["34","32","24","31","41","21"] for Sept. 3 base_list_diff=["12", "13", "14", "23", "24", "34"] # order for the telescope differential quantities port_diff = np.array([[0,1,1],[0,2,1],[0,3,1],[1,2,1],[1,3,1],[2,3,1]]) # ports of the differential quantities, with sign port_sc = np.zeros((nbase,3),dtype=np.int) port_ft = np.zeros((nbase,3),dtype=np.int) for baseline in range(0,nbase): port_sc[baseline,:2] = phaseftsc.ports_sc[baseline*phaseftsc.nregion/nbase,:]-1 # SC ports, from 0 to 3 if port_sc[baseline,0]>port_sc[baseline,1]: port_sc[baseline,0],port_sc[baseline,1] = port_sc[baseline,1],port_sc[baseline,0] # swap port_sc[baseline,2] = -1 # mark the sign inversion else: port_sc[baseline,2] = +1 port_ft[baseline,:2] = phaseftsc.ports_ft[baseline*phaseftsc.nregion/nbase,:]-1 # FT ports, from 0 to 3 if port_ft[baseline,0]>port_ft[baseline,1]: port_ft[baseline,0],port_ft[baseline,1] = port_ft[baseline,1],port_ft[baseline,0] # swap port_ft[baseline,2] = -1 # mark the sign inversion else: port_ft[baseline,2] = +1 order_sc = np.zeros(nbase,dtype=np.int) sign_sc = np.zeros(nbase,dtype=np.int) order_ft = np.zeros(nbase,dtype=np.int) sign_ft = np.zeros(nbase,dtype=np.int) order_diff_ft = np.zeros(nbase,dtype=np.int) order_diff_sc = np.zeros(nbase,dtype=np.int) for base in range(0,nbase): for base_diff in range(0,nbase): if (port_sc[base,0:2]==port_diff[base_diff,0:2]).all(): order_sc[base] = base_diff sign_sc[base] = port_sc[base_diff,2] order_diff_sc[base_diff] = base if (port_ft[base,0:2]==port_diff[base_diff,0:2]).all(): order_ft[base] = base_diff sign_ft[base] = port_ft[base_diff,2] order_diff_ft[base_diff] = base # Signs for pre-3 sept. data if phaseftsc.header['MJD-OBS'] < 57267: # (2 Sept. 2015) sign_sc[4]=-1 sign_ft[:]=-1 # The FT signs are opposite to the SC sign_ft[4]=+1 else:# Signs for post-3 sept. data (signs of baselines 13 and 23 are inverted for FT) sign_sc = np.array([+1,+1,-1,+1,-1,-1]) sign_ft = np.array([-1,+1,-1,+1,+1,+1]) # DIFF > order_sc > SC # SC > order_diff_sc > DIFF #['12', '13', '14', '23', '24', '34'] #signs before #SC: [ 1 1 -1 1 -1 -1] #FT: [-1 -1 -1 -1 1 1] #signs after #SC: [ 1 1 -1 1 -1 -1] #FT: [-1 1 -1 1 1 1] #============================================================================== # Time scales of SC and FT #============================================================================== # METROLOGY time scale: # The metrology time scale is defined as the reference (shift = 0) # dt_met = phaseftsc.time_met[1]-phaseftsc.time_met[0] # time between each frame of the sc time_shift_met = 0.0 # No time shift for MET (reference time scale) time_met = phaseftsc.time_met + time_shift_met nframe_met = phaseftsc.time_met.shape[0] # SC time scale: dt_sc = phaseftsc.time_sc[1]-phaseftsc.time_sc[0] # time between each frame of the sc dit_sc = phaseftsc.header['HIERARCH ESO DET2 SEQ1 DIT'] # actual exposure time of SC # buffer_fddl = dt_sc - phaseftsc.header['HIERARCH ESO INS FDDL DELAY'] # delay of FDDL motion (for smoothing) # To bring back the prepoc.time_sc scale to the beginning of the cycle if phaseftsc.header['MJD-OBS'] < 57267: # (2 Sept. 2015) phaseftsc.time_sc = phaseftsc.time_sc + frame_shift*dt_sc # Time shifts to synchronize the SC in seconds time_sc = np.zeros((noutputs,phaseftsc.nframe_sc),dtype='d') # shifted time scale of SC for output in range(0,noutputs): # in DETECTOR OUTPUT order time_sc[output,:] = phaseftsc.time_sc[:] - dt_sc/2.0 + delay_sc[output] + dit_sc/2.0 # phaseftsc.time_sc is defined in the FITS files as the center of the inter-frame >>> should be changed # *** BEWARE of change of SC timescale definition on 09/02/16 *** (NOT IMPLEMENTED) # Avant le temps de la pause était au milieu de la période soit: # t_sc=start_time+row*periode+periode/2 # Maintenant, le temps est décalé de WINDOW/2 c'est donc : # t_sc=start_time+row*periode+(periode-window)/2 # C'est donc plus proche de la réalité. # # WINDOW, c'est la largeur de la fenêtre rajoutée après l’intégration, avant de passer au cycle suivant, donc le temps ou le détecteur est inactif. # Normalement DELAY+WINDOW=periode # # Pour la longueur d'onde, avant on devait prendre en compte le DELAY pour extraire dans le signal FT et MET uniquement # les frames pendant l'intégration SC, maintenant, on prend simplement les frames entre t_sc-DIT_SC/2 et t_sc+DIT_SC/2 # FT time scale: # dt_ft = phaseftsc.time_ft[1]-phaseftsc.time_ft[0] # time between each frame of the sc time_ft = np.zeros((nbase,phaseftsc.nframe_ft,phaseftsc.nwave_ft),dtype='d') # shifted time scale of FT time_shift_ft = delay_ft time_ft = phaseftsc.time_ft + time_shift_ft print("*** Parameters:") print(" Inter-frame time of SC = %.4f seconds" % dt_sc) print(" Additional time delays of each SC output (seconds):") print(delay_sc) print(" Additional delay on FT time scale = %+.3f milliseconds\n Total FT delay wrt MET time stamp = %+.3f milliseconds\n" % (1000.*delay_ft,1000.*time_shift_ft)) logfile.write("*** Parameters:\n") logfile.write(" Inter-frame time of SC = %.4f seconds\n" % dt_sc) logfile.write(" Additional time delays of each SC output (seconds):\n") for baseline in range(0,nbase): # in DETECTOR order logfile.write(delay_sc) logfile.write("\n") logfile.write(" Additional delay on FT time scale = %+.3f milliseconds\n Total FT delay wrt MET time stamp = %+.3f milliseconds\n\n" % (1000.*delay_ft, 1000.*time_shift_ft)) # print " Total delay of MET time scale = %+.3f seconds\n" % time_shift_met print(" Fiber dispersion slope = %.5f\n" % slope_k) print(" Smoothing of the FT resampled signal = %i" % boxsmooth_ft) print(" Smoothing of the MET signal = %i" % boxsmooth_met) print(" MET signal from fiber coupler only = %r\n" % fiber_coupler_only) # logfile.write(" Total delay of MET time scale = %+.3f seconds\n\n" % time_shift_met) logfile.write(" Fiber dispersion slope = %.5f\n\n" % slope_k) logfile.write(" Smoothing of the FT resampled signal = %i\n" % boxsmooth_ft) logfile.write(" Smoothing of the MET signal = %i\n" % boxsmooth_met) logfile.write(" MET signal from fiber coupler only = %r\n\n" % fiber_coupler_only) if 'HIERARCH ESO INS MLC WAVELENG' in phaseftsc.header: Lambda_met=phaseftsc.header['HIERARCH ESO INS MLC WAVELENG']/1000. # Metrology wavelength in microns (1.908 microns) else: Lambda_met = 1.908287 #============================================================================== # Average voltages of the FDDLs (FDDL table) over the sequence #============================================================================== print("*** Average voltages of FDDLs:") print(" FT FDDLs (V): %.4f %.4f %.4f %.4f (nominal: 1.7716 2.1313 2.2658 2.0524)" % (np.mean(phaseftsc.fddl_ft_pos[:,0]),np.mean(phaseftsc.fddl_ft_pos[:,1]), np.mean(phaseftsc.fddl_ft_pos[:,2]),np.mean(phaseftsc.fddl_ft_pos[:,3]))) print(" SC FDDLs (V): %.4f %.4f %.4f %.4f (nominal: 1.8323 2.4473 2.6608 0.6833)\n" % (np.mean(phaseftsc.fddl_sc_pos[:,0]),np.mean(phaseftsc.fddl_sc_pos[:,1]), np.mean(phaseftsc.fddl_sc_pos[:,2]),np.mean(phaseftsc.fddl_sc_pos[:,3]))) logfile.write(" Average voltages of FDDLs:\n") logfile.write(" FT FDDLs (V): %.4f %.4f %.4f %.4f (nominal: 1.7716 2.1313 2.2658 2.0524)\n" % (np.mean(phaseftsc.fddl_ft_pos[:,0]),np.mean(phaseftsc.fddl_ft_pos[:,1]), np.mean(phaseftsc.fddl_ft_pos[:,2]),np.mean(phaseftsc.fddl_ft_pos[:,3]))) logfile.write(" SC FDDLs (V): %.4f %.4f %.4f %.4f (nominal: 1.8323 2.4473 2.6608 0.6833)\n\n" % (np.mean(phaseftsc.fddl_sc_pos[:,0]),np.mean(phaseftsc.fddl_sc_pos[:,1]), np.mean(phaseftsc.fddl_sc_pos[:,2]),np.mean(phaseftsc.fddl_sc_pos[:,3]))) # #============================================================================== # # Fiber coupler unwrapped metrology OPL per telescope (in microns) # #============================================================================== # # unwrapped_ft = np.zeros((ntel,nframe_met),dtype='d') # unwrapped_sc = np.zeros((ntel,nframe_met),dtype='d') # unwrapped_met_ft = np.zeros((nbase,nframe_met),dtype='d') # unwrapped_met_sc = np.zeros((nbase,nframe_met),dtype='d') # # if fiber_coupler_only == True: # Use only the fiber coupler diodes # print "*** Unwrapping of the metrology OPL per telescope using ONLY the fiber coupler signals" # for tel in range(0,ntel): # phase = np.arctan2(preproc.metrology_volt[:,2*tel+64], preproc.metrology_volt[:,2*tel+65]) # unwrapped_ft[tel] = np.unwrap(phase)*Lambda_met/(2*np.pi) ## unwrapped_ft = phaseftsc.met_opl # for tel in range(0,ntel): # phase = np.arctan2(preproc.metrology_volt[:,2*tel+72], preproc.metrology_volt[:,2*tel+73]) # unwrapped_sc[tel] = np.unwrap(phase)*Lambda_met/(2*np.pi) # else: # # metrology vectors for the 4 telescopes of each combiner FROM THE FIBER COUPLER AND TELESCOPES # print "*** Unwrapping of the metrology OPL per telescope using the telescope diodes AND fiber coupler signals" # ndiode = 4 # number of diodes per telescope + 1 for the fiber coupler diode # phase_tel = np.zeros((ndiode,nframe_met),dtype='d') # phase_fc = np.zeros(nframe_met,dtype='d') # # for tel in range(0,ntel): # # for diode in range(0,ndiode): # 4 telescope diodes for FT, unwrap each diode signal individually before averaging # phase_tel[diode,:] = np.unwrap(np.arctan2(preproc.metrology_volt[:,(8*tel)+2*diode], preproc.metrology_volt[:,(8*tel+1)+2*diode])) # # mean phase for the telescope diodes: # mean_phase_tel = np.mean(phase_tel[0:(ndiode-1),:],axis=0) # # phase_fc = np.unwrap(np.arctan2(preproc.metrology_volt[:,(2*tel)+64], preproc.metrology_volt[:,(2*tel)+65])) # # residual phase telescope - fiber coupler: # diff_phase = mean_phase_tel - phase_fc # # clean_unwrap=4 # threshold = 1.5*np.pi # diff_phase = clean_unwrap_phase(diff_phase,clean_unwrap,threshold) # # #plt.plot(diff_phase/(2*np.pi),label="FT"+str(tel+1)) # # Smoothing of the diff_phase signal: # diff_phase = np.convolve(diff_phase, np.ones((boxsmooth_telmet,))/boxsmooth_telmet, mode='same') # #plt.plot(diff_phase/(2*np.pi),label="FT"+str(tel+1)) # unwrapped_ft[tel] = (phase_fc + diff_phase)*Lambda_met/(2*np.pi) # # for diode in range(0,ndiode): # SC # phase_tel[diode,:] = np.unwrap(np.arctan2(preproc.metrology_volt[:,(8*tel+32)+2*diode], preproc.metrology_volt[:,(8*tel+33)+2*diode])) # mean_phase_tel = np.mean(phase_tel[0:(ndiode-1),:],axis=0) # # phase_fc = np.unwrap(np.arctan2(preproc.metrology_volt[:,(2*tel)+72], preproc.metrology_volt[:,(2*tel)+73])) # diff_phase = (mean_phase_tel - phase_fc) # # plt.plot(diff_phase/(2*np.pi),label="SC"+str(tel+1)) # diff_phase = clean_unwrap_phase(diff_phase,clean_unwrap,threshold) # # plt.plot(diff_phase/(2*np.pi),label="SC"+str(tel+1)) # # Smoothing of the diff_phase signal: # diff_phase = np.convolve(diff_phase, np.ones((boxsmooth_telmet,))/boxsmooth_telmet, mode='same') # unwrapped_sc[tel] = (phase_fc + diff_phase)*Lambda_met/(2*np.pi) # # for baseline in range(0,nbase): # with the sign and baseline numbering convention of the DIFFERENTIAL signals # unwrapped_met_ft[order_ft[baseline],:] = (unwrapped_ft[port_ft[baseline,1]]-unwrapped_ft[port_ft[baseline,0]]) # unwrapped_met_sc[order_sc[baseline],:] = (unwrapped_sc[port_sc[baseline,1]]-unwrapped_sc[port_sc[baseline,0]]) # # # #============================================================================== # # Integration/smoothing of the SC metrology using the SC integration windows # #============================================================================== # if integrate_met == True: # print "*** Resampling/smoothing of the SC metrology using the SC integration windows" # # box_sc_met = int(dit_sc/(2.*dt_met)) # Box half-size for integration of the SC MET signals # resamp_sc_met = np.zeros((nbase,phaseftsc.nframe_sc),dtype='d') # # Final smoothing to match the actual shape of the SC steps (in number of MET steps, 50 steps = 100 ms) # smooth_sc_steps = int(np.abs(buffer_fddl)/(2.0*dt_met)) # # # meanwave = int(phaseftsc.nwave_sc/2.0) # for baseline in range(0,nbase): # in DIFF order # for frame_sc in range(0,phaseftsc.nframe_sc): # index_sc_met = np.min(np.where(time_met >= (phaseftsc.time_sc[frame_sc]))) # # corresponding index in the metrology sequence # # resamp_sc_met[baseline,frame_sc] = np.mean(unwrapped_met_sc[baseline,(index_sc_met-box_sc_met):(index_sc_met+box_sc_met)]) # resamp_sc_met[baseline,frame_sc] = np.mean(unwrapped_met_sc[baseline,(index_sc_met-box_sc_met):(index_sc_met+box_sc_met)]) # # # re-interpolation at the frequency of the metrology # near_sc_met = interp1d((phaseftsc.time_sc[:]), # resamp_sc_met[baseline,:],kind='nearest',bounds_error=False, fill_value=0.0) # # Gaussian smoothing # near_sc_met_int = near_sc_met(time_met) # #unwrapped_met_sc[baseline,:] = gaussian_filter1d(near_sc_met_int,smooth_sc_steps) # unwrapped_met_sc[baseline,:]=near_sc_met_int # # # indicator of exposure start on metrology signal ## indexmet_sc = np.zeros((phaseftsc.nframe_sc),dtype=np.int) ## for time in range(0,phaseftsc.nframe_sc): ## indexmet_sc[time] = int((phaseftsc.time_sc[time]-dit_sc/2.)//dt_met) ## for baseline in range(0,nbase): # with the sign and baseline numbering convention of the DIFFERENTIAL signals ## unwrapped_met_sc[baseline,indexmet_sc] += 0.5 ## for time in range(0,phaseftsc.nframe_sc): ## indexmet_sc[time] = int((phaseftsc.time_sc[time]+dit_sc/2.)//dt_met) ## for baseline in range(0,nbase): # with the sign and baseline numbering convention of the DIFFERENTIAL signals ## unwrapped_met_sc[baseline,indexmet_sc] += 0.3 #============================================================================== # Computation of the differential metrology signal per baseline #============================================================================== print("*** Computation of the differential metrology signal per baseline") # for baseline in range(0,nbase): # diff_met[baseline,:] = unwrapped_met_sc[baseline,:] - unwrapped_met_ft[baseline,:] diff_met = np.zeros((nbase,nframe_met),dtype='d') for baseline in range(0,nbase): diff_met[baseline,:] = (phaseftsc.met_dopl[port_diff[baseline,0],:]-phaseftsc.met_dopl[port_diff[baseline,1],:])*port_diff[baseline,2] plt.close('all') #============================================================================== # Binning of the HR spectra by a factor binwidth before computing the phases #============================================================================== if ('HIGH' in phaseftsc.header['HIERARCH ESO INS SPEC RES']) & (binwidth > 1): # HR case, binning by binwidth print("*** Binning of the HR spectra by a factor %i for phase unwrapping" % binwidth) nwave_rebin = phaseftsc.nwave_sc//binwidth # integer division spect_sc = np.zeros((phaseftsc.nregion,phaseftsc.nframe_sc,nwave_rebin),dtype='d') for region in range(0,phaseftsc.nregion): for timestep in range(0,phaseftsc.nframe_sc): spectrum = phaseftsc.spectrum_sc[region,timestep,:] spect_sc[region,timestep,:] = np.mean(spectrum[:(spectrum.size // binwidth) * binwidth].reshape(-1, binwidth),axis=1) else: # MR case, no binning nwave_rebin = phaseftsc.nwave_sc spect_sc = np.copy(phaseftsc.spectrum_sc) #============================================================================== # Time shift of the SC individual spectrum signals #============================================================================== # correction of the offset and interpolation of each output to bring them to the common phaseftsc.time_sc scale # with this shift, we have a common SC time scale with the FT and the MET. The SC DIT is centered on phaseftsc.time_sc ticks. for output in range(0,noutputs): for wave in range(0,nwave_rebin): spect_sc[output,:,wave] = np.interp(phaseftsc.time_sc, time_sc[output,:], spect_sc[output,:,wave]) #============================================================================== # Determination of the phase of the SC as a function of time using ellipse fitting #============================================================================== print("*** Computation of the SC phases using ellipse fitting") # self.spectrum_sc = output, time, wavelength unwrapped_phi_sc_s = np.zeros((nbase,phaseftsc.nframe_sc,nwave_rebin),dtype='d') unwrapped_phi_sc_p = np.zeros((nbase,phaseftsc.nframe_sc,nwave_rebin),dtype='d') import scipy.optimize as optimization def ell(x, a, b, c, d, e): return np.sqrt(np.square(a*x[0,:] + b*x[1,:] +c) + np.square(d*x[1,:] + e)) phase_s = np.zeros((nbase,phaseftsc.nframe_sc,phaseftsc.nwave_sc),dtype='d') phase_p = np.zeros((nbase,phaseftsc.nframe_sc,phaseftsc.nwave_sc),dtype='d') x0 = [0.5]*5 # initial guess for baseline in range(0,nbase): for wave in range(0,nwave_rebin): X = spect_sc[0+8*baseline,:,wave]-spect_sc[2+8*baseline,:,wave] Y = spect_sc[4+8*baseline,:,wave]-spect_sc[6+8*baseline,:,wave] # normalization by average radius to avoid instability meanrad = np.sqrt(np.square(np.mean(X))+np.square(np.mean(Y))) X /= meanrad Y /= meanrad ellipse,ellipsecov = optimization.curve_fit(ell, [X,Y], [1.0]*phaseftsc.nframe_sc, x0) X2 = ellipse[0]*X + ellipse[1]*Y + ellipse[2] Y2 = ellipse[3]*Y + ellipse[4] phase_s[baseline,:,wave] = np.arctan2(X2,Y2) X = spect_sc[1+8*baseline,:,wave]-spect_sc[3+8*baseline,:,wave] Y = spect_sc[5+8*baseline,:,wave]-spect_sc[7+8*baseline,:,wave] meanrad = np.sqrt(np.square(np.mean(X))+np.square(np.mean(Y))) X /= meanrad Y /= meanrad ellipse,ellipsecov = optimization.curve_fit(ell, [X,Y], [1.0]*phaseftsc.nframe_sc, x0) X2 = ellipse[0]*X + ellipse[1]*Y + ellipse[2] Y2 = ellipse[3]*Y + ellipse[4] phase_p[baseline,:,wave] = np.arctan2(X2,Y2) for baseline in range(0,nbase): for wave in range(0,nwave_rebin): unwrapped_phi_sc_s[order_sc[baseline],:,wave] = np.unwrap(phase_s[baseline,:,wave])*sign_sc[baseline] # unwrap over time WITH SIGN AND DIFF order unwrapped_phi_sc_p[order_sc[baseline],:,wave] = np.unwrap(phase_p[baseline,:,wave])*sign_sc[baseline] # unwrap over time WITH SIGN AND DIFF order unwrapped_phi_sc_s[baseline,:,:] = np.unwrap(unwrapped_phi_sc_s[baseline,:,:],axis=1) # unwrapping over wavelength (to remove remaining jumps) unwrapped_phi_sc_p[baseline,:,:] = np.unwrap(unwrapped_phi_sc_p[baseline,:,:],axis=1) #============================================================================== # Resampling of the SC signal at the frequency of the MET (that is faster) #============================================================================== print("*** Resampling of the SC signal at the MET frequency using nearest neighbor") phi_sc_s_int = np.zeros((nbase,phaseftsc.nframe_met,nwave_rebin),dtype='d') phi_sc_p_int = np.zeros((nbase,phaseftsc.nframe_met,nwave_rebin),dtype='d') # Smoothing to match the actual shape of the SC steps (in number of MET steps, 50 steps = 100 ms) # smooth_sc_steps = int(np.abs(dt_sc-delay_fddl)/(2.0*dt_met)) smooth_sc_steps = 1 for baseline in range(0,nbase): # in DIFF order for wave in range(0,nwave_rebin): # Nearest neighbor interpolation near_s = interp1d(phaseftsc.time_sc, unwrapped_phi_sc_s[baseline,:,wave],kind='nearest',bounds_error=False, fill_value=0.0) s_int = np.unwrap(near_s(time_met)) near_p = interp1d(phaseftsc.time_sc, unwrapped_phi_sc_p[baseline,:,wave],kind='nearest',bounds_error=False, fill_value=0.0) p_int = np.unwrap(near_p(time_met)) # Gaussian smoothing phi_sc_s_int[baseline,:,wave]=gaussian_filter1d(s_int,smooth_sc_steps) phi_sc_p_int[baseline,:,wave]=gaussian_filter1d(p_int,smooth_sc_steps) phi_sc_s_int[baseline,:,:]=np.unwrap(phi_sc_s_int[baseline,:,:],axis=1) # unwrap in wavelength to avoid any residual jump phi_sc_p_int[baseline,:,:]=np.unwrap(phi_sc_p_int[baseline,:,:],axis=1) #============================================================================== # Determination of the phase of the FT as a function of time using ellipse fitting #============================================================================== print("*** Computation of the FT phases using ellipse fitting") # self.spectrum_sc = output, time, wavelength unwrapped_phi_ft_s = np.zeros((nbase,phaseftsc.nframe_ft,phaseftsc.nwave_ft),dtype='d') unwrapped_phi_ft_p = np.zeros((nbase,phaseftsc.nframe_ft,phaseftsc.nwave_ft),dtype='d') x0 = [0.5]*5 # initial guess phase_s = np.zeros((nbase,phaseftsc.nframe_ft,phaseftsc.nwave_ft),dtype='d') phase_p = np.zeros((nbase,phaseftsc.nframe_ft,phaseftsc.nwave_ft),dtype='d') for baseline in range(0,nbase): for wave in range(0,phaseftsc.nwave_ft): # X = phaseftsc.spectrum_ft[0+8*baseline,:,wave]-phaseftsc.spectrum_ft[2+8*baseline,:,wave] # Y = phaseftsc.spectrum_ft[4+8*baseline,:,wave]-phaseftsc.spectrum_ft[6+8*baseline,:,wave] X = phaseftsc.spectrum_ft[0+8*baseline,:,wave]-phaseftsc.spectrum_ft[2+8*baseline,:,wave] Y = phaseftsc.spectrum_ft[4+8*baseline,:,wave]-phaseftsc.spectrum_ft[6+8*baseline,:,wave] meanrad = np.sqrt(np.square(np.mean(X))+np.square(np.mean(Y))) X /= meanrad Y /= meanrad ellipse,ellipsecov = optimization.curve_fit(ell, [X,Y], [1.0]*phaseftsc.nframe_ft, x0) X2 = ellipse[0]*X + ellipse[1]*Y + ellipse[2] Y2 = ellipse[3]*Y + ellipse[4] phase_s[baseline,:,wave] = np.arctan2(X2,Y2) # X = phaseftsc.spectrum_ft[1+8*baseline,:,wave]-phaseftsc.spectrum_ft[3+8*baseline,:,wave] # Y = phaseftsc.spectrum_ft[5+8*baseline,:,wave]-phaseftsc.spectrum_ft[7+8*baseline,:,wave] X = phaseftsc.spectrum_ft[1+8*baseline,:,wave]-phaseftsc.spectrum_ft[3+8*baseline,:,wave] Y = phaseftsc.spectrum_ft[5+8*baseline,:,wave]-phaseftsc.spectrum_ft[7+8*baseline,:,wave] meanrad = np.sqrt(np.square(np.mean(X))+np.square(np.mean(Y))) X /= meanrad Y /= meanrad ellipse,ellipsecov = optimization.curve_fit(ell, [X,Y], [1.0]*phaseftsc.nframe_ft, x0) X2 = ellipse[0]*X + ellipse[1]*Y + ellipse[2] Y2 = ellipse[3]*Y + ellipse[4] phase_p[baseline,:,wave] = np.arctan2(X2,Y2) for baseline in range(0,nbase): for wave in range(0,phaseftsc.nwave_ft): # DETECTOR order unwrapped_phi_ft_s[order_ft[baseline],:,wave] = np.unwrap(phase_s[baseline,:,wave])*sign_ft[baseline] # unwrap over time WITH SIGN CORRECTION and DIFF order unwrapped_phi_ft_p[order_ft[baseline],:,wave] = np.unwrap(phase_p[baseline,:,wave])*sign_ft[baseline] # unwrap over time WITH SIGN CORRECTION and DIFF order unwrapped_phi_ft_s[baseline,:,:] = np.unwrap(unwrapped_phi_ft_s[baseline,:,:],axis=1) # unwrapping over wavelength (to remove remaining jumps) unwrapped_phi_ft_p[baseline,:,:] = np.unwrap(unwrapped_phi_ft_p[baseline,:,:],axis=1) #============================================================================== # Resampling of the FT signal at the frequency of the MET (that is slower) #============================================================================== print("*** Resampling of the FT signal at the MET frequency with %i pixel AD HOC smoothing" % boxsmooth_ft) phi_ft_s_int = np.zeros((nbase,phaseftsc.nframe_met,phaseftsc.nwave_ft),dtype='d') phi_ft_p_int = np.zeros((nbase,phaseftsc.nframe_met,phaseftsc.nwave_ft),dtype='d') for baseline in range(0,nbase): for wave in range(0,phaseftsc.nwave_ft): phi_ft_s_int[baseline,:,wave] = np.interp(time_met, time_ft, unwrapped_phi_ft_s[baseline,:,wave]) phi_ft_p_int[baseline,:,wave] = np.interp(time_met, time_ft, unwrapped_phi_ft_p[baseline,:,wave]) phi_ft_s_int[baseline,:,wave] = np.convolve(phi_ft_s_int[baseline,:,wave], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') phi_ft_p_int[baseline,:,wave] = np.convolve(phi_ft_p_int[baseline,:,wave], np.ones((boxsmooth_ft,))/boxsmooth_ft, mode='same') phi_ft_s_int[baseline,:,:]=np.unwrap(phi_ft_s_int[baseline,:,:],axis=1) # unwrap in wavelength to avoid any residual jump phi_ft_p_int[baseline,:,:]=np.unwrap(phi_ft_p_int[baseline,:,:],axis=1) #============================================================================== # Smoothing of the differential MET signals #============================================================================== print("*** Smoothing of the MET signal with %i pixel AD HOC window\n" % boxsmooth_met) for baseline in range(0,nbase): diff_met[baseline,:]=np.convolve(diff_met[baseline,:], np.ones((boxsmooth_met,))/boxsmooth_met, mode='same') #============================================================================== # Average wavelength determination for FT and SC #============================================================================== print("*** Computation of the average wavelength of the SC and FT per baseline") avg_wave_sc_s = np.zeros((nbase,2)) # wavelength, error bar avg_wave_sc_p = np.zeros((nbase,2)) avg_wave_ft_s = np.zeros((nbase,2)) avg_wave_ft_p = np.zeros((nbase,2)) opd_shift_s = np.zeros((nbase,2)) opd_shift_p = np.zeros((nbase,2)) phi_sc_avg_s= np.zeros((nbase,phaseftsc.nframe_met)) phi_sc_avg_p= np.zeros((nbase,phaseftsc.nframe_met)) phi_ft_avg_s= np.zeros((nbase,phaseftsc.nframe_met)) phi_ft_avg_p= np.zeros((nbase,phaseftsc.nframe_met)) import scipy.optimize as optimization x0 = [2.2,2.2,0.0] # initial guess # minfitwave = 0 # minimum wavelength channel for the computation of the average opd over wavelengths # maxfitwave = phaseftsc.nwave_sc # maximum wavelength channel minfitwavesc = 0 # minimum wavelength channel for the computation of the average opd over wavelengths maxfitwavesc = phaseftsc.nwave_sc # maximum wavelength channel minfitwaveft = 0 # minimum wavelength channel for the computation of the average opd over wavelengths maxfitwaveft = phaseftsc.nwave_ft # maximum wavelength channel for baseline in range(0,nbase): phi_sc_avg_s[baseline,:] = np.mean(phi_sc_s_int[baseline,:,minfitwavesc:maxfitwavesc],axis=1) phi_sc_avg_p[baseline,:] = np.mean(phi_sc_p_int[baseline,:,minfitwavesc:maxfitwavesc],axis=1) phi_ft_avg_s[baseline,:] = np.mean(phi_ft_s_int[baseline,:,minfitwaveft:maxfitwaveft],axis=1) phi_ft_avg_p[baseline,:] = np.mean(phi_ft_p_int[baseline,:,minfitwaveft:maxfitwaveft],axis=1) def lin2(x, a, b, c): return (a+x[2])*x[0] + (b+x[2])*x[1] + c def fitlambda_mean(phi_sc,phi_ft,diff_met,lambda_met): ydata = diff_met wave_sc = np.zeros(2,dtype='d') wave_ft = np.zeros(2,dtype='d') opd_shift = np.zeros(2,dtype='d') x0 = [2.2-lambda_met,2.2-lambda_met,-2.0] # initial guess xdata = [phi_sc[:]/(2.0*np.pi),-1.0*phi_ft[:]/(2.0*np.pi),lambda_met] # NOTE THE MINUS SIGN to fit SC-FT opt2 = optimization.curve_fit(lin2, xdata, ydata, x0) wave_sc = [opt2[0][0]+lambda_met, np.sqrt(opt2[1][0,0])] wave_ft = [opt2[0][1]+lambda_met, np.sqrt(opt2[1][1,1])] opd_shift =[opt2[0][2], np.sqrt(opt2[1][2,2])] return wave_sc,wave_ft,opd_shift # ******************************************************************************* # **** Do not consider the first 6 seconds and the last 6 seconds in the fit **** # ******************************************************************************* # They are affected by spurious effects, maybe coming from the interpolations. mintime = np.min(np.where(time_met>6.2)) maxtime = phaseftsc.nframe_met - mintime # FOR TESTS: # mintime = np.min(np.where(time_met>30)) # maxtime = phaseftsc.nframe_met-np.min(np.where(time_met>6.2)) # Linear fit of [SC-FT] vs [METSC - METFT] for baseline in range(0,nbase): # DIFF order avg_wave_sc_s[baseline,:],avg_wave_ft_s[baseline,:],opd_shift_s[baseline,:] = \ fitlambda_mean(phi_sc_avg_s[baseline,mintime:maxtime], phi_ft_avg_s[baseline,mintime:maxtime], diff_met[baseline,mintime:maxtime],Lambda_met) avg_wave_sc_p[baseline,:],avg_wave_ft_p[baseline,:],opd_shift_p[baseline,:] = \ fitlambda_mean(phi_sc_avg_p[baseline,mintime:maxtime], phi_ft_avg_p[baseline,mintime:maxtime], diff_met[baseline,mintime:maxtime],Lambda_met) # Force the wavelength of each baseline to the average of all baselines # ********************************************************************* # This gives a better stability of the derived wavelength scales grand_mean_wave_sc = np.mean([avg_wave_sc_s[:,0],avg_wave_sc_p[:,0]]) avg_wave_sc_s[:,0] = grand_mean_wave_sc avg_wave_sc_p[:,0] = grand_mean_wave_sc # grand_mean_wave_ft = np.mean([avg_wave_ft_s[:,0],avg_wave_ft_p[:,0]]) # avg_wave_ft_s[:,0] = grand_mean_wave_ft # avg_wave_ft_p[:,0] = grand_mean_wave_ft #============================================================================== # Display and plot of the average wavelengths #============================================================================== print(" Average wavelengths of SC baselines: ") logfile.write(" Average wavelengths of SC baselines: \n") for baseline in range(0,nbase): # order of the SC baselines on the detector print(" > %s-S = %.5f +/- %.5f mic %s-P = %.5f +/- %.5f mic" % \ (base_list_sc[baseline],avg_wave_sc_s[order_sc[baseline],0],avg_wave_sc_s[order_sc[baseline],1], base_list_sc[baseline],avg_wave_sc_p[order_sc[baseline],0],avg_wave_sc_p[order_sc[baseline],1])) logfile.write(" > %s-S = %.5f +/- %.5f mic %s-P = %.5f +/- %.5f mic\n" % \ (base_list_sc[baseline],avg_wave_sc_s[order_sc[baseline],0],avg_wave_sc_s[order_sc[baseline],1], base_list_sc[baseline],avg_wave_sc_p[order_sc[baseline],0],avg_wave_sc_p[order_sc[baseline],1])) print(" Average of all SC baselines: S: %.5f mic P: %.5f mic" % (np.mean(avg_wave_sc_s[:,0]),np.mean(avg_wave_sc_p[:,0]))) logfile.write(" Average of all SC baselines: S: %.5f mic P: %.5f mic\n" % (np.mean(avg_wave_sc_s[:,0]),np.mean(avg_wave_sc_p[:,0]))) print(" Average wavelengths of FT baselines:") logfile.write(" Average wavelengths of FT baselines:\n") for baseline in range(0,nbase): # order of the FT baselines on the detector print(" > %s-S = %.5f +/- %.5f mic %s-P = %.5f +/- %.5f mic" % \ (base_list_ft[baseline],avg_wave_ft_s[order_ft[baseline],0],avg_wave_ft_s[order_ft[baseline],1], base_list_ft[baseline],avg_wave_ft_p[order_ft[baseline],0],avg_wave_ft_p[order_ft[baseline],1])) logfile.write(" > %s-S = %.5f +/- %.5f mic %s-P = %.5f +/- %.5f mic\n" % \ (base_list_ft[baseline],avg_wave_ft_s[order_ft[baseline],0],avg_wave_ft_s[order_ft[baseline],1], base_list_ft[baseline],avg_wave_ft_p[order_ft[baseline],0],avg_wave_ft_p[order_ft[baseline],1])) print(" Average of all FT baselines: S: %.5f mic P: %.5f mic" % (np.mean(avg_wave_ft_s[:,0]),np.mean(avg_wave_ft_p[:,0]))) logfile.write(" Average of all FT baselines: S: %.5f mic P: %.5f mic\n" % (np.mean(avg_wave_ft_s[:,0]),np.mean(avg_wave_ft_p[:,0]))) print(" dOPD fit constant term:") logfile.write(" dOPD fit constant term:\n") for baseline in range(0,nbase): # DIFF order print(" > %s-S = %+.5f +/- %.5f mic %s-P = %+.5f +/- %.5f mic" % \ (base_list_diff[baseline],opd_shift_s[baseline,:][0],opd_shift_s[baseline,:][1], \ base_list_diff[baseline],opd_shift_p[baseline,:][0],opd_shift_p[baseline,:][1])) logfile.write(" > %s-S = %+.5f +/- %.5f mic %s-P = %+.5f +/- %.5f mic\n" % \ (base_list_diff[baseline],opd_shift_s[baseline,:][0],opd_shift_s[baseline,:][1], \ base_list_diff[baseline],opd_shift_p[baseline,:][0],opd_shift_p[baseline,:][1])) plt.figure(figsize=(20,10)) x=list(range(0,nbase)) plt.subplot(121) plt.plot(x,avg_wave_sc_s[order_sc[:],0],label="SC S") plt.plot(x,avg_wave_sc_p[order_sc[:],0],label="SC P") plt.xticks(x,base_list_sc[:]) plt.ylabel('Mean wavelength (mic)') plt.xlabel('Baseline') plt.margins(0.2) plt.legend() plt.subplot(122) plt.plot(x,avg_wave_ft_s[order_ft[:],0],label="FT S") plt.plot(x,avg_wave_ft_p[order_ft[:],0],label="FT P") plt.xticks(x,base_list_ft[:]) plt.ylabel('Mean wavelength (mic)') plt.xlabel('Baseline') plt.margins(0.2) plt.legend() pdf.savefig() #============================================================================== # Precise determination of the FT time delays #============================================================================== if fit_ft_delay: # set to True to compute the FT time delays print("*** Determination of the time delay for FT through residual minimization with dMET") logfile.write("\n*** Determination of the time delay for FT through residual minimization with dMET\n") avg_unwrapped_phi_ft_s = np.zeros((nbase,phaseftsc.nframe_ft)) avg_unwrapped_phi_ft_p = np.zeros((nbase,phaseftsc.nframe_ft)) avg_phi_ft_s_int = np.zeros((nbase,phaseftsc.nframe_met)) avg_phi_ft_p_int = np.zeros((nbase,phaseftsc.nframe_met)) # Compute average of FT unwrapped phase over wavelengths avg_unwrapped_phi_ft_s[:,:] = np.mean(unwrapped_phi_ft_s[:,:,:],axis=2) avg_unwrapped_phi_ft_p[:,:] = np.mean(unwrapped_phi_ft_p[:,:,:],axis=2) std_residual_dOPD_s = np.zeros((nbase,2*nsteps+1)) std_residual_dOPD_p = np.zeros((nbase,2*nsteps+1)) delay_ft_s_fit = np.zeros(nbase,dtype='d') delay_ft_p_fit = np.zeros(nbase,dtype='d') for baseline in range(0,nbase): # in DIFF order for delta_ft in range(-nsteps,nsteps+1): near_s = interp1d(time_ft[:] + delta_ft*step_ft, avg_unwrapped_phi_ft_s[baseline,:], kind='nearest',bounds_error=False, fill_value=0.0) s_int = np.unwrap(near_s(time_met)) near_p = interp1d(time_ft[:] + delta_ft*step_ft, avg_unwrapped_phi_ft_p[baseline,:], kind='nearest',bounds_error=False, fill_value=0.0) p_int = np.unwrap(near_p(time_met)) # Gaussian smoothing # avg_phi_ft_s_int[baseline,:]=gaussian_filter1d(s_int,smooth_sc_steps) # avg_phi_ft_p_int[baseline,:]=gaussian_filter1d(p_int,smooth_sc_steps) avg_phi_ft_s_int[baseline,:]=s_int avg_phi_ft_p_int[baseline,:]=p_int std_residual_dOPD_s[baseline,delta_ft+nsteps] = \ np.std((phi_sc_avg_s[baseline,mintime:maxtime]*avg_wave_sc_s[baseline,0] - \ avg_phi_ft_s_int[baseline,mintime:maxtime]*avg_wave_ft_s[baseline,0])/(2*np.pi) - \ diff_met[baseline,mintime:maxtime]) std_residual_dOPD_p[baseline,delta_ft+nsteps] = \ np.std((phi_sc_avg_p[baseline,mintime:maxtime]*avg_wave_sc_p[baseline,0] - \ avg_phi_ft_p_int[baseline,mintime:maxtime]*avg_wave_ft_p[baseline,0])/(2*np.pi) - \ diff_met[baseline,mintime:maxtime]) delay_ft_s_fit[baseline] = (np.argmin(std_residual_dOPD_s[baseline,:])-nsteps)*step_ft + time_shift_ft delay_ft_p_fit[baseline] = (np.argmin(std_residual_dOPD_p[baseline,:])-nsteps)*step_ft + time_shift_ft print(" Best fit FT time offsets:") for baseline in range(0,nbase): # in DETECTOR order print(" > Baseline %s-S: %+.3f ms %s-P: %+.3f ms" % \ (base_list_ft[baseline], delay_ft_s_fit[order_ft[baseline]]*1000., base_list_ft[baseline], delay_ft_p_fit[order_ft[baseline]]*1000.)) logfile.write(" Best fit FT time offsets:\n") for baseline in range(0,nbase): logfile.write(" > Baseline %s-S: %+.3f ms %s-P: %+.3f ms\n" % \ (base_list_ft[baseline], delay_ft_s_fit[order_ft[baseline]]*1000., base_list_ft[baseline], delay_ft_p_fit[order_ft[baseline]]*1000.)) #============================================================================== # Verification of the fit quality #============================================================================== residual_s = np.zeros(nbase,dtype='d') residual_p = np.zeros(nbase,dtype='d') print(" RMS of residuals dOPD-dMET baseline:") logfile.write(" RMS of residuals dOPD-dMET baseline:\n") for baseline in range(0,nbase): # S & P polarizations #[0]: # for only one baseline # DIFF baseline order residual_s[baseline] = np.std((np.mean(phi_sc_s_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_sc_s[baseline,0] - np.mean(phi_ft_s_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_ft_s[baseline,0])/(2*np.pi) + opd_shift_s[baseline,0]-diff_met[baseline,mintime:maxtime]) residual_p[baseline] = np.std((np.mean(phi_sc_p_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_sc_p[baseline,0] - np.mean(phi_ft_p_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_ft_p[baseline,0])/(2*np.pi) + opd_shift_p[baseline,0]-diff_met[baseline,mintime:maxtime]) for baseline in range(0,nbase): # SC baseline order print(" > %s-S = %.2f nm %s-P = %.2f nm" % \ (base_list_sc[baseline], 1000.*residual_s[order_sc[baseline]], base_list_sc[baseline], 1000.*residual_p[order_sc[baseline]])) logfile.write(" > %s-S = %.2f nm %s-P = %.2f nm\n" % \ (base_list_sc[baseline], 1000.*residual_s[order_sc[baseline]], base_list_sc[baseline], 1000.*residual_p[order_sc[baseline]])) print(" Global average dOPD-dMET RMS residual: S = %.2f nm P = %.2f nm\n" % (1000.*np.mean(residual_s),1000.*np.mean(residual_p))) logfile.write(" Global average dOPD-dMET RMS residual: S = %.2f nm P = %.2f nm\n\n" % (1000.*np.mean(residual_s),1000.*np.mean(residual_p))) if True: for baseline in range(0,nbase): # S & P polarizations # replace by [0] for only one baseline in DIFF order plt.figure(figsize=(20, 10)) ax1=plt.subplot(221) ax1.plot(time_met[mintime:maxtime], (np.mean(phi_sc_s_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_sc_s[baseline,0] - np.mean(phi_ft_s_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_ft_s[baseline,0])/(2*np.pi) + opd_shift_s[baseline,0],label='SC-FT '+base_list_diff[baseline]+'-S') ax1.plot(time_met[mintime:maxtime], diff_met[baseline,mintime:maxtime], label='MET_SC-MET_FT') plt.ylabel('dOPD base %s-S (mic)' % base_list_diff[baseline]) plt.legend() ax2=plt.subplot(222, sharex=ax1, sharey=ax1) ax2.plot(time_met[mintime:maxtime], (np.mean(phi_sc_p_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_sc_p[baseline,0] - np.mean(phi_ft_p_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_ft_p[baseline,0])/(2*np.pi) + opd_shift_p[baseline,0],label='SC-FT '+base_list_diff[baseline]+'-P') ax2.plot(time_met[mintime:maxtime], diff_met[baseline,mintime:maxtime], label='MET_SC-MET_FT') plt.ylabel('dOPD base %s-P (mic)' % base_list_diff[baseline]) plt.legend() ax3=plt.subplot(223, sharex=ax1) ax3.plot(time_met[mintime:maxtime], (np.mean(phi_sc_s_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_sc_s[baseline,0] - np.mean(phi_ft_s_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_ft_s[baseline,0])/(2*np.pi) + opd_shift_s[baseline,0]-diff_met[baseline,mintime:maxtime], color='r') plt.ylabel('Residual dOPD-dMET S (mic)') plt.xlabel('Time (s)') ax4=plt.subplot(224, sharex=ax1, sharey=ax3) ax4.plot(time_met[mintime:maxtime], (np.mean(phi_sc_p_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_sc_p[baseline,0] - np.mean(phi_ft_p_int[baseline,mintime:maxtime,:],axis=1)*avg_wave_ft_p[baseline,0])/(2*np.pi) + opd_shift_p[baseline,0]-diff_met[baseline,mintime:maxtime], color='r') plt.ylabel('Residual dOPD-dMET P (mic)') plt.xlabel('Time (s)') if baseline==0: pdf.savefig() # all figures are rather large due to the number of samples and clog the PDF file def lin1met(x, a, b): return (a+x[1])*x[0] + b def fitlambda_color(phi,opd,lambda_met): wave = np.zeros((phi.shape[1],2),dtype='d') opd_shift = np.zeros((phi.shape[1],2),dtype='d') x0 = [2.2,0.0] # initial guess for index in range(0,phi.shape[1]): xdata = [phi[:,index]/(2.0*np.pi),lambda_met] opt1 = optimization.curve_fit(lin1met, xdata, opd, x0) wave[index,:]=[opt1[0][0]+lambda_met,np.sqrt(opt1[1][0,0])] opd_shift[index,:]=[opt1[0][1],np.sqrt(opt1[1][1,1])] return wave,opd_shift #============================================================================== # Determination of the wavelength scales of the SC and FT #============================================================================== print("*** Computation of the wavelength scales of the SC and FT per baseline") opd_sc_s = np.zeros((nbase,phaseftsc.nframe_met),dtype='d') opd_sc_p = np.zeros((nbase,phaseftsc.nframe_met),dtype='d') opd_ft_s = np.zeros((nbase,phaseftsc.nframe_met),dtype='d') opd_ft_p = np.zeros((nbase,phaseftsc.nframe_met),dtype='d') wave_sc_s = np.zeros((nbase,nwave_rebin,2),dtype='d') wave_sc_p = np.zeros((nbase,nwave_rebin,2),dtype='d') wave_ft_s = np.zeros((nbase,phaseftsc.nwave_ft,2),dtype='d') wave_ft_p = np.zeros((nbase,phaseftsc.nwave_ft,2),dtype='d') opd_shift_sc_s = np.zeros((nbase,nwave_rebin,2),dtype='d') opd_shift_sc_p = np.zeros((nbase,nwave_rebin,2),dtype='d') opd_shift_ft_s = np.zeros((nbase,phaseftsc.nwave_ft,2),dtype='d') opd_shift_ft_p = np.zeros((nbase,phaseftsc.nwave_ft,2),dtype='d') for baseline in range(0,nbase): # SC - FT = diff_MET(SC-FT), DIFF order opd_sc_s[baseline,:] = phi_sc_avg_s[baseline,:]*avg_wave_sc_s[baseline,0]/(2.0*np.pi) opd_sc_p[baseline,:] = phi_sc_avg_p[baseline,:]*avg_wave_sc_p[baseline,0]/(2.0*np.pi) opd_ft_s[baseline,:] = phi_ft_avg_s[baseline,:]*avg_wave_ft_s[baseline,0]/(2.0*np.pi) opd_ft_p[baseline,:] = phi_ft_avg_p[baseline,:]*avg_wave_ft_p[baseline,0]/(2.0*np.pi) # opd_sc_s[baseline,:] = diff_met[baseline,:] + phi_ft_avg_s[baseline,:]*avg_wave_ft_s[baseline,0]/(2.0*np.pi) # opd_sc_p[baseline,:] = diff_met[baseline,:] + phi_ft_avg_p[baseline,:]*avg_wave_ft_p[baseline,0]/(2.0*np.pi) # # opd_ft_s[baseline,:] = phi_sc_avg_s[baseline,:]*avg_wave_sc_s[baseline,0]/(2.0*np.pi) - diff_met[baseline,:] # opd_ft_p[baseline,:] = phi_sc_avg_p[baseline,:]*avg_wave_sc_p[baseline,0]/(2.0*np.pi) - diff_met[baseline,:] wave_sc_s[baseline,:,:],opd_shift_sc_s[baseline,:,:] = \ fitlambda_color(phi_sc_s_int[baseline,mintime:maxtime,:], opd_sc_s[baseline,mintime:maxtime], Lambda_met) wave_sc_p[baseline,:,:],opd_shift_sc_p[baseline,:,:] = \ fitlambda_color(phi_sc_p_int[baseline,mintime:maxtime,:], opd_sc_p[baseline,mintime:maxtime], Lambda_met) wave_ft_s[baseline,:,:],opd_shift_ft_s[baseline,:,:] = \ fitlambda_color(phi_ft_s_int[baseline,mintime:maxtime,:], opd_ft_s[baseline,mintime:maxtime], Lambda_met) wave_ft_p[baseline,:,:],opd_shift_ft_p[baseline,:,:] = \ fitlambda_color(phi_ft_p_int[baseline,mintime:maxtime,:], opd_ft_p[baseline,mintime:maxtime], Lambda_met) if False: # Verification of the fit quality as a function of color baseline = 0 plt.figure(figsize=(20,10)) ax1 = plt.subplot(221, sharex=ax1) indices = [0,phaseftsc.nwave_sc//2,phaseftsc.nwave_sc-1] for channel in indices: ax1.plot(time_met[mintime:maxtime],(phi_sc_s_int[baseline,mintime:maxtime,channel]*wave_sc_s[baseline,channel,0]/(2*np.pi)+ opd_shift_sc_s[baseline,channel,0]),label='channel '+str(channel)) ax1.plot(time_met[mintime:maxtime],opd_sc_s[baseline,mintime:maxtime],label='dMET') plt.title('Color dependent fit of SC baseline '+base_list_diff[baseline]+'-S') plt.xlabel('Time (s)') plt.ylabel('OPD (microns)') plt.legend() ax2 = plt.subplot(222, sharex=ax1) # for channel in [50,100,150,200]: # ax2.plot(time_met[mintime:maxtime],(phi_sc_s_int[baseline,mintime:maxtime,channel]*wave_sc_s[baseline,channel,0]/(2*np.pi)+ # opd_shift_sc_s[baseline,channel,0])-opd_sc_s[baseline,mintime:maxtime]) # plt.legend() # plt.ylabel('Residual SC OPD-MET S (mic)') # plt.xlabel('Time (s)') ax2.plot(time_met[mintime:maxtime],((phi_sc_s_int[baseline,mintime:maxtime,phaseftsc.nwave_sc-1]*wave_sc_s[baseline,phaseftsc.nwave_sc-1,0]- \ phi_sc_s_int[baseline,mintime:maxtime,0]*wave_sc_s[baseline,0,0])/(2*np.pi)+ \ opd_shift_sc_s[baseline,phaseftsc.nwave_sc-1,0]-opd_shift_sc_s[baseline,0,0])) plt.legend() plt.ylabel('Diff. SC Channel %i-%i S (mic)'%(phaseftsc.nwave_sc-1,0)) plt.xlabel('Time (s)') ax3 = plt.subplot(223, sharex=ax1) indices = [0,phaseftsc.nwave_ft//2,phaseftsc.nwave_ft-1] for channel in indices: ax3.plot(time_met[mintime:maxtime],(phi_ft_s_int[baseline,mintime:maxtime,channel]*wave_ft_s[baseline,channel,0]/(2*np.pi)+ opd_shift_ft_s[baseline,channel,0]),label='channel '+str(channel)) ax3.plot(time_met[mintime:maxtime],opd_ft_s[baseline,mintime:maxtime],label='dMET') plt.title('Color dependent fit of FT baseline '+base_list_diff[baseline]+'-S') plt.xlabel('Time (s)') plt.ylabel('OPD (microns)') plt.legend() pdf.savefig() ax4 = plt.subplot(224, sharex=ax1) ax4.plot(time_met[mintime:maxtime],((phi_ft_s_int[baseline,mintime:maxtime,phaseftsc.nwave_ft-1]*wave_ft_s[baseline,phaseftsc.nwave_ft-1,0]- \ phi_ft_s_int[baseline,mintime:maxtime,0]*wave_ft_s[baseline,0,0])/(2*np.pi)+ \ opd_shift_ft_s[baseline,phaseftsc.nwave_ft-1,0]-opd_shift_ft_s[baseline,0,0])) plt.ylabel('Diff. FT Channel %i-%i S (mic)'%(phaseftsc.nwave_ft-1,0)) plt.xlabel('Time (s)') plt.legend() pdf.savefig() # dOPD fit constant terms plt.figure(figsize=(20,10)) ax1=plt.subplot(141) for baseline in range(0,nbase): # DETECTOR order ax1.plot(opd_shift_sc_s[order_sc[baseline],:,0],label=base_list_sc[baseline]+"-S") plt.ylabel("SC dOPD fit constant term (microns)") plt.xlabel("SC pixel") plt.margins(0.2) plt.legend() ax2=plt.subplot(142) for baseline in range(0,nbase): # DETECTOR order ax2.plot(opd_shift_sc_p[order_sc[baseline],:,0],label=base_list_sc[baseline]+"-P") plt.ylabel("SC dOPD fit constant term (microns)") plt.xlabel("SC pixel") plt.margins(0.2) plt.legend() ax3=plt.subplot(143) for baseline in range(0,nbase): # DETECTOR order ax3.plot(opd_shift_ft_s[order_ft[baseline],:,0],label=base_list_ft[baseline]+"-S") plt.ylabel("FT dOPD fit constant term (microns)") plt.xlabel("FT pixel") plt.margins(0.2) plt.legend() ax4=plt.subplot(144) for baseline in range(0,nbase): # DETECTOR order ax4.plot(opd_shift_ft_p[order_ft[baseline],:,0],label=base_list_ft[baseline]+"-P") plt.ylabel("FT dOPD fit constant term (microns)") plt.xlabel("FT pixel") plt.margins(0.2) plt.legend() pdf.savefig() #============================================================================== # Write the resulting wavelength tables to a FITS table file #============================================================================== print("*** Saving the resulting wavelength tables in FITS file") now = datetime.datetime.now() prihdr = fits.Header() prihdr['DATE'] = now.strftime("%Y-%m-%dT%H:%M:%S") prihdr['COMMENT'] = "GRAVITY wavelength table from the testbed python code." prihdu = fits.PrimaryHDU(header=prihdr) hdulist = fits.HDUList(prihdu) col_sc_s = [] col_sc_p = [] col_ft_s = [] col_ft_p = [] for baseline in range(0,nbase): # DETECTOR order col_sc_s.append(fits.Column(name=base_list_sc[baseline], format='E', unit='microns', array=wave_sc_s[order_sc[baseline],:,0])) cols_sc_s = fits.ColDefs(col_sc_s) tbhdu_s = fits.BinTableHDU.from_columns(cols_sc_s,name='WAVELENGTH_SC_S') hdulist.append(tbhdu_s) for baseline in range(0,nbase): col_sc_p.append(fits.Column(name=base_list_sc[baseline], format='E', unit='microns', array=wave_sc_p[order_sc[baseline],:,0])) cols_sc_p = fits.ColDefs(col_sc_p) tbhdu_p = fits.BinTableHDU.from_columns(cols_sc_p,name='WAVELENGTH_SC_P') hdulist.append(tbhdu_p) for baseline in range(0,nbase): col_ft_s.append(fits.Column(name=base_list_ft[baseline], format='E', unit='microns', array=wave_ft_s[order_ft[baseline],:,0])) cols_ft_s = fits.ColDefs(col_ft_s) tbhdu_s = fits.BinTableHDU.from_columns(cols_ft_s,name='WAVELENGTH_FT_S') hdulist.append(tbhdu_s) for baseline in range(0,nbase): col_ft_p.append(fits.Column(name=base_list_ft[baseline], format='E', unit='microns', array=wave_ft_p[order_ft[baseline],:,0])) cols_ft_p = fits.ColDefs(col_ft_p) tbhdu_p = fits.BinTableHDU.from_columns(cols_ft_p,name='WAVELENGTH_FT_P') hdulist.append(tbhdu_p) hdulist.writeto(results_dir+filename+'-WAVETABLE.fits',clobber=True) #============================================================================== # Wavelength scale display #============================================================================== avg_wave_sc_s = np.mean(wave_sc_s[:,:,0],axis=0) avg_wave_sc_p = np.mean(wave_sc_p[:,:,0],axis=0) avg_wave_ft_s = np.mean(wave_ft_s[:,:,0],axis=0) avg_wave_ft_p = np.mean(wave_ft_p[:,:,0],axis=0) avg_wave_sc = np.mean([avg_wave_sc_s,avg_wave_sc_p],axis=0) avg_wave_ft = np.mean([avg_wave_ft_s,avg_wave_ft_p],axis=0) # Computation of a model fit on the SC wavelength scale x0 = [0.01, 0.01, 2.2] # initial guess def lin1(x, a, b): return a*x + b def quad(x, a, b, c): return a*x**2 + b*x + c xdata = list(range(0,nwave_rebin)) #opt1 = optimization.curve_fit(lin1, xdata, avg_wave_sc, x0) #[slope_sc_linfit,zp_sc_linfit]=[opt1[0][0],opt1[0][1]] opt1 = optimization.curve_fit(quad, xdata, avg_wave_sc, x0) [a_sc,b_sc,c_sc]=[opt1[0][0],opt1[0][1],opt1[0][2]] xdata = list(range(0,phaseftsc.nwave_ft)) # opt1 = optimization.curve_fit(lin1, xdata, avg_wave_ft, x0) # [slope_ft_linfit,zp_ft_linfit]=[opt1[0][0],opt1[0][1]] opt1 = optimization.curve_fit(quad, xdata, avg_wave_ft, x0) [a_ft,b_ft,c_ft]=[opt1[0][0],opt1[0][1],opt1[0][2]] scalesc = np.array(list(range(0,nwave_rebin))) # modelsc = slope_sc_linfit*scalesc + zp_sc_linfit modelsc = a_sc*scalesc**2 + b_sc*scalesc + c_sc scaleft = np.array(list(range(0,phaseftsc.nwave_ft))) # modelft = slope_ft_linfit*scaleft + zp_ft_linfit modelft = a_ft*scaleft**2 + b_ft*scaleft + c_ft # DIFF > order_sc > SC # SC > order_diff_sc > DIFF # Plots for the SC of the wavelength solution and the residual with respect to a model law in DETECTOR order plt.figure(figsize=(20, 10)) ax1=plt.subplot(221,title='SC wavelength scale polar. S') for baseline in range(0,nbase): # DETECTOR order ax1.plot(list(range(0,nwave_rebin)), wave_sc_s[order_sc[baseline],:,0], label=base_list_sc[baseline]+"-S") ax1.plot(scalesc, modelsc, label='QuadModelSP') plt.ylabel('Wavelength (mic)') plt.legend() ax2=plt.subplot(223, sharex=ax1) for baseline in range(0,nbase): ax2.errorbar(list(range(0,nwave_rebin)),wave_sc_s[order_sc[baseline],:,0]-modelsc, yerr=wave_sc_s[order_sc[baseline],:,1]) plt.xlabel('SC Wavelength channel') plt.ylabel('WaveSC S - Model[avg(WaveSP)] (mic)') plt.margins(0.2) plt.legend() ax3=plt.subplot(222,title='SC wavelength scale polar. P', sharex=ax1, sharey=ax1) for baseline in range(0,nbase): ax3.plot(list(range(0,nwave_rebin)), wave_sc_p[order_sc[baseline],:,0], label=base_list_sc[baseline]+"-P") ax3.plot(scalesc, modelsc, label='QuadModelSP') plt.ylabel('Wavelength (mic)') plt.legend() ax4=plt.subplot(224, sharex=ax1) for baseline in range(0,nbase): ax4.errorbar(list(range(0,nwave_rebin)),wave_sc_p[order_sc[baseline],:,0]-modelsc, yerr=wave_sc_p[order_sc[baseline],:,1]) plt.xlabel('SC Wavelength channel') plt.ylabel('WaveSC P - Model[avg(WaveSP)] (mic)') plt.margins(0.2) plt.legend() pdf.savefig() # Plots for the FT of the wavelength solution and the residual with respect to a model law plt.figure(figsize=(20, 10)) ax1=plt.subplot(221,title='FT wavelength scale polar. S') for baseline in range(0,nbase): ax1.plot(list(range(0,phaseftsc.nwave_ft)), wave_ft_s[order_ft[baseline],:,0], label=base_list_ft[baseline]+"-S") ax1.plot(scaleft, modelft, label='QuadModelSP') plt.ylabel('Wavelength (mic)') plt.legend() ax2=plt.subplot(223, sharex=ax1) for baseline in range(0,nbase): ax2.errorbar(list(range(0,phaseftsc.nwave_ft)),wave_ft_s[order_ft[baseline],:,0]-modelft, yerr=wave_ft_s[order_ft[baseline],:,1]) plt.xlabel('FT Wavelength channel') plt.ylabel('WaveFT S - Model[avg(baselineSP)] (mic)') plt.margins(0.2) plt.legend() ax3=plt.subplot(222,title='FT wavelength scale polar. P', sharex=ax1, sharey=ax1) for baseline in range(0,nbase): ax3.plot(list(range(0,phaseftsc.nwave_ft)), wave_ft_p[order_ft[baseline],:,0], label=base_list_ft[baseline]+"-P") ax3.plot(scaleft, modelft, label='QuadModel') plt.ylabel('Wavelength (mic)') plt.legend() ax4=plt.subplot(224, sharex=ax1) for baseline in range(0,nbase): ax4.errorbar(list(range(0,phaseftsc.nwave_ft)),wave_ft_p[order_ft[baseline],:,0]-modelft, yerr=wave_ft_p[order_ft[baseline],:,1]) plt.xlabel('FT Wavelength channel') plt.ylabel('WaveFT P - Model[avg(baselineSP)] (mic)') plt.margins(0.2) plt.legend() pdf.savefig() print(' Average +/- RMS wavelength shifts with respect to average model scale for SC:') logfile.write(' Average +/- RMS wavelength shifts with respect to average model scale for SC:\n') for baseline in range(0,nbase): # DETECTOR order print(' > %s-S: shift = %+.2f +/- %.2f nm %s-P: shift = %+.2f +/- %.2f nm' % \ (base_list_sc[baseline], 1000.*np.mean(wave_sc_s[order_sc[baseline],:,0]-modelsc), 1000.*np.std(wave_sc_s[order_diff_sc[baseline],:,0]-modelsc), base_list_sc[baseline], 1000.*np.mean(wave_sc_p[order_sc[baseline],:,0]-modelsc), 1000.*np.std(wave_sc_p[order_sc[baseline],:,0]-modelsc))) logfile.write(' > %s-S: shift = %+.2f +/- %.2f nm %s-P: shift = %+.2f +/- %.2f nm\n' % \ (base_list_sc[baseline], 1000.*np.mean(wave_sc_s[order_sc[baseline],:,0]-modelsc), 1000.*np.std(wave_sc_s[order_sc[baseline],:,0]-modelsc), base_list_sc[baseline], 1000.*np.mean(wave_sc_p[order_sc[baseline],:,0]-modelsc), 1000.*np.std(wave_sc_p[order_sc[baseline],:,0]-modelsc))) print(' Average +/- RMS wavelength shifts with respect to average model scale for FT:') logfile.write(' Average +/- RMS wavelength shifts with respect to average model scale for FT:\n') for baseline in range(0,nbase): print(' > %s-S: shift = %+.2f +/- %.2f nm %s-P: shift = %+.2f +/- %.2f nm' % \ (base_list_ft[baseline], 1000.*np.mean(wave_ft_s[order_ft[baseline],:,0]-modelft), 1000.*np.std(wave_ft_s[order_ft[baseline],:,0]-modelft), base_list_ft[baseline], 1000.*np.mean(wave_ft_p[order_ft[baseline],:,0]-modelft), 1000.*np.std(wave_ft_p[order_ft[baseline],:,0]-modelft))) logfile.write(' > %s-S: shift = %+.2f +/- %.2f nm %s-P: shift = %+.2f +/- %.2f nm\n' % \ (base_list_ft[baseline], 1000.*np.mean(wave_ft_s[order_ft[baseline],:,0]-modelft), 1000.*np.std(wave_ft_s[order_ft[baseline],:,0]-modelft), base_list_ft[baseline], 1000.*np.mean(wave_ft_p[order_ft[baseline],:,0]-modelft), 1000.*np.std(wave_ft_p[order_ft[baseline],:,0]-modelft))) #============================================================================== # Comparison with the Argon wavelength scale for the SC (MR and HR only) #============================================================================== if (('MEDIUM' in phaseftsc.header['HIERARCH ESO INS SPEC RES']) or ('HIGH' in phaseftsc.header['HIERARCH ESO INS SPEC RES'])): print("\n*** Comparison with the Argon wavelength scale for the SC") # Starting X values for the spectra in the Argon lamp images # startx_mr = 53 # startx_hr = 38 # # data_path = os.path.dirname(__file__) # if ('MEDIUM' in phaseftsc.header['HIERARCH ESO INS SPEC RES']): # MR case # gravi_ref_wave = fits.open(os.path.join(data_path,'gravity_SC_MR.fits'), mode='readonly') # # gravi_ref_wave = fits.open(os.path.join(data_path,'gravity_SC_MR_20150908.fits'), mode='readonly') # startx = startx_mr # if ('HIGH' in phaseftsc.header['HIERARCH ESO INS SPEC RES']): # HR case # gravi_ref_wave = fits.open(os.path.join(data_path,'gravity_SC_HR.fits'), mode='readonly') # startx = startx_hr # ref_wave_channel = np.zeros((phaseftsc.nregion,phaseftsc.nwave_sc),dtype='d') ref_wave_channel_rebin = np.zeros((phaseftsc.nregion,nwave_rebin),dtype='d') ref_wave_sc_s = np.zeros((nbase,nwave_rebin),dtype='d') ref_wave_sc_p = np.zeros((nbase,nwave_rebin),dtype='d') # for region in range(0,phaseftsc.nregion): # ref_wave_channel[region,:] = gravi_ref_wave['SC_MODEL'].data['W%02i'%region][startx-1:startx-1+phaseftsc.nwave_sc] # first value: phaseftsc.header['HIERARCH ESO DET2 FRAM STRX'] # second value: STARTX value in PROFILE files, table PROFILE_DATA # third value = phaseftsc.nwave_rebin (~235) ref_wave_channel = gravi_fit_argon_mr(910,67,nwave_rebin) # Binning of the HR reference spectra (if necessary) before comparison if 'HIGH' in phaseftsc.header['HIERARCH ESO INS SPEC RES']: # HR case, binning by binwidth print("*** Binning of the reference SC HR scale by a factor %i for comparison" % binwidth) # nwave_rebin = phaseftsc.nwave_sc//binwidth # integer division for region in range(0,phaseftsc.nregion): spectrum = ref_wave_channel[region,:] ref_wave_channel_rebin[region,:] = np.mean(spectrum[:(spectrum.size // binwidth) * binwidth].reshape(-1, binwidth),axis=1) else: ref_wave_channel_rebin = np.copy(ref_wave_channel) # Organization of baselines: # 0 1 2 3 4 5 # base_list_diff=["12", "13", "14", "23", "24", "34"] # base_list_sc= ["14", "34", "13", "24", "12", "23"] for baseline in range(0,nbase): # The baseline signals are sorted as on the SC detector ref_wave_sc_s[baseline,:] = np.mean([ref_wave_channel_rebin[0+8*baseline,:],ref_wave_channel_rebin[2+8*baseline,:], ref_wave_channel_rebin[4+8*baseline,:],ref_wave_channel_rebin[6+8*baseline,:]],axis=0) ref_wave_sc_p[baseline,:] = np.mean([ref_wave_channel_rebin[1+8*baseline,:],ref_wave_channel_rebin[3+8*baseline,:], ref_wave_channel_rebin[5+8*baseline,:],ref_wave_channel_rebin[7+8*baseline,:]],axis=0) residual_argon_s = np.zeros((nbase,nwave_rebin),dtype='d') residual_argon_p = np.zeros((nbase,nwave_rebin),dtype='d') slope_argon_s = np.zeros(nbase,dtype='d') slope_argon_p = np.zeros(nbase,dtype='d') # ---------------------------------- # Correction for dispersion slope # ---------------------------------- for baseline in range(0,nbase): # The baseline signals are sorted as on the SC detector wave_sc_s[baseline,:,0] = wave_sc_s[baseline,:,0]-slope_k*(wave_sc_s[baseline,:,0]-Lambda_met) wave_sc_p[baseline,:,0] = wave_sc_p[baseline,:,0]-slope_k*(wave_sc_p[baseline,:,0]-Lambda_met) wave_ft_s[baseline,:,0] = wave_ft_s[baseline,:,0]-slope_k*(wave_ft_s[baseline,:,0]-Lambda_met) wave_ft_p[baseline,:,0] = wave_ft_p[baseline,:,0]-slope_k*(wave_ft_p[baseline,:,0]-Lambda_met) x0 = [0.1,0.0] for baseline in range(0,nbase): # in DETECTOR order residual_argon_s[baseline,:] = wave_sc_s[order_sc[baseline],:,0] - ref_wave_sc_s[baseline,:] residual_argon_p[baseline,:] = wave_sc_p[order_sc[baseline],:,0] - ref_wave_sc_p[baseline,:] opt1 = optimization.curve_fit(lin1, ref_wave_sc_s[baseline,:], wave_sc_s[order_sc[baseline],:,0], x0) slope_argon_s[baseline] = opt1[0][0] opt1 = optimization.curve_fit(lin1, ref_wave_sc_p[baseline,:], wave_sc_p[order_sc[baseline],:,0], x0) slope_argon_p[baseline] = opt1[0][0] from matplotlib.pyplot import cm color=cm.rainbow(np.linspace(0,1,nbase)) plt.figure(figsize=(20, 10)) ax1=plt.subplot(121) plt.title('Comparison Fourier vs. Argon SC wavelength scales - S polar.') for baseline,c in zip(list(range(0,nbase)),color): # plot in the order of the baselines on the SC detector ax1.errorbar(list(range(0,nwave_rebin)), residual_argon_s[baseline,:], yerr=wave_sc_s[order_sc[baseline],:,1], label=base_list_sc[baseline]+"-S", c=c) plt.xlabel('Wavelength channel') plt.ylabel('Diff. Fourier - Argon (mic)') ax1.set_ylim([-0.005,+0.015]) plt.margins(0.2) plt.legend() ax2=plt.subplot(122, sharex=ax1, sharey=ax1) plt.title('Comparison Fourier vs. Argon SC wavelength scales - P polar.') for baseline,c in zip(list(range(0,nbase)),color): ax2.errorbar(list(range(0,nwave_rebin)), residual_argon_p[baseline,:], yerr=wave_sc_p[order_sc[baseline],:,1], label=base_list_sc[baseline]+"-P", c=c) plt.xlabel('Wavelength channel') plt.ylabel('Diff. Fourier - Argon (mic)') plt.margins(0.2) plt.legend() pdf.savefig() print(' Average +/- RMS wavelength shifts and slopes [Fourier-Argon] for SC:') logfile.write(' Average +/- RMS wavelength shifts and slopes [Fourier-Argon] for SC:\n') for baseline in range(0,nbase): print(' > %s-S: shift = %+.2f +/- %.2f nm slope = %+.2f nm/mic %s-P: shift = %+.2f +/- %.2f nm slope = %+.2f nm/mic' % \ (base_list_sc[baseline], 1000.*np.mean(residual_argon_s[baseline,:]), 1000.*np.std(residual_argon_s[baseline,:]), (slope_argon_s[baseline]-1)*1000., base_list_sc[baseline], 1000.*np.mean(residual_argon_p[baseline,:]), 1000.*np.std(residual_argon_p[baseline,:]), (slope_argon_p[baseline]-1)*1000.)) logfile.write(' > %s-S: shift = %+.2f +/- %.2f nm slope = %+.2f nm/mic %s-P: shift = %+.2f +/- %.2f nm slope = %+.2f nm/mic\n' % \ (base_list_sc[baseline], 1000.*np.mean(residual_argon_s[baseline,:]), 1000.*np.std(residual_argon_s[baseline,:]), (slope_argon_s[baseline]-1)*1000., base_list_sc[baseline], 1000.*np.mean(residual_argon_p[baseline,:]), 1000.*np.std(residual_argon_p[baseline,:]), (slope_argon_p[baseline]-1)*1000.)) print('\n Grand average shift [Fourier-Argon] S-polar. = %+.2f +/- %.2f nm slope = %+.2f +/- %.2f nm/mic' % \ (np.mean(residual_argon_s)*1000.,np.std(residual_argon_s)*1000., (np.mean(slope_argon_s)-1)*1000.,np.std(slope_argon_s)*1000.)) print(' Grand average shift [Fourier-Argon] P-polar. = %+.2f +/- %.2f nm slope = %+.2f +/- %.2f nm/mic' % \ (np.mean(residual_argon_p)*1000.,np.std(residual_argon_p)*1000., (np.mean(slope_argon_p)-1)*1000.,np.std(slope_argon_p)*1000.)) logfile.write('\n Grand average shift [Fourier-Argon] S-polar. = %+.2f +/- %.2f nm slope = %+.2f +/- %.2f nm/mic\n' % \ (np.mean(residual_argon_s)*1000.,np.std(residual_argon_s)*1000., (np.mean(slope_argon_s)-1)*1000.,np.std(slope_argon_s)*1000.)) logfile.write(' Grand average shift [Fourier-Argon] P-polar. = %+.2f +/- %.2f nm slope = %+.2f +/- %.2f nm/mic\n' % \ (np.mean(residual_argon_p)*1000.,np.std(residual_argon_p)*1000., (np.mean(slope_argon_p)-1)*1000.,np.std(slope_argon_p)*1000.)) print("-"*80) pdf.close() logfile.close() #============================================================================== # MAIN PROGRAM #============================================================================== if __name__ == '__main__': # if len(sys.argv) == 3 : # filename=sys.argv[1] # filename2=sys.argv[2] # # if filename == '' : # filename=raw_input(" Enter PREPROC file name (without .fits) : ") # if filename2 == '' : # filename2=raw_input(" Enter PHSE_SC_FT file name (without .fits) : ") # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-06-01T10-22-48-MR', # 'phase_sc_ft.2015-06-01T10-22-48-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=0, frame_shift=1, binwidth=1) # # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-08-06T15-36-59-MR', # 'phase_sc_ft.2015-08-06T15-36-59-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=0, frame_shift=1, binwidth=1) # # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-08-14T07-13-40-MR', # 'phase_sc_ft.2015-08-14T07-13-40-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=0, frame_shift=1, binwidth=1) # # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-08-16T19-28-51-MR', # 'phase_sc_ft.2015-08-16T19-28-51-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=0, frame_shift=0, binwidth=1) # # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-08-18T15-23-50-MR', # 'phase_sc_ft.2015-08-18T15-23-50-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=0, frame_shift=0, binwidth=1) # # # High spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-08-15T08-22-03-HR', # 'phase_sc_ft.2015-08-15T08-22-03-HR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=0, frame_shift=1, binwidth=5) # # # Low spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-08-14T06-40-05-LR', # 'phase_sc_ft.2015-08-14T06-40-05-LR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=0, frame_shift=1, binwidth=1) # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-09-03T13-05-25-MR', # 'phase_sc_ft.2015-09-03T13-05-25-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=1, frame_shift=0, binwidth=1) # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE # gravi_wavelength('gravi_preproc_GRAVITY.2015-09-11T04-34-38-MR', # 'phase_sc_ft.2015-09-11T04-34-38-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=1, frame_shift=0, binwidth=1) # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE WITH 1mm FIBER STRECHING ON TEL 1 # gravi_wavelength('gravi_preproc_GRAVITY.2015-09-16T00-50-06-MR', # 'phase_sc_ft.2015-09-16T00-50-06-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=1, frame_shift=0, binwidth=1) # # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE WITH 1mm FIBER STRECHING ON TEL 1 # gravi_wavelength('gravi_preproc_GRAVITY.2015-09-16T03-13-46-MR', # 'phase_sc_ft.2015-09-16T03-13-46-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=1, frame_shift=0, binwidth=1) # NEW SINGLE PHASEFTSC file processing # Medium spectral resolution UNCOMMENT THE FOLLOWING LINES TO ACTIVATE WITH 1mm FIBER STRECHING ON TEL 1 # gravi_wavelength('phase_sc_ft.2015-09-16T03-13-46-MR', delay_ft, boxsmooth_ft, boxsmooth_met, # fiber_coupler_only, window_sc=1, frame_shift=0, binwidth=1) # Medium spectral resolution with scanning problem gravi_wavelength('phase_sc_ft.2016-03-27-MR', delay_ft, boxsmooth_ft, boxsmooth_met, fiber_coupler_only, window_sc=1, frame_shift=0, binwidth=1) #============================================================================== # Masking of the segments when the FDDLs are moving #============================================================================== # tmin = np.zeros(phaseftsc.nframe_sc) # tmax = np.zeros(phaseftsc.nframe_sc) # mask_d = np.zeros(phaseftsc.nframe_met) # for tick_sc in range(0,phaseftsc.nframe_sc-1): # tmin[tick_sc] = phaseftsc.time_sc[tick_sc] + dit_sc/2.0 # tmax[tick_sc] = phaseftsc.time_sc[tick_sc+1] - dit_sc/2.0 # for i in range(int(tmin[tick_sc]//dt_met),int(tmax[tick_sc]//dt_met)): # mask_d[i] = 1 # d = np.where(mask_d == 0) # #print d # plt.figure() # t = time_met[d] # r = np.squeeze(diff_met[0,d]) # print t.shape # print r.shape # plt.plot(t,r) #