Commit bb6415e7 authored by Arthur Vigan's avatar Arthur Vigan
Browse files

Setup for use in cluster

parent 1ec96875
import numpy as np
import logging
from scipy.ndimage import rotate
from scipy.ndimage import shift
from scipy.interpolate import interp2d
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
if __name__ == '__main__':
# timing
start_time = time.time()
# find out how process was started
manual = True
if sys.argv[0] == 'MAIN_multiproc.py':
manual = False
# paths
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
#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 = 'JQ1'
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.463, 0.641, 0.820, 3500]
IZJ_band = [0.810, 1.090, 1.369, 3500]
HK_band = [1.45, 1.950, 2.45, 3500]
# 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, 7000]
J_band = [1.046, 1.185, 1.324, 7000]
H_band = [1.435, 1.625, 1.815, 7000]
K_band = [1.951, 2.210, 2.469, 7000]
Z_high = [0.829, 0.866, 0.902, 17000]
#J_high=[1.190, 1.245, 1.300, 17000] # not actually present in harmoni
H_high = [1.541, 1.609, 1.676, 17000]
#K_high=[2.093, 2.190, 2.287, 17000] # old definition ; not correct anymore
K1_high = [2.024, 2.113, 2.202, 17000]
K2_high = [2.197, 2.294, 2.390, 17000]
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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP1_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_1'
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_NoADC.fits')
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_2'
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')
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_2'
#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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_3'
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')
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_2'
#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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = -30 #mas
FPM = 'FPM_3'
#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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP1_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_1'
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_NoADC.fits')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_2'
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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 7 #mas
FPM = 'FPM_2'
#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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_3'
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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP1_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_1'
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_NoADC.fits')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_2'
#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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 7 #mas
FPM = 'FPM_2'
#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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_3'
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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_H_SP2_ADC.fits')
MASK_OFFSET = 7 #mas
FPM = 'FPM_2'
#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')
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = 0 #mas
FPM = 'FPM_3'
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')
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_2'
else:
if ADC_flag == 0:
MASK = fits.getdata(COM_DIR + 'FPM_HK_SP2_NoADC.fits')
else:
MASK = fits.getdata(COM_DIR + 'FPM_K_SP2_ADC.fits')
MASK_OFFSET = -30 #mas
FPM = 'FPM_3'
#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!')
stop
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('Missing Segments: YES')
else:
log.info('Missing Segments: NO')
if (eeltphase_flag == 1):
log.info('EELT phase: YES')
else:
log.info('EELT phase: 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[MASK == 0] = 1e-4
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))
MISSSEG = fits.getdata(COM_DIR + 'MissingSegments_1132.fits')
MISSSEG = miss_flag*MISSSEG
ORIGINAL_PUPIL = M1pupil1132_1024
M1pupil1132_1024 = M1pupil1132_1024*(1-MISSSEG)
P_M4 = EELTNLS(1024, 100, 100, 110) ### EELTNLS not written yet!!! ###
P_M4 = np.pad(P_M4, [54, 54], 'constant', constant_values=(0, 0))
#P_M4, P_M4_col = MAKE_ELT(COM_DIR, 1132, 1132/1024*38.542, 0, 0)
x, X, Y, R, T = VECT(N, respix*N)
P = M1pupil1132_1024
log.info('Star DEC: {:05.2f} deg'.format(DEC))
log.info('Star first HA: {:05.2f} h ; last HA: {:05.2f} h'.format(np.min(HA_vect), np.max(HA_vect)))
log.info('Star initial Zenith angle: {:05.2f} deg ; final Zenith angle: {:05.2f} deg'.format(90-ELEVATION(DEC, np.min(HA_vect)), 90-ELEVATION(DEC, np.max(HA_vect))))
ROT_ANG_vect = ELEVATION(DEC, HA_vect)
ROT_ANG_diff = np.mean(abs(np.diff(ROT_ANG_vect)))
log.info('Mean rotation angle between two exposures: {:05.2f} deg'.format(ROT_ANG_diff))
N_act = 80
t = 0.5
pix_per_LD = 2
M = 2*np.ceil(pix_per_LD*MAS2LD(600, BigD, lambda_vect[0]))
OWA_W = 64
D_tel=D
if (CONFIG == 'Ideal'):
RMS_opt_vect = 10e-9#rms_main[0]
D_opt_vect = D #;%[D, D, D, D, D, D, D, D];
Z_opt_vect = 218e3 #;%[218e3];%[500e3, 500e3, 500e3, 500e3, 500e3, 500e3, 500e3, 500e3];
rot_vect = 1 #;%[0 0 0 1 0 0 0 0];
com_vect = 1 #;%[1 1 1 1 0 0 0 0];
elif (CONFIG == 'Relay_PDR_ZELDA_REQUIREMENT'):
RMS_opt_vect = RMS_gamma*[82.1, 20.5, 30.0, 76.9, 615.4, 95.7, 49.2, 13.7, 12.1, 12.7, 13.7, 14.8, 6.1, 6.1, 6.1, 6.1]
D_opt_vect = [D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D]
Z_opt_vect = [120.9e3, 478e3, 323e3, 129.5e3, 12.3e3, 104.9e3, 197.2e3, 713e3, 814e3, 768e3, 716e3, 656e3, 2.1e3, 744e3, 892e3, 2460e3]
rot_vect = [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
com_vect = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
elif (CONFIG == 'Relay_PDR_ZELDA_GOAL'):
RMS_opt_vect = RMS_gamma*[43.6, 10.9, 15.9, 40.9, 326.9, 50.3, 26.2, 7.3, 6.4, 6.7, 7.3, 7.9, 2.9, 2.9, 2.9, 2.9]
D_opt_vect = [D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D]
Z_opt_vect = [120.9e3, 478e3, 323e3, 129.5e3, 12.3e3, 104.9e3, 197.2e3, 713e3, 814e3, 768e3, 716e3, 656e3, 2.1e3, 744e3, 892e3, 2460e3]
rot_vect = [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
com_vect = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
elif (CONFIG == 'Relay_PDR_ZELDA_OLD'): # %Last optics reprensent SCAO dichroic, whose conj. height is unknown for the moment, but should be a bit higher than the last fold (505 km)
RMS_opt_vect = RMS_gamma*[36.6, 11.49, 15.91, 35, 59.72, 39.88, 25.31, 8.1, 7.21, 7.57, 8.06, 8.72, 5, 11, 5, 5]
D_opt_vect = [D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D]
Z_opt_vect = [120.9e3, 478e3, 323e3, 129.5e3, 12.3e3, 104.9e3, 197.2e3, 713e3, 814e3, 768e3, 716e3, 656e3, 2.1e3, 744e3, 892e3, 2460e3]
rot_vect = [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
com_vect = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
elif (CONFIG == 'Relay_PDR_BASIC_REQUIREMENT'):
RMS_opt_vect = RMS_gamma*[25.1, 6.4, 9.5, 23.7, 246.2, 29.3, 15.6, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30]
D_opt_vect = [D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D]
Z_opt_vect = [120.9e3, 478e3, 323e3, 129.5e3, 12.3e3, 104.9e3, 197.2e3, 713e3, 2460e3, 1000e3, 1000e3, 1000e3, 1000e3, 1000e3, 1000e3, 1000e3, 1000e3]
rot_vect = [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
com_vect = [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
else:
RMS_opt_vect = RMS_gamma*[36.6, 11.49, 15.91, 35, 59.72, 39.88, 25.31, 8.1, 7.21, 7.57, 8.06, 8.72, 5, 11, 5, 5]
D_opt_vect = [D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D]
Z_opt_vect = [120.9e3, 478e3, 323e3, 129.5e3, 12.3e3, 104.9e3, 197.2e3, 713e3, 814e3, 768e3, 716e3, 656e3, 2.1e3, 744e3, 892e3, 2460e3]
rot_vect = [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
com_vect = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
N_optics = len(rot_vect)
#%the phase screen is initially 1.2 times larger than the actual phase
#%screen to allow a 10% shift
SCREEN = np.zeros([1358, 1358, N_optics])
log.info('Reading the optical surfaces...')
#disp('---------------------------------------------')
#disp(['Optics #' char(9) 'Conj. height (km)' char(9) 'Rotation?' char(9) 'Com path?'])
for k in range(N_optics):
#% READING THE SURFACE TEMPLATES TO MAKE SURE THAT COMPARISONS WILL NOT
#% DEPEND ON SURFACE CHOICES
#%SCREEN(:, :, k)=RMS_opt_vect(k)*load(['SCREEN_TEMPLATE_1nmRMS_Nbr' num2str(k) '.dat']);
SCREEN[:, :, k] = RMS_opt_vect[k]*fits.getdata(COM_DIR + 'SCREEN_NEW_N={0}.fits'.format(k+1))*1e-9
#disp('---------------------------------------------')
log.info('Optics loaded.')
#disp(D_opt_vect)
#disp(Z_opt_vect)
# save import parameters
fits.writeto(path_directory + 'PARANG_VECT.fits', PA_vect, overwrite=True)
fits.writeto(path_directory + 'HA_VECT.fits', HA_vect, overwrite=True)
fits.writeto(path_directory + 'WAVELENGTH_VECT.fits', lambda_vect, overwrite=True)
# _ _ _ ____
# | | | | /\ | | / __ \
# | |__| | / \ | | | | | |
# | __ | / /\ \ | | | | | |
# | | | |/ ____ \| |___| |__| |
# |_| |_/_/ \_\______\____/
#
dsp = fits.getdata(COM_DIR + seeing + '_Res-SpatialPSD_Scaled.fits')
dsp = dsp*DSP_gamma**2
size_dsp = len(dsp[0])
dsp_array= dsp
DiamFTO = 2#4
dim = 200#size_dsp/DiamFTO
fmax = (1024/1132)*size_dsp/2 # size_dsp/2
f1D, Xf, Yf, Rf, Tf = VECT(size_dsp, 2*fmax)
fudge_df = 1#(1024/1132)**2
df = np.abs(f1D[1]-f1D[0])*fudge_df
# _____ _____ __ __ _ _ _ _______ _____ _ _ _____ ________ _______ ____ _____ _ _ _____ ______ _____
# / ____|_ _| \/ | | | | | /\|__ __|_ _| \ | |/ ____| | ____\ \ / / __ \ / __ \ / ____| | | | __ \| ____|/ ____|
# | (___ | | | \ / | | | | | / \ | | | | | \| | | __ | |__ \ V /| |__) | | | | (___ | | | | |__) | |__ | (___
# \___ \ | | | |\/| | | | | | / /\ \ | | | | | . ` | | |_ | | __| > < | ___/| | | |\___ \| | | | _ /| __| \___ \
# ____) |_| |_| | | | |__| | |____ / ____ \| | _| |_| |\ | |__| | | |____ / . \| | | |__| |____) | |__| | | \ \| |____ ____) |
# |_____/|_____|_| |_|\____/|______/_/ \_\_| |_____|_| \_|\_____| |______/_/ \_\_| \____/|_____/ \____/|_| \_\______|_____/
#
SR = np.zeros((len(HA_vect), N_LD))