Added Strehl ratio fits

hrm-hc/MAIN.py

@@ -17,7 +17,7 @@ start_time = time.time()
 #HARMONI_DIR = '/Users/avigan/data/HARMONI-HC/'
 #COM_DIR     = '/Users/avigan/data/HARMONI-HC/COMMON_FILES/'
 HARMONI_DIR = '/Users/carlotal/harmoni-hc/hrm-hc'
-COM_DIR     = '/Users/carlotal/harmoni-hc/hrm-hc/COMMON_FILES/'
+COM_DIR     = '/Users/carlotal/harmoni-hc/hrm-hc/common_files/'
 
 # Artifical dispersion corrector (dispersion is gamma times what it should be)
 gamma = 1#%0.1;%1;
@@ -91,7 +91,7 @@ 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)
+BAND = 'TEST_K'  # %Choose between the bands defined hereafter (use a valid string)
 
 ### Bands ###
 
@@ -575,6 +575,8 @@ df = np.abs(f1D[1]-f1D[0])*fudge_df
 # |_____/|_____|_|  |_|\____/|______/_/    \_\_|  |_____|_| \_|\_____|     |______/_/ \_\_|     \____/|_____/ \____/|_|  \_\______|_____/ 
 #                                                                                                                                                                                                                                                                                 
 
+SR = np.zeros((len(HA_vect),N_LD))
+
 log.info('Computing the different exposures...')
 for i in range(len(HA_vect)):
     try:
@@ -648,6 +650,9 @@ for i in range(len(HA_vect)):
             WF_CP_0 = np.squeeze(WF_CP_0)
             WF_NCPA_0 = np.squeeze(WF_NCPA_0)
             WF_0=(WF_CP_0-WF_corr_0+EELT_phase_0+WF_NCPA_0)
+            
+            DSP_V,DSP_I = DSP(WF_0,BigD/D,np.transpose(ORIGINAL_PUPIL),0,1000)
+            log.info('HA={0}: Aberration level (0-1000 lambda/D) at {1} nm: {2} nm'.format(HA_0, lambda_vect[k]*1e9, 1e9*DSP_V))
 
             #DSP_V, DSP_I = DSP(WF_CP_0+EELT_phase_0+WF_NCPA_0, 1132/1024, PUPILLE_0, 2, 30)
             #log.info('HA={0}: Aberration level pre correction (2-30 lambda/D) at {1} nm: {2} nm'.format(HA_0, lambda_vect[k]*1e9, 1e9*DSP_V))
@@ -699,10 +704,10 @@ for i in range(len(HA_vect)):
             
             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
             FTO_tel_0_ForSR = ft(I_s_0_ForSR,1,2*OWA,2*OWA,DiamFTO,DiamFTO,1,len(dsp),len(dsp))
-            SR = np.sum(np.abs(FTO_tel_0_ForSR*FTO_atm))/np.sum(np.abs(FTO_tel_0_ForSR))
+            SR[i,k] = np.sum(np.abs(FTO_tel_0_ForSR*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: {:4.2f}'.format(SR))
+            log.info(' => Strehl at {:4.2f}nm: {:4.2f}'.format(lambda_vect[k]*1e9,SR[i,k]))
 
             #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
@@ -804,6 +809,8 @@ for i in range(len(HA_vect)):
     fits.setval(image_filename_nomask, 'APODIZER', value=APODIZER)
     fits.setval(image_filename_nomask, 'FPM', value=FPM)
 
+SR_filename = path_directory + 'SR.fits'.format(i)
+fits.writeto(SR_filename, SR)
 log.info('Total time: {0:.2f} min'.format((time.time()-start_time)/60))
 
 #save(['DATA_Aberrations_CONFIG=' CONFIG '_DEC=' num2str(DEC) 'deg_TimeBetweenExp=' num2str(TIMEBETWEENEXPOSURES) 'min_MeanTimeAfterMeridian=' num2str(HA_offset) 'h.dat'], 'DATA_AB')