try: import pyfits except: from astropy.io import fits as pyfits import numpy as np import matplotlib.pyplot as plt import os def smooth(y,n): w = np.ones(n) return np.convolve (w/w.sum(),y,mode='same') def wrap(y): return np.angle (np.exp (1.j * y)) #file = "phase_sc_ft_normal.fits" #file = "phase_sc_ft_1er.fits" #file = "phase_sc_ft_200um_1period.fits" file= "/data/1/2016-12-09/reduced/GRAVITY.2016-12-09T15-49-03_debug.fits" file= "/data/1/2016-12-09/reduced/GRAVITY.2016-12-09T15-41-52_debug.fits" file= "/data/1/2016-12-10/reduced/GRAVITY.2016-12-10T03-32-11_debug.fits" file= "/data/1/2016-12-10/reduced_reference/GRAVITY.2016-12-10T03-32-11_debug.fits" # file medium comb with pb # ********************* # 2016-12-11 # # LOW SPLIT # OK #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T05-40-41_debug.fits" # OK #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T05-33-33_debug.fits" # OK # MED SPLIT OK #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T06-20-18_debug.fits" # OK file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T06-13-09_debug.fits" # OK # HIGH SPLIT #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T06-55-41_debug.fits" # ~BAD #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T06-50-34_debug.fits" # ~BAD # HIGH COMB BAD #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T07-31-35_debug.fits" # BAD #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T07-26-27_debug.fits" # BAD # MED COMB BAD BAD file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T08-08-41_debug.fits" # BAD file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T08-02-39_debug.fits" # BAD # LOW COMB #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T08-42-16_debug.fits" # OK #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T08-36-14_debug.fits" # OK # 2016-10-14 #file="/data/1/ESO_archive_night/2016-11-13/reduced/GRAVI.2016-11-14T11:15:43.004_debug.fits" file="/data/1/ESO_archive_night/2016-11-13/reduced/GRAVI.2016-11-14T11:10:03.985_debug.fits" # small instability on all tel #file="/data/1/ESO_archive_night/2016-11-13/GRAVI.2016-11-14T09:30:03.151_debug.fits" # ********************* # 2016-09-22 # # MED COMB #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T09-26-38_debug.fits" # BAD #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T09-32-24_debug.fits" #file= "/data/1/2016-12-11/reduced/GRAVITY.2016-12-11T08-02-39_debug.fits" # BAD # LOW COMB #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T06-25-23_debug.fits" #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T06-30-54_debug.fits" # LOW SPLIT #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T07-03-08_debug.fits" #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T07-09-09_debug.fits" # MED SPLIT #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T07-37-19_debug.fits" #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T07-43-38_debug.fits" # HIGH SPLIT #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T08-11-52_debug.fits" #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T08-17-03_debug.fits" # HIGH COMB #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T08-53-01_debug.fits" #file= "/data/1/P2VMs_September/reduced/GRAVITY.2016-09-22T08-57-46_debug.fits" file= "/data/1/2016-12-14/GRAVITY.2016-12-14T20-26-13_debug.fits" file= "/data/1/2016-12-15/GRAVITY.2016-12-15T11-33-52_debug.fits" file= "/data/1/2016-12-15/GRAVITY.2016-12-15T11-28-39_debug.fits" def plot_debug (file): hdulist = pyfits.open(file) tel = [[0,1],[0,2],[0,3],[1,2],[1,3],[2,3]] region_sc=np.array([[8,9,10,11],[20,21,22,23],[4,5,6,7],[16,17,18,19],[0,1,2,3],[12,13,14,15]]) region_ft=np.array([[20,21,22,23],[12,13,14,15],[16,17,18,19],[4,5,6,7],[8,9,10,11],[0,1,2,3]]) if hdulist['IMAGING_DETECTOR_SC'].data.size == 48: region_sc=region_sc*2 if hdulist['IMAGING_DETECTOR_FT'].data.size == 48: region_ft=region_ft*2 lbd_met = 1.908 lbd_sc = 2.24 lbd_ft = 2.2 tpi = np.pi * 2 plt.clf() phi=[] plt.figure(10) plt.clf() plt.figure(11) plt.clf() for base in range(6): plt.figure(base) plt.clf() opd_sc = hdulist['OPD_FT'].data['OPD_SC'][base::6] # time_sc = hdulist['OPD_FT'].data['TIME'][base::6] * 1e-6 opd_ft = hdulist['OPD_FT'].data['OPD'][base::6] time_ft = hdulist['OPD_FT'].data['TIME'][base::6] * 1e-6 opd_met = hdulist['OPD_FT'].data['PHASE_MET_FC'][base::6] / tpi * lbd_met * 1e-6 # time_met = hdulist['OPD_FT'].data['TIME'][tel[base][1]::4] * 1e-6 plt.plot (time_ft, opd_sc, 'r,', label="OPD SC") plt.plot (time_ft, opd_ft, 'g-', label="OPD FT") plt.plot (time_ft, opd_met-opd_met.mean(), 'b,', label="OPD MET") # opd_met = np.interp (time_ft, time_met, smooth(opd_met,10)) #opd_ft = np.interp (time_sc, time_ft, smooth(opd_ft,10)) # opd_sc = np.interp (time_ft, time_sc, opd_sc) #plt.plot (time_sc, opd_sc, 'r-') #plt.plot (time_sc, opd_met_sc, 'r--') #plt.plot (time_sc, opd_ft, 'b-') #plt.plot (time_sc, opd_met_ft, 'b--') #plt.plot (time_sc, opd_met, 'g-') # residus=(opd_sc - opd_met - opd_ft)[opd_sc != 0] plt.plot (time_ft[opd_sc != 0], residus-residus.mean(), 'y,', label="Residus") plt.legend() print(np.std(residus)) phi.append(residus-np.mean(residus)) #plt.show() #plt.plot (time_sc, wrap (opd_sc/lbd_sc*tpi)) #plt.plot (time_sc, opd_ft/lbd_ft*tpi) #plt.plot (time_sc, opd_met/lbd_met*tpi) # opd_met = np.interp (time_sc, time_met, opd_met) # opd_ft = np.interp (time_sc, time_ft, opd_ft) # Residual phases phi_a = np.angle (np.exp (tpi*1.j * (residus)/lbd_sc)) print(np.std(phi_a)) #phi.append(phi_a-np.mean(phi_a)) # phi_b = np.angle (np.exp (tpi*1.j * (opd_sc - opd_met_sc)/lbd_sc)) # plt.plot (time_sc, phi_a-np.mean(phi_a), "-") # plt.plot (time_sc, phi_b, "g") # plt.show() # plt.plot (tphase * 1e6, phi_unwrap) plt.figure(10) spectrum_SC = hdulist['SPECTRUM_DATA_SC'] spectrum_sum=np.zeros((spectrum_SC.data["DATA"+str(region_ft[base][0]+1)].shape[0])) for out in range(4) : spectrum_sum+=spectrum_SC.data["DATA"+str(region_sc[base][out]+1)].sum(axis=1) plt.plot(spectrum_sum, label="Flux base"+str(tel[base][0]+1)+str(tel[base][1]+1)) plt.figure(11) spectrum_FT = hdulist['SPECTRUM_DATA_FT'] spectrum_sum=np.zeros((spectrum_FT.data["DATA"+str(region_ft[base][0]+1)].shape[0])) for out in range(4) : spectrum_sum+=spectrum_FT.data["DATA"+str(region_ft[base][out]+1)].sum(axis=1) plt.plot(spectrum_sum, label="Flux base"+str(tel[base][0]+1)+str(tel[base][1]+1)) plt.legend() hdulist.close() plot_debug (file)