Skip to content

GitLab

Commit 70831f39 authored by Arthur Vigan's avatar Arthur Vigan
Browse files

Save the PSF and MASK shift, minor code changes 

The PSF and MASK shifts are necessary for the HARMONI-HC photometry
code because fake planet injection is performed *before* masking the
data with the FPM. Otherwise the signal of the planet is not masked
properly in the data.
parent 9b86cbe3
Hide whitespace changes
Inline Side-by-side
hrm-hc/MAIN.py
View file @ 70831f39
......@@ -604,7 +604,9 @@ 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))
SHIFT_MASK = np.zeros((len(HA_vect), N_LD))
SHIFT_PSF = np.zeros((len(HA_vect), N_LD))
log.info('Computing the different exposures...')
for i in range(len(HA_vect)):
......@@ -645,7 +647,7 @@ for i in range(len(HA_vect)):
WF_opt_sensing_0 = np.squeeze(WF_opt_sensing_0)
WF_opt_sensing_NCPA_0 = np.squeeze(WF_opt_sensing_NCPA_0)
WF_corr_0=np.real(M4Magic(COM_DIR,EELT_phase_0+WF_opt_sensing_0, PUPILLE_0, -pix_shift_0[0], -pix_shift_0[1]))
WF_corr_0 = np.real(M4Magic(COM_DIR, EELT_phase_0+WF_opt_sensing_0, PUPILLE_0, -pix_shift_0[0], -pix_shift_0[1]))
#DSP_V, DSP_I = DSP(EELT_phase_0+WF_opt_sensing_0-WF_corr_0+WF_opt_sensing_NCPA_0, 1132/1024, PUPILLE_0, 2, 30)
#log.info('Residual aberrations at HA={0} h: {1} nm'.format(HA_0, 1e9*DSP_V))
......@@ -678,9 +680,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)
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)
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)
......@@ -692,10 +694,12 @@ for i in range(len(HA_vect)):
if ADC_flag == 1:
#alpha_disp_0 = MAS2LD(1000*DISPER_ADC(1.45,lambda_vect[k]*1e6,90-ELEVATION(DEC,HA_0),z_min,z_max),38.542,lambda_vect[k])-MAS2LD(1000*DISPER_ADC(1.45,lambda_vect[-1]*1e6,90-ELEVATION(DEC,HA_0),z_min,z_max),38.542,lambda_vect[k])
alpha_disp_0 = MAS2LD(1000*DISPER_ADC(lambda_vect[0]*1e6,lambda_vect[k]*1e6,90-ELEVATION(DEC,HA_0),z_min,z_max),38.542,lambda_vect[k])
alpha_disp_0 = MAS2LD(1000*DISPER_ADC(lambda_vect[0]*1e6, lambda_vect[k]*1e6, 90-ELEVATION(DEC, HA_0), z_min, z_max),
38.542, lambda_vect[k])
else:
#alpha_disp_0 = MAS2LD(1000*DISPER(1.45,lambda_vect[k]*1e6,90-ELEVATION(DEC,HA_0)),38.542,lambda_vect[k])-MAS2LD(1000*DISPER(1.45,lambda_vect[-1]*1e6,90-ELEVATION(DEC,HA_0)),38.542,lambda_vect[k])
alpha_disp_0 = MAS2LD(1000*DISPER(lambda_vect[0]*1e6,lambda_vect[k]*1e6,90-ELEVATION(DEC,HA_0)),38.542,lambda_vect[k])
alpha_disp_0 = MAS2LD(1000*DISPER(lambda_vect[0]*1e6, lambda_vect[k]*1e6, 90-ELEVATION(DEC, HA_0)),
38.542, lambda_vect[k])
#log.info('Differential dispersion angle with 2.45um: {0}'.format(alpha_disp_0))
x, X, Y, R, T = VECT(len(SP), 1132/1024)
......@@ -703,17 +707,17 @@ for i in range(len(HA_vect)):
E_0 = 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, N_im, N_im)
#Commented to save time ; halo only
I_0=np.abs(E_0)**2
I_0 = np.abs(E_0)**2
#HALO
FT2DSP_norm=1e-18*(2*pi/lambda_vect[k])**2
FT2DSP_norm = 1e-18*(2*pi/lambda_vect[k])**2
# Normalisation by maximum of phase autocorrelation (in 0) = phase variance
# then PHASE STRUCTURE FUNCTION (Roddier) = 2*(autocol(0, 0) - autocol(u, v))
# then Atmospheric Optical Transfert function
FT_PST = RENORM_FT_DSP*FT2DSP_norm
#Rotation of FT_PST to simulate changing wind direction
FT_PST = rotate(FT_PST, Wind_Angle, order=1, reshape = False)
FT_PST = rotate(FT_PST, Wind_Angle, order=1, reshape=False)
PH_STR = 2*(FT_PST[int(size_dsp//2+1), int(size_dsp//2+1)]-FT_PST)
......@@ -731,12 +735,12 @@ for i in range(len(HA_vect)):
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_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[i,k] = np.sum(np.abs(FTO_tel_0_ForSR*FTO_atm))/np.sum(np.abs(FTO_tel_0_ForSR))
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[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 at {:4.2f}nm: {:4.2f}'.format(lambda_vect[k]*1e9,SR[i,k]))
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
......@@ -762,8 +766,8 @@ for i in range(len(HA_vect)):
# alpha_shift = -sum(I1D*V1D)
# shift_value[k] = alpha_shift/(x_mask[1]-x_mask[0])
I_temp[k, :, :] = I_s_post_0
I_temp_nomask[k, :, :] = I_s_post_0
I_temp[k] = I_s_post_0
I_temp_nomask[k] = I_s_post_0
# if exhaustive_writing_flag==1:
# fits.writeto(path_directory + 'PSF_HALO_OFF_Nexp{0}_Nlbd{1}.fits'.format(i, k), I_0)
......@@ -776,9 +780,9 @@ for i in range(len(HA_vect)):
I_temp = APPLY_CROSSTALK(I_temp)
if ADC_flag == 1:
shift_pix_FPM = 0.5*(LD2MAS(IWA_SP,38.542,lambda_vect[N_LD-1])-LD2MAS(IWA_SP,38.542,lambda_vect[0]) + 1000*DISPER_ADC(lambda_vect[0]*1e6,lambda_vect[N_LD-1]*1e6,z_max,z_min,z_max))
DISPER_MAX = DISPER_ADC(lambda_vect[0]*1e6,lambda_vect[N_LD-1]*1e6,z_max,z_min,z_max)
DISPER_Z = DISPER_ADC(lambda_vect[0]*1e6,lambda_vect[N_LD-1]*1e6,90-ELEVATION(DEC,HA_0),z_min,z_max)
shift_pix_FPM = 0.5*(LD2MAS(IWA_SP, 38.542, lambda_vect[N_LD-1])-LD2MAS(IWA_SP, 38.542, lambda_vect[0]) + 1000*DISPER_ADC(lambda_vect[0]*1e6, lambda_vect[N_LD-1]*1e6, z_max, z_min, z_max))
DISPER_MAX = DISPER_ADC(lambda_vect[0]*1e6, lambda_vect[N_LD-1]*1e6, z_max, z_min, z_max)
DISPER_Z = DISPER_ADC(lambda_vect[0]*1e6, lambda_vect[N_LD-1]*1e6, 90-ELEVATION(DEC, HA_0), z_min, z_max)
shift_pix_FPM *= (DISPER_MAX-DISPER_Z)/DISPER_MAX
# Mask offset ; non-zero when mask is *not* the one that should be ideally be used (required because of limited number of slots in FPM wheel)
......@@ -787,59 +791,66 @@ for i in range(len(HA_vect)):
# Conversion from mas to pixel
shift_pix_FPM /= MASperPIXEL
else:
shift_pix_FPM = 0.5*(LD2MAS(IWA_SP,38.542,lambda_vect[N_LD-1])-LD2MAS(IWA_SP,38.542,lambda_vect[0]) + 1000*DISPER(lambda_vect[0]*1e6,lambda_vect[N_LD-1]*1e6,z_max))
DISPER_MAX = DISPER(lambda_vect[0]*1e6,lambda_vect[N_LD-1]*1e6,z_max)
DISPER_Z = DISPER(lambda_vect[0]*1e6,lambda_vect[N_LD-1]*1e6,90-ELEVATION(DEC,HA_0))
shift_pix_FPM = 0.5*(LD2MAS(IWA_SP, 38.542, lambda_vect[N_LD-1])-LD2MAS(IWA_SP, 38.542, lambda_vect[0]) + 1000*DISPER(lambda_vect[0]*1e6, lambda_vect[N_LD-1]*1e6, z_max))
DISPER_MAX = DISPER(lambda_vect[0]*1e6, lambda_vect[N_LD-1]*1e6, z_max)
DISPER_Z = DISPER(lambda_vect[0]*1e6, lambda_vect[N_LD-1]*1e6, 90-ELEVATION(DEC, HA_0))
shift_pix_FPM *= (DISPER_MAX-DISPER_Z)/DISPER_MAX
shift_pix_FPM /= MASperPIXEL
#print('FPM offset: {:f}'.format(shift_pix_FPM))
#MASK_temp = MASK
MASK_temp = shift(MASK, [-shift_pix_FPM,0], output=None, order=1, mode='constant', cval=0.0, prefilter=True)
MASK_temp = shift(MASK, [-shift_pix_FPM, 0], output=None, order=1, mode='constant', cval=0.0, prefilter=True)
for k in range(N_LD):
for k in range(N_LD):
if ADC_flag == 1:
shift_value = 1000*DISPER_ADC(lambda_vect[0]*1e6,lambda_vect[k]*1e6,90-ELEVATION(DEC,HA_0),z_min,z_max)
shift_value = 1000*DISPER_ADC(lambda_vect[0]*1e6, lambda_vect[k]*1e6, 90-ELEVATION(DEC, HA_0), z_min, z_max)
shift_value /= MASperPIXEL
else:
shift_value = 1000*DISPER(lambda_vect[0]*1e6,lambda_vect[k]*1e6,90-ELEVATION(DEC,HA_0))
shift_value = 1000*DISPER(lambda_vect[0]*1e6, lambda_vect[k]*1e6, 90-ELEVATION(DEC, HA_0))
shift_value /= MASperPIXEL
#print(shift_value*MASperPIXEL)
SHIFT_MASK[i, k] = -shift_pix_FPM
SHIFT_PSF[i, k] = shift_value
I_temp[k,:,:] = I_temp[k,:,:]*np.flipud(MASK_temp)
I_temp[k,:,:] = shift(I_temp[k,:,:], [shift_value,0], output=None, order=1, mode='constant', cval=0.0, prefilter=True)
I_temp[k] = I_temp[k, :, :]*np.flipud(MASK_temp)
I_temp[k] = shift(I_temp[k, :, :], [shift_value, 0], output=None, order=1, mode='constant', cval=0.0, prefilter=True)
image_filename = path_directory + 'PSF_HALO_ON_masked_centered_Nexp{0:04d}.fits'.format(i)
image_filename_nomask = path_directory + 'PSF_HALO_ON_notmasked_centered_Nexp{0:04d}.fits'.format(i)
fits.writeto(image_filename, I_temp)
fits.setval(image_filename, 'ORIGIN', value='ESO-ARMAZONES')
fits.setval(image_filename, 'TELESCOP', value='ESO-ELT')
fits.setval(image_filename, 'INSTRUME', value='HARMONI')
fits.setval(image_filename, 'EXPTIME', value='{:.3f}'.format(T_exp*3600))
fits.setval(image_filename, 'AIRMASS', value='{:.3f}'.format(1/np.cos(np.pi/180*(90-ELEVATION(DEC*1.0,HA_0)))))
fits.setval(image_filename, 'AIRMASS', value='{:.3f}'.format(1/np.cos(np.pi/180*(90-ELEVATION(DEC*1.0, HA_0)))))
fits.setval(image_filename, 'PI-COI', value='NIRANJAN')
fits.setval(image_filename, 'DISPELEM', value=BAND)
fits.setval(image_filename, 'APODIZER', value=APODIZER)
fits.setval(image_filename, 'FPM', value=FPM)
image_filename_nomask = path_directory + 'PSF_HALO_ON_notmasked_centered_Nexp{0:04d}.fits'.format(i)
fits.writeto(image_filename_nomask, I_temp_nomask)
fits.setval(image_filename_nomask, 'ORIGIN', value='ESO-ARMAZONES')
fits.setval(image_filename_nomask, 'TELESCOP', value='ESO-ELT')
fits.setval(image_filename_nomask, 'INSTRUME', value='HARMONI')
fits.setval(image_filename_nomask, 'EXPTIME', value='{:.3f}'.format(T_exp*3600))
fits.setval(image_filename_nomask, 'AIRMASS', value='{:.3f}'.format(1/np.cos(np.pi/180*(90-ELEVATION(DEC*1.0,HA_0)))))
fits.setval(image_filename_nomask, 'AIRMASS', value='{:.3f}'.format(1/np.cos(np.pi/180*(90-ELEVATION(DEC*1.0, HA_0)))))
fits.setval(image_filename_nomask, 'PI-COI', value='NIRANJAN')
fits.setval(image_filename_nomask, 'DISPELEM', value=BAND)
fits.setval(image_filename_nomask, 'APODIZER', value=APODIZER)
fits.setval(image_filename_nomask, 'FPM', value=FPM)
SR_filename = path_directory + 'SR.fits'.format(i)
shift_filename = path_directory + 'SHIFT_PSF.fits'
fits.writeto(shift_filename, SHIFT_PSF)
shift_filename = path_directory + 'SHIFT_MASK.fits'
fits.writeto(shift_filename, SHIFT_MASK)
SR_filename = path_directory + 'SR.fits'
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')
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment