Commit 88885e86 authored by Alexis Carlotti's avatar Alexis Carlotti
Browse files

Merge branch 'avigan' into 'master'

Merge branch avigan into master

See merge request !6
parents b2e4a9f6 7a661099
......@@ -830,7 +830,7 @@ for i in range(len(HA_vect)):
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)
image_filename_nomask = path_directory + 'PSF_HALO_ON_notmasked_notcentered_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')
......
import numpy as np
import logging
from scipy.ndimage import rotate
from scipy.ndimage import shift
from scipy.interpolate import interp2d
import multiprocessing as mp
from numpy import pi
import os
import sys
from time import localtime, strftime
import time
import matplotlib.pyplot as plt
from astropy.io import fits
from utils_basic import minterp2d, DISPER, VECT, VECT1, BEAMSHIFT, ELEVATION, PARANG, ft, MAKE_ELT, PSF, LD2MAS, MAS2LD, M4Magic, DSP, EELTNLS, APPLY_CROSSTALK, DISPER_ADC
def compute_image_at_wavelength(i, k):
#st = time.time()
WF_CP_0, WF_NCPA_0 = BEAMSHIFT(N, SCREEN, Pupil_shift_vect_0, D_opt_vect, Z_opt_vect, lambda_wfs, lambda_vect[k], 90-ROT_ANG_vect[i], ROT_ANG_vect[0], rot_flag*rot_vect, com_vect, gamma)
#log.info(time.time()-st)
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:.1f} nm: {2:.1f} 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))
#DSP_V, DSP_I = DSP(WF_0, 1132/1024, PUPILLE_0, 2, 30)
#log.info('HA={0}: Aberration level (2-30 lambda/D) at {1} nm: {2} nm'.format(HA_0, lambda_vect[k]*1e9, 1e9*DSP_V))
OWA = MAS2LD(FOV, 38.542, lambda_vect[k])/2
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])
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])
#log.info('Differential dispersion angle with 2.45um: {0}'.format(alpha_disp_0))
x, X, Y, R, T = VECT(len(SP), 1132/1024)
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
#HALO
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)
PH_STR = 2*(FT_PST[int(size_dsp//2+1), int(size_dsp//2+1)]-FT_PST)
# Optical transfer function of the atmosphere (long exposure model)
FTO_atm = np.exp(-PH_STR/2.)
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_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))
log.info('Wavelength {:4d}/{:4d}'.format(k+1, N_LD))
log.info(' => Strehl at {:4.2f}nm: {:4.2f}'.format(lambda_vect[k]*1e9, SR))
#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
#photons collected by the aperture. Careful: this normalization
#takes into account the central obscuration, and the apodizer
#throughput.
I_s_post_0 = I_s_post_0/PUP_NORM_INTENSITY*(2*OWA/N_im)**2
#[u, U, V, Ru, Tu]=VECT(N_im, 2*OWA);
#REG_100=(Ru<MAS2LD(125, 38.542, lambda_vect(k)*1e9)).*(Ru>MAS2LD(75, 38.542, lambda_vect(k)*1e9));
#REG_200=(Ru<MAS2LD(225, 38.542, lambda_vect(k)*1e9)).*(Ru>MAS2LD(175, 38.542, lambda_vect(k)*1e9));
#I_stat_nohalo=I_0/max(I_0(:));
#I_stat=I_s_post_0/max(I_s_post_0(:));
#disp(['5 sigma static Contrast at 100 mas - NO HALO: ' num2str(5*std(I_stat_nohalo(REG_100==1)))])
#disp(['5 sigma static Contrast at 200 mas - NO HALO: ' num2str(5*std(I_stat_nohalo(REG_200==1)))])
#disp(['Static Contrast -- std -- w/ E2E model: ' num2str(5*std(I_stat(REG==1)))])
#fits.writeto(path_directory + 'PSF_HALO_ON_Nexp{0:03d}_Nlbd{1:04d}.fits'.format(i, k), I_s_post_0)
# I1D = I_s_post_0[:, 107]
# I1D = I1D/sum(I1D)
# V1D = Y_mask[:, 107]
# 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
# if exhaustive_writing_flag==1:
# fits.writeto(path_directory + 'PSF_HALO_OFF_Nexp{0}_Nlbd{1}.fits'.format(i, k), I_0)
# WF_small_0_f = interp2d(x, x, WF_0, kind='linear')
# WF_small_0 = WF_small_0_f(x2, x2)
# WF_small_0 = WF_small_0*P_small_0
# fits.writeto(path_directory + 'WF_Nexp{0}_Nlbd{1}.fits'.format(i, k), WF_small_0)
return k, I_s_post_0, SR
if __name__ == '__main__':
# timing
start_time = time.time()
# find out how process was started
manual = True
if sys.argv[0] == 'MAIN_mp.py':
manual = False
HARMONI_DIR = '/Users/avigan/Work/HARMONI/Simulations/harmoni-hc/hrm-hc/'
COM_DIR = '/Users/avigan/Work/HARMONI/Simulations/harmoni-hc/hrm-hc/COMMON_FILES/'
# HARMONI_DIR = '/Users/carlotal/harmoni-hc/hrm-hc'
# 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;
#Artifical rotation suppressor (sky rotates, optics do not) ; 0 means no
#rotation ; 1 means rotation ; partial rotation is not possible
rot_flag = 1
#Flag for missing segments (1 = missing segments ; 0 = no missing segments)
miss_flag = 1
# Configuration of M1 amplitude errors: 1 & 2 refer to two different missing segment configurations
# A & B refer to two different reflectivity errors configurations
# In both cases there is a uniform distribution of reflectivity for each segment that ranges between 100% and 95%.
# This 95% value comes from paper "Towards a Stray Light Analysis of the ELT" from R. Holzlöhner at ESO
# Reflectivity errors are randomly distributed over M1, i.e., I have not assumed any specific recoating strategy (spiral, etc.)
ELT_amplitude_errors = 'CONFIG_1A'
#Flag for E-ELT M1 phase aberrations
eeltphase_flag = 1
#Flag for cross-talk and diffusion ; FWHM=1pix across x ; 2pix across lambda ; 0.5% diffusion across lambda
crosstalk_flag = 0
#Flag for writing HALO_off fits & WF fits
exhaustive_writing_flag = 0
# Fixed Dispersion Correction (optimized for 5-50 deg zenith range)
# ONLY used when forming the image since the ADC is not seen by (most) of HARMONI optics
ADC_flag = 1
#Zenith angle range [deg] !!! imposed by instrument !!! DO NOT CHANGE
z_min = 5
z_max = 50
#Seeing conditions - options: JQ1, JQ2, JQ3, MED
if manual:
seeing = 'MED'
else:
seeing = sys.argv[1]
#Sensing wavelength, ZELDA [m]
lambda_wfs = 1.175e-6
#CONFIG='Ideal';
#CONFIG='Relay_PDR_ZELDA_REQUIREMENT';
#CONFIG='Relay_PDR_ZELDA_GOAL';
CONFIG = 'Relay_PDR_ZELDA_OLD' #KEEP USING THIS ONE FOR NOW
RMS_gamma = 1 #%0.6;%0.3;
# Gamma factor to increase the strengh of the AO residuals (for DSP_gamma >1)
DSP_gamma = 1
if manual:
APODIZER = 'SP_1'
else:
APODIZER = sys.argv[2]
#disp(['Observation time will be: ' num2str(Delta_HA) 'h'])
#disp(['Observation starts ' num2str(HA_offset) 'h before meridian'])
#Declination of host star [degree]
DEC = -15
# -20 -> 5° zenith angle
# -10 -> 15°
# 0 -> 25°
# 10 -> 35°
# 20 -> 45°
#Mean Hour angle [hours]
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
#Exposure time of individual exposures [hours]
T_exp = 1/60
#Number of exposures
N_exp = np.int(T_obs/T_exp)# %60 exposures per hour = 1 per minute
##############################################################################
if manual:
BAND = 'TEST_H' # %Choose between the bands defined hereafter (use a valid string)
else:
BAND = sys.argv[3]
### Bands ###
VIS_band = [0.458, 0.639, 0.820, 3100]
IZJ_band = [0.811, 1.090, 1.369, 3355]
HK_band = [1.450, 1.950, 2.450, 3355]
# Test band with R=7000 resolution, but only 43 wavelengths
H_small_band = [1.620, 1.625, 1.630, 7000]
IZ_band = [0.830, 0.940, 1.050, 7104]
J_band = [1.046, 1.185, 1.324, 7104]
H_band = [1.435, 1.625, 1.815, 7104]
K_band = [1.951, 2.210, 2.469, 7104]
Z_high = [0.828, 0.865, 0.902, 17385]
H_high = [1.538, 1.608, 1.678, 17385]
K1_high = [2.017, 2.109, 2.201, 17385]
K2_high = [2.199, 2.300, 2.400, 17385]
if BAND == 'TEST_H':
N_LD = 3
LBD_min = H_band[0]*1e-6
LBD_max = H_band[2]*1e-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
LBD_max = K_band[2]*1e-6
if APODIZER == 'SP_1':
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_NoADC.fits')
FPM = 'FPM_K'
else:
# MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_ADC.fits')
# MASK_OFFSET = 0 # mas
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 7 # mas
FPM = 'FPM_H'
# MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_ADC.fits')
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_NoADC.fits')
FPM = 'FPM_K'
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = -2#0 # mas
FPM = 'FPM_K'
elif BAND == 'HK':
N_LD = round((HK_band[2]-HK_band[0])/(HK_band[1]/HK_band[3]))
LBD_min = HK_band[0]*1e-6
LBD_max = HK_band[2]*1e-6
if APODIZER == 'SP_1':
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_HK_SP1_NoADC.fits')
FPM = 'FPM_HK'
else:
# MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_ADC.fits')
# MASK_OFFSET = -22 # mas
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = -15 # mas
FPM = 'FPM_H'
# MASK = fits.getdata(COM_DIR + 'FPM_HK_SP1_ADC.fits')
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_HK_SP2_NoADC.fits')
FPM = 'FPM_HK'
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = -30 # mas
FPM = 'FPM_K'
# MASK = fits.getdata(COM_DIR + 'FPM_HK_SP2_ADC.fits')
elif BAND == 'H':
N_LD = round((H_band[2]-H_band[0])/(H_band[1]/H_band[3]))
LBD_min = H_band[0]*1e-6
LBD_max = H_band[2]*1e-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_H_SP2_ADC.fits')
MASK_OFFSET = 0 # mas
FPM = 'FPM_H'
elif BAND == 'K':
N_LD = round((K_band[2]-K_band[0])/(K_band[1]/K_band[3]))
LBD_min = K_band[0]*1e-6
LBD_max = K_band[2]*1e-6
if APODIZER == 'SP_1':
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_NoADC.fits')
FPM = 'FPM_K'
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 7 # mas
FPM = 'FPM_H'
# MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_ADC.fits')
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_NoADC.fits')
FPM = 'FPM_K'
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = 0 # mas
FPM = 'FPM_K'
elif BAND == 'H_high':
N_LD = round((H_high[2]-H_high[0])/(H_high[1]/H_high[3]))
LBD_min = H_high[0]*1e-6
LBD_max = H_high[2]*1e-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_H_SP2_ADC.fits')
MASK_OFFSET = 0 # mas
FPM = 'FPM_H'
#elif BAND == 'K_high':
# N_LD = round((K_high[2]-K_high[0])/(K_high[1]/K_high[3]))
# LBD_min = K_high[0]*1e-6
# LBD_max = K_high[2]*1e-6
elif BAND == 'K1_high':
N_LD = round((K1_high[2]-K1_high[0])/(K1_high[1]/K1_high[3]))
LBD_min = K1_high[0]*1e-6
LBD_max = K1_high[2]*1e-6
if APODIZER == 'SP_1':
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_NoADC.fits')
FPM = 'FPM_K'
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 7 # mas
FPM = 'FPM_H'
# MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_ADC.fits')
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_NoADC.fits')
FPM = 'FPM_K'
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = 0 # mas
FPM = 'FPM_K'
elif BAND == 'K2_high':
N_LD = round((K2_high[2]-K2_high[0])/(K2_high[1]/K2_high[3]))
LBD_min = K2_high[0]*1e-6
LBD_max = K2_high[2]*1e-6
if APODIZER == 'SP_1':
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_NoADC.fits')
FPM = 'FPM_K'
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 7 # mas
FPM = 'FPM_H'
# MASK = fits.getdata(COM_DIR + 'FPM_K_SP1_ADC.fits')
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_NoADC.fits')
FPM = 'FPM_K'
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = 0 # mas
FPM = 'FPM_K'
else:
#log.info('Invalid band ; switching to test band: 3 wavelengths (1.650, 1.925, 2.200 um)')
N_LD = 3
LBD_min = 1.65e-6
LBD_max = 2.20e-6
if APODIZER == 'SP_1':
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_HK_SP1_NoADC.fits')
FPM = 'FPM_HK'
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
# MASK = fits.getdata(COM_DIR + 'FPM_HK_SP1_ADC.fits')
MASK_OFFSET = -15 # mas
FPM = 'FPM_H'
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_HK_SP2_NoADC.fits')
FPM = 'FPM_HK'
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = -30 # mas
FPM = 'FPM_K'
# MASK = fits.getdata(COM_DIR + 'FPM_HK_SP2_ADC.fits')
# Science wavelength vector
if N_LD == 1:
lambda_vect = np.array([np.mean((LBD_min, LBD_max))])
else:
lambda_vect = np.linspace(LBD_min, LBD_max, N_LD)
#On sky field rotation necessary for ADI ; NOT A PLANET SIMULATION
Rho_planet = 100
Sky_Rot = np.arctan(LD2MAS(1, 38.542, 2400e-9)/Rho_planet)*180/pi
#Vector of hour angles
HA_vect = np.linspace(Mid_HA-T_obs/2, Mid_HA+T_obs/2, N_exp)
#Parallactic angles
PA_vect = PARANG(DEC, HA_vect)
#Time computed to give an idea of the minimum time required for ADI
# INDICATIVE ONLY ; DOES NOT HAVE AN IMPACT ON THE CODE
ADI_time_FIRST = 60*Sky_Rot/np.abs((PARANG(DEC, np.min(HA_vect)-1/60)-PARANG(DEC, np.min(HA_vect)+1/60))/(2/60))
ADI_time_MID = 60*Sky_Rot/np.abs((PARANG(DEC, Mid_HA-1/60)-PARANG(DEC, Mid_HA+1/60))/(2/60))
ADI_time_LAST = 60*Sky_Rot/np.abs((PARANG(DEC, np.max(HA_vect)-1/60)-PARANG(DEC, np.max(HA_vect)+1/60))/(2/60))
# Pupil shift model
Pupil_shift = 1.6/100
Pupil_apodizer_differential_shift = 0.2/100
Pupil_shift_period = 1# %period of pupil shift in hours
Pupil_shift_vect = Pupil_shift*np.cos(2*pi*HA_vect/Pupil_shift_period)# %simple model for pupil shift ; not sure at this time how it will actually move
Pupil_shift_diff_vect = Pupil_apodizer_differential_shift*np.cos(2*pi*HA_vect/Pupil_shift_period)# %simple model for pupil shift ; not sure at this time how it will actually move
#Simulating the evolution of the wind direction as the rotation of the wind with the paralactic angle
Wind_Angle_vect = PARANG(DEC, HA_vect)
# _ __ _ _ _ _
# | | / _| | | | | | | (_)
# ___ _ __ __| | ___ | |_ _ __ __ _ _ __ __ _ _ __ ___ ___| |_ ___ _ __ ___ ___| | ___ ___| |_ _ ___ _ __
# / _ \ '_ \ / _` | / _ \| _| | '_ \ / _` | '__/ _` | '_ ` _ \ / _ \ __/ _ \ '__| / __|/ _ \ |/ _ \/ __| __| |/ _ \| '_ \
# | __/ | | | (_| | | (_) | | | |_) | (_| | | | (_| | | | | | | __/ || __/ | \__ \ __/ | __/ (__| |_| | (_) | | | |
# \___|_| |_|\__,_| \___/|_| | .__/ \__,_|_| \__, _|_||_| |_|\___|\__\___|_| |___/\___|_|\___|\___|\__|_|\___/|_| |_|
# | |
# |_|
path_directory = HARMONI_DIR + '/HC_seeing={}_apo={}_band={}_{}/'.format(seeing, APODIZER, BAND, strftime("%d-%b-%Y-%H-%M-%S", localtime()))
if (os.path.isdir(path_directory) is False):
os.mkdir(path_directory)
else:
print('Old directory exists!')
logging.basicConfig(filename=path_directory + 'simulation_log.txt', filemode='w',
format='%(message)s', level=logging.DEBUG)
log = logging.getLogger('harmoni-hc')
#if not log.handlers:
handler = logging.StreamHandler(sys.stdout)
log.addHandler(handler)
log.info('DATA MAKER')
log.info('Number of exposures: {0}'.format(N_exp))
log.info('Total observing time: {0} min'.format(T_obs*60))
log.info('Time between two exposures: {0} min'.format(T_obs/N_exp*60))
if (gamma == 1):
log.info('Dispersion: YES')
elif (gamma == 0):
log.info('Dispersion: NO')
else:
log.info('Dispersion: {0}'.format(gamma))
if (rot_flag == 1):
log.info('Rotation: YES')
else:
log.info('Rotation: NO')
if (miss_flag == 1):
log.info('M1 Amplitude errors: YES')
else:
log.info('M1 Amplitude errors: NO')
if (eeltphase_flag == 1):
log.info('M1 phase errors: YES')
else:
log.info('M1 phase errors: NO')
if (ADC_flag == 1):
log.info('Fixed ADC: YES')
else:
log.info('Fixed ADC: NO')
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')
N = len(EELT_phase)
PUP = 1-1*(EELT_phase == 0)
if (APODIZER == 'SP_1'):
IWA_SP = 5
A = fits.getdata(COM_DIR + 'HSP1.fits')
elif (APODIZER == 'SP_2'):
IWA_SP = 7
A = fits.getdata(COM_DIR + 'HSP2.fits')
else:
log.info('Invalid apodizer ; switching to SP_1')
A = fits.getdata('HSP1.fits')
log.info('Apodizer: ' + APODIZER)
log.info('Band: ' + BAND)
log.info('FPM: ' + FPM) # this parameter is set as a function of the apodizer & band choice
#EA = ft(A, 1, 1132/1024, 1132/1024, 566/2, 566/2, 1, 566, 566)
#A_big = -np.real(ft(EA, 1, 566/2, 566/2, 1132/1024, 1132/1024, 1, 1132, 1132))
#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)
MASK[np.abs(Y_mask) > 610/2] = 0
fits.writeto(path_directory + 'FOCAL_PLANE_MASK.fits', MASK)
start_dir = os.getcwd()
#os.chdir(path_directory)
Ideal_PSF = PSF(ft(SP, 1, 1132/1024, 1132/1024, MAS2LD(FOV, 38.542, 1650e-9), MAS2LD(FOV, 38.542, 1650e-9), 1, N_im, N_im))
fits.writeto(path_directory + 'Ideal_PSF.fits', Ideal_PSF, overwrite=True)
#os.chdir(start_dir)
respix = 0.0376387
diam_allglass = 36.905
diam_ext = 39.146
diam_obsc = 11.056
diam_obsc_allglass = 11.213
cobsc_allglass = diam_obsc_allglass / diam_allglass
entrance_pup = 38.542
D = 38.542 #37 #why not 38.542 ??? check difference
BigD = 38.542*1132/1024
M1pupil1132_1024 = fits.getdata(COM_DIR + 'PUPIL_1132_1024D38.542.fits')
PUP_NORM_INTENSITY = np.sum(M1pupil1132_1024**2*(1132/1024/1132)*(1132/1024/1132))
if ELT_amplitude_errors == 'CONFIG_1A':
MISSSEG = fits.getdata(COM_DIR + 'ELT_1a_MS=7_RE=5.fits')
elif ELT_amplitude_errors == 'CONFIG_1B':
MISSSEG = fits.getdata(COM_DIR + 'ELT_1b_MS=7_RE=5.fits')
elif ELT_amplitude_errors == 'CONFIG_2A':
MISSSEG = fits.getdata(COM_DIR + 'ELT_2a_MS=7_RE=5.fits')
elif ELT_amplitude_errors == 'CONFIG_2B':
MISSSEG = fits.getdata(COM_DIR + 'ELT_2b_MS=7_RE=5.fits')