small changes to MAIN.py

hrm-hc/MAIN.py

@@ -87,8 +87,8 @@ Mid_HA = 0.2 #%0%30/60;
 ### Observation time & exposure time are constrained by dsp data #############
 ### !!! DO NOT CHANGE THAT FOR THE MOMENT !!! ################################
 
-#Total observation time [hours] ; must be 2h to correctly use dsp data 
-T_obs = 2#2 #;%1/60;%30/60;%1/60
+#Total observation time [hours] ; must be 2h or 4h to correctly use dsp data 
+T_obs = 2# 4 #2 #;%1/60;%30/60;%1/60
 
 #Exposure time of individual exposures [hours] 
 T_exp = 1/60
@@ -97,7 +97,13 @@ T_exp = 1/60
 N_exp = np.int(T_obs/T_exp)# %60 exposures per hour = 1 per minute
 ##############################################################################
 
-BAND = 'TEST_H'  # %Choose between the bands defined hereafter (use a valid string)
+# FOV %%%
+Min_wavelength_for_Nyquist = 1450e-9#*1.2
+MASperPIXEL                = LD2MAS(0.5, 38.542, Min_wavelength_for_Nyquist)
+N_im                       = 215
+FOV                        = N_im * MASperPIXEL
+
+BAND = 'TEST_Nyquist'  # %Choose between the bands defined hereafter (use a valid string)
 
 ### Bands ###
 
@@ -140,6 +146,28 @@ if BAND == 'TEST_H':
             MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
             MASK_OFFSET = 0  # mas
             FPM = 'FPM_H'
+elif BAND == 'TEST_Nyquist':
+    N_LD = 3
+    LBD_min = H_band[0]*1e-6
+    LBD_max = 1.45e-6
+    if APODIZER == 'SP_1':
+        if ADC_flag == 0:
+            MASK = fits.getdata(COM_DIR + 'FPM_H_SP1_NoADC.fits')
+            FPM = 'FPM_H'
+        else:
+            MASK = fits.getdata(COM_DIR + 'FPM_H_SP1_ADC.fits')
+            MASK_OFFSET = 0  # mas
+            FPM = 'FPM_H'
+    else:
+        if ADC_flag == 0:
+            MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_NoADC.fits')
+            FPM = 'FPM_H'
+        else:
+            # MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_ADC.fits')
+            # MASK_OFFSET = -7 #mas
+            MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
+            MASK_OFFSET = 0  # mas
+            FPM = 'FPM_H'
 elif BAND == 'TEST_K':
     N_LD = 3
     LBD_min = K_band[0]*1e-6
@@ -410,7 +438,6 @@ else:
 log.info('Chosen configuration: ' + CONFIG)
 log.info('Required time for ADI at First HA, Mid HA, and Last HA: {:05.2f}-{:05.2f}-{:05.2f} min' .format(ADI_time_FIRST, ADI_time_MID, ADI_time_LAST))
 
-
 log.info('Lambda vector: {:05.1f}nm - {:05.1f}nm ; Number of spectral elements={:5.0f}'.format(np.round(np.min(lambda_vect)*1e9), np.round(np.max(lambda_vect)*1e9), N_LD))
 
 EELT_phase = eeltphase_flag*fits.getdata(COM_DIR + 'M1phase_1132_1024D38.542.fits')
@@ -438,12 +465,6 @@ log.info('FPM: ' + FPM)  # this parameter is set as a function of the apodizer &
 #SP    = A_big
 SP = np.rot90(A)
 
-# FOV %%%
-Min_wavelength_for_Nyquist = 1450e-9
-MASperPIXEL                = LD2MAS(0.5, 38.542, Min_wavelength_for_Nyquist)
-N_im                       = 215
-FOV                        = N_im * MASperPIXEL
-
 # Focal Plane Mask transmission & detector geometry
 #MASK = 1-((1-MASK)*(1-1e-4)) #not necessary anymore (hardcoded in FPM)
 x_mask, X_mask, Y_mask, R_mask, T_mask = VECT(len(MASK), FOV)
@@ -584,7 +605,11 @@ fits.writeto(path_directory + 'WAVELENGTH_VECT.fits', lambda_vect, overwrite=Tru
 # |_|  |_/_/    \_\______\____/ 
 #                                                             
 
-dsp = fits.getdata(COM_DIR + seeing + '_Res-SpatialPSD_Scaled.fits')
+if T_obs == 2:
+    dsp = fits.getdata(COM_DIR + seeing + '_Res-SpatialPSD_Scaled.fits')
+elif T_obs == 4:
+    dsp = fits.getdata(COM_DIR + seeing + '_4h_Res-SpatialPSD_Scaled.fits')
+    
 dsp = dsp*DSP_gamma**2
 size_dsp = len(dsp[0])
 dsp_array= dsp
@@ -604,7 +629,7 @@ df = np.abs(f1D[1]-f1D[0])*fudge_df
 # |_____/|_____|_|  |_|\____/|______/_/    \_\_|  |_____|_| \_|\_____|     |______/_/ \_\_|     \____/|_____/ \____/|_|  \_\______|_____/ 
 #                                                                                                                                                                                                                                                                                 
 
-SR = np.zeros((len(HA_vect), N_LD))
+SR = np.zeros((len(HA_vect), N_LD,2))
 SHIFT_MASK = np.zeros((len(HA_vect), N_LD))
 SHIFT_PSF  = np.zeros((len(HA_vect), N_LD))
 
@@ -733,14 +758,19 @@ for i in range(len(HA_vect)):
 
             I_s_0 = np.abs(ft(np.transpose(SP)*MISSSEG*np.exp(2*1j*pi*WF_0/lambda_vect[k])*np.exp(2*1j*pi*alpha_disp_0*Y), 1, BigD/D, BigD/D, 2*OWA, 2*OWA, 1, len(dsp), len(dsp)))**2
             FTO_tel_0 = ft(I_s_0, 1, 2*OWA, 2*OWA, DiamFTO, DiamFTO, 1, len(dsp), len(dsp))
-            I_s_post_0 = np.abs(ft(FTO_tel_0*FTO_atm, 1, DiamFTO, DiamFTO, 2*OWA, 2*OWA, 1, N_im, N_im))
-            
+            I_s_post_0 = np.abs(ft(FTO_tel_0*FTO_atm, 1, DiamFTO, DiamFTO, 2*OWA, 2*OWA, 1, N_im, N_im))        
             I_s_0_ForSR = np.abs(ft(np.transpose(ORIGINAL_PUPIL), 1, BigD/D, BigD/D, 2*OWA, 2*OWA, 1, len(dsp), len(dsp)))**2
+            I_s_0_ForSR_2 = np.abs(ft(MISSSEG*np.exp(2*1j*pi*WF_0/lambda_vect[k])*np.exp(2*1j*pi*alpha_disp_0*Y), 1, BigD/D, BigD/D, 2*OWA, 2*OWA, 1, len(dsp), len(dsp)))**2
+            
             FTO_tel_0_ForSR = ft(I_s_0_ForSR, 1, 2*OWA, 2*OWA, DiamFTO, DiamFTO, 1, len(dsp), len(dsp))
-            SR[i, k] = np.sum(np.abs(FTO_tel_0_ForSR*FTO_atm))/np.sum(np.abs(FTO_tel_0_ForSR))
+            FTO_tel_0_ForSR_2 = ft(I_s_0_ForSR_2, 1, 2*OWA, 2*OWA, DiamFTO, DiamFTO, 1, len(dsp), len(dsp))
+            
+            SR[i, k, 0] = np.sum(np.abs(FTO_tel_0_ForSR*FTO_atm))/np.sum(np.abs(FTO_tel_0_ForSR))
+            SR[i, k, 1] = np.sum(np.abs(FTO_tel_0_ForSR_2*FTO_atm))/np.sum(np.abs(FTO_tel_0_ForSR))
             
             log.info('Wavelength {:4d}/{:4d} - {:.2f} min left'.format(k+1, N_LD, t_left))
-            log.info(' => Strehl at {:4.2f}nm: {:4.2f}'.format(lambda_vect[k]*1e9, SR[i, k]))
+            log.info(' => Strehl ideal telescope at {:4.2f}nm: {:4.2f}'.format(lambda_vect[k]*1e9, SR[i, k, 0]))
+            log.info(' => Strehl tel+instrument at {:4.2f}nm: {:4.2f}'.format(lambda_vect[k]*1e9, SR[i, k, 1]))
 
             #Normalization of I_s_post_0 to reflect a unit energy arriving on
             #the pupil, i.e., we just have to multiply the PSF by the number of