diff --git a/ADI_processing_old.py b/ADI_processing_old.py new file mode 100644 index 0000000..76a64b6 --- /dev/null +++ b/ADI_processing_old.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python3 + +import vip_hci +import numpy as np +from astropy.io import fits +import matplotlib.pyplot as plt + + +out_dir = str('./output_files/') + +#### inputs psfs and cube +psf_noMask = vip_hci.fits.open_fits(out_dir + 'test_PSF_ELT.fits') # non-coronagraphic PSF for normalization (just a ELT PSF) +psf = vip_hci.fits.open_fits(out_dir +'test_PSF_OFFAXIS.fits') # off-axis PSF (non-normalized) (No mask--with apodizer and LS) +#cube = vip_hci.fits.open_fits(out_dir + 'Heidelberg_100_PSF_cube_VC.fits') # coronagraphic PSFs cube (non-normalized) +#cube = vip_hci.fits.open_fits(out_dir + 'Heidelberg_1000_PSF_cube_VC.fits') # coronagraphic PSFs cube (non-normalized) +#cube = vip_hci.fits.open_fits(out_dir + 'Linz_100_PSF_cube_VC.fits') # coronagraphic PSFs cube (non-normalized) +cube = vip_hci.fits.open_fits(out_dir + 'test_PSF_cube_VC.fits')[0:10] # coronagraphic PSFs cube (non-normalized) + + +length,x,y = cube.shape + +psf= psf/np.sum(psf_noMask) # off-axis PSF normalized wrt the total flux of the non-coronagraphic PSF +cube = cube/np.sum(psf_noMask) # cube normalized wrt the total flux of the non-coronagraphic PSF + +angs = np.linspace(-160,160,length) + +exp_t = 60*60/length # exposure time in sec for 1 hr sequence +Lmag = 5 # choose the star magnitude here +star_ref = exp_t * 3.4539e9 # flux in ADU corresponding to a 5th mag star in the NACO-Lprime filter (METIS with no coronagraph) +star_flux = star_ref * 10**(-0.4*(Lmag-5)) # actual star flux for the chosen magnitude +bckg = exp_t * 229.78e3 # flux in ADU/pix corresponding to the L-band background in the NACO-Lprime filter + +#### pixelscale (arcsec/pix) +psc_simu = 0.005 # plate scale (arcsec/pix) in the simulations +psc_metis = 0.00525 # plate scale (arcsec/pix) of the METIS L-band detector + + + +# to calculate the transmission +trans = 0.34 # if VC then use '0.81' and '0.34' in case of RAVC + +# Resample the disk PSF and cube with the pixelscale of the Metis instrument +cube_metis = vip_hci.preproc.cube_px_resampling(cube, psc_simu/psc_metis) +cube_metis = cube_metis.clip(min=0) # to avoid negative values due to the resampling +#cube_metis[cube_metis < 0] = 0 # to avoid negative values due to the resampling +psf_metis = vip_hci.preproc.frame_px_resampling(psf, psc_simu/psc_metis) + + +xpsf = int((psf_metis.shape[0])/2) +ypsf = int((psf_metis.shape[0])/2) + +w = 19 +fit = vip_hci.var.fit_2dgaussian(psf_metis[xpsf-w:xpsf+w+1, ypsf-w:ypsf+w+1], True, (w,w), debug=False, full_output=True) + +# cropping the off-axis PSF +fwhmx = fit['fwhm_x'] +fwhmy = fit['fwhm_y'] +xcen = fit['centroid_x'] +ycen = fit['centroid_y'] + +psf_cen = vip_hci.preproc.frame_shift(psf_metis, w-xcen, w-ycen) +psf_crop = psf_cen[xpsf-w:xpsf+w+1, ypsf-w:ypsf+w+1] + +# Derive the FWHM +fwhm = np.mean([fwhmx,fwhmy]) + +# Adding the stellar flux and the background +psf_phot = psf_crop * star_flux +cube_phot = cube_metis * star_flux +cube_phot = cube_phot + bckg*trans # need to multiply the backgound fo the transmission + +ran = np.random.randn(cube_phot.shape[0],cube_phot.shape[1],cube_phot.shape[2]) +cube_phot = cube_phot + ran * np.sqrt(cube_phot) + +starphot = vip_hci.metrics.aperture_flux(psf_phot,[w],[w],fwhm) + +# ADI, ADI-PCA +cube_out_adi_mean100ms, cube_der_adi_mean100ms, image_adi_mean100ms = vip_hci.medsub.median_sub(cube_phot, angs, full_output=True) # ADI + +ravc2_mean100ms_cc_adi = vip_hci.metrics.contrast_curve(cube_phot, angs, psf_crop, fwhm, psc_metis, starphot[0], vip_hci.medsub.median_sub, nbranch=1, debug=False, plot=False) # ADI contrast curves + +x = ravc2_mean100ms_cc_adi.get('distance_arcsec') +y1 = ravc2_mean100ms_cc_adi.get('sensitivity_gaussian') +y2 = ravc2_mean100ms_cc_adi.get('sensitivity_student') + +#np.savetxt('Heidelberg_1000_y1.txt', y1, delimiter=',') +#np.savetxt('Heidelberg_1000_x.txt', x, delimiter=',') + +#np.savetxt('Linz_1000_y1.txt', y1, delimiter=',') +#np.savetxt('Linz_1000_x.txt', x, delimiter=',') + +# +plt.figure(1, figsize=(12,8)) +ax = plt.axes() +ax.semilogy(x, y1) +ax.semilogy(x,y2) +ax.grid() +ax.legend() +ax.set_ylabel("Contrast") +ax.set_xlabel("Arcsec") + diff --git a/example_coronagraph_psf.py b/example_coronagraph_psf.py new file mode 100644 index 0000000..a3479f9 --- /dev/null +++ b/example_coronagraph_psf.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 + +# ============================================================================= +# Example file for creating a coronagraphic/non-coronagraphic METIS PSF +# ============================================================================= + +# Required libraries +import matplotlib.pyplot as plt # for plotting simulated PSFs +import numpy as np +from astropy.io import fits + +from heeps.config import download_from_gdrive +from heeps import metis_hci +from heeps.config import conf, pre_sim + +pre_sim(conf) + +""" +Default simulation cofiguration defined in "read_config.py" can be +overridden here by updating the dictionary +""" +conf['WAVELENGTH'] = 3.80*10**-6 +conf['CHARGE'] = 2 # charge is modified here +conf['MODE'] = 'VC' +conf['STATIC_NCPA'] = False + +# 2-D pahse screen +atm_screen = 0 + +# getting multi-cube phase screen +download_from_gdrive(conf['GDRIVE_ID'], conf['INPUT_DIR'], conf['ATM_SCREEN_CUBE']) +atm_screen = fits.getdata(conf['INPUT_DIR']+ conf['ATM_SCREEN_CUBE'])[0] + +TILT = np.array(conf['TILT_2D']) +#TILT = np.random.randn(conf['TILT_CUBE'],2) + +psf = metis_hci(conf, atm_screen, TILT) + + +plt.figure() +plt.imshow(psf**0.05) +plt.colorbar() + +plt.show() + diff --git a/heeps/__init__.py b/heeps/__init__.py new file mode 100644 index 0000000..166cb64 --- /dev/null +++ b/heeps/__init__.py @@ -0,0 +1,16 @@ +from __future__ import (absolute_import) + +__version__ = "0.1.0" + + +from . import abberations +from . import adi +from . import config +from . import coronagraphs +from . import detector +from . import fits +from . import pupil + +from .metis_hci import metis_hci, coronagraphs + + diff --git a/heeps/abberations/__init__.py b/heeps/abberations/__init__.py new file mode 100644 index 0000000..63caefb --- /dev/null +++ b/heeps/abberations/__init__.py @@ -0,0 +1,7 @@ +from __future__ import absolute_import + +from .island_effect_piston import island_effect_piston +from .atmosphere import atmosphere +from .wavefront_abberations import wavefront_abberations +from .static_ncpa import static_ncpa + diff --git a/heeps/abberations/atmosphere.py b/heeps/abberations/atmosphere.py new file mode 100644 index 0000000..f607072 --- /dev/null +++ b/heeps/abberations/atmosphere.py @@ -0,0 +1,34 @@ +import numpy as np +import cv2 +import proper +import math +from astropy.io import fits + +def atmosphere(wf, npupil, atm_screen, Debug_print, Debug): + + n = int(proper.prop_get_gridsize(wf)) + lamda=proper.prop_get_wavelength(wf) #save the wavelength value [m] into lamda + + screen = atm_screen # read the phase screen + if (Debug_print == True): + print ("screen", screen.shape) + screen_pixels = (screen.shape)[0] # size of the phase screen + scaling_factor_screen = float(npupil)/float(screen_pixels) # scaling factor the phase screen to the simulation pupil size + if (Debug_print == True): + print ("scaling_factor_screen: ", scaling_factor_screen) + screen_scale = cv2.resize(screen.astype(np.float32), (0,0), fx=scaling_factor_screen, fy=scaling_factor_screen, interpolation=cv2.INTER_LINEAR) # scale the the phase screen to the simulation pupil size + if (Debug_print == True): + print ("screen_resample", screen_scale.shape) + screen_large = np.zeros((n,n)) # define an array of n-0s, where to insert the screen + if (Debug_print == True): + print("n: ", n) + print("npupil: ", npupil) + screen_large[int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2),int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2)] =screen_scale # insert the scaled screen into the 0s grid + + + lambda2=lamda/(1e-6) # need to use lambda in microns + screen_atm = np.exp(1j*screen_large/lambda2*2*math.pi) + proper.prop_multiply(wf, screen_atm) # multiply the atm screen to teh wavefront + + + return wf diff --git a/heeps/abberations/island_effect_piston.py b/heeps/abberations/island_effect_piston.py new file mode 100644 index 0000000..dfffbb5 --- /dev/null +++ b/heeps/abberations/island_effect_piston.py @@ -0,0 +1,50 @@ +import numpy as np +import cv2 +import proper +import math +from astropy.io import fits +import os + +def island_effect_piston(wf, npupil, Island_Piston, path, Debug_print, Debug): + + n = int(proper.prop_get_gridsize(wf)) + lamda=proper.prop_get_wavelength(wf) #save the wavelength value [m] into lamda + + PACKAGE_PATH = os.path.abspath(os.path.join(__file__, os.pardir)) + petal1 = fits.getdata(PACKAGE_PATH+'/1024_pixelsize5mas_petal1_243px.fits') + petal2 = fits.getdata(PACKAGE_PATH+'/1024_pixelsize5mas_petal2_243px.fits') + petal3 = fits.getdata(PACKAGE_PATH+'/1024_pixelsize5mas_petal3_243px.fits') + petal4 = fits.getdata(PACKAGE_PATH+'/1024_pixelsize5mas_petal4_243px.fits') + petal5 = fits.getdata(PACKAGE_PATH+'/1024_pixelsize5mas_petal5_243px.fits') + petal6 = fits.getdata(PACKAGE_PATH+'/1024_pixelsize5mas_petal6_243px.fits') + + piston_petal1 = Island_Piston[0]*petal1 + piston_petal2 = Island_Piston[1]*petal2 + piston_petal3 = Island_Piston[2]*petal3 + piston_petal4 = Island_Piston[3]*petal4 + piston_petal5 = Island_Piston[4]*petal5 + piston_petal6 = Island_Piston[5]*petal6 + + piston = piston_petal1 + piston_petal2 + piston_petal3 + piston_petal4 + piston_petal5 + piston_petal6 + + piston_pixels = (piston.shape)[0]## fits file size + scaling_factor = float(npupil)/float(piston_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + if (Debug_print==True): + print ("scaling_factor: ", scaling_factor) + piston_scale = cv2.resize(piston.astype(np.float32), (0,0), fx=scaling_factor, fy=scaling_factor, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + if (Debug_print==True): + print ("piston_resample", piston_scale.shape) + piston_large = np.zeros((n,n)) # define an array of n-0s, where to insert the pupuil + if (Debug_print==True): + print("n: ", n) + print("npupil: ", npupil) + piston_large[int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2),int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2)] =piston_scale # insert the scaled pupil into the 0s grid + if (Debug == 1): + fits.writeto(path+'piston_phase.fits', piston_large, overwrite=True) # fits file for the screen + + + lambda2=lamda/(1e-6) # need to use lambda in microns + piston_phase = np.exp(1j*piston_large/lambda2*2*math.pi) + proper.prop_multiply(wf, piston_phase) # multiply the atm screen to teh wavefront + + return diff --git a/heeps/abberations/static_ncpa.py b/heeps/abberations/static_ncpa.py new file mode 100644 index 0000000..17d498c --- /dev/null +++ b/heeps/abberations/static_ncpa.py @@ -0,0 +1,18 @@ +import numpy as np +import cv2 +import proper + +def static_ncpa(wf, npupil, screen): + + n = int(proper.prop_get_gridsize(wf)) + + screen_pixels = (screen.shape)[0] # size of the phase screen + sf = float(npupil)/float(screen_pixels) # scaling factor the phase screen to the simulation pupil size + + screen_scale = cv2.resize(screen.astype(np.float32), (0,0), fx = sf, fy = sf, interpolation=cv2.INTER_LINEAR) # scale the the phase screen to the simulation pupil size + screen_large = np.zeros((n,n)) # define an array of n-0s, where to insert the screen + screen_large[int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2),int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2)] =screen_scale # insert the scaled screen into the 0s grid + + proper.prop_add_phase(wf, screen_large) + + return wf diff --git a/heeps/abberations/wavefront_abberations.py b/heeps/abberations/wavefront_abberations.py new file mode 100644 index 0000000..d33a631 --- /dev/null +++ b/heeps/abberations/wavefront_abberations.py @@ -0,0 +1,48 @@ +import numpy as np +import proper +from astropy.io import fits + +from .island_effect_piston import island_effect_piston +from .atmosphere import atmosphere +from .static_ncpa import static_ncpa + + +def wavefront_abberations(wfo, conf,atm_screen, TILT): + + diam = conf['DIAM'] + gridsize = conf['GRIDSIZE'] + pixelsize = conf['PIXEL_SCALE'] + wavelength = conf['WAVELENGTH'] + lamda = proper.prop_get_wavelength(wfo) + Debug = conf['DEBUG'] + Debug_print = conf['DEBUG_PRINT'] + + Island_Piston = np.array(conf['ISLAND_PISTON']) +# atm_screen = fits.getdata(input_dir + '/' + atm_screen) + beam_ratio = pixelsize*4.85e-9/(wavelength/diam) + npupil = np.ceil(gridsize*beam_ratio) # compute the pupil size --> has to be ODD (proper puts the center in the up right pixel next to the grid center) + + if npupil % 2 == 0: + npupil = npupil + 1 + + if (isinstance(atm_screen, (list, tuple, np.ndarray)) == True) and (atm_screen.ndim >= 2): # when the atmosphere is present + atmosphere(wfo, npupil, atm_screen, Debug_print, Debug) + + if (all(v == 0 for v in Island_Piston) == False): # when the piston is present + island_effect_piston(wfo, npupil, Island_Piston, Debug_print, Debug) + + if conf['STATIC_NCPA'] == True: + filename = conf['IMG_LM_SCAO'] + phase_screen = fits.getdata(conf['INPUT_DIR'] + filename) + phase_screen = np.nan_to_num(phase_screen) + phase_screen *= 10**-9 # scale the wavelenth to nm + static_ncpa(wfo, npupil, phase_screen) + + if (TILT.any != 0.): # when tip/tilt + if (Debug_print==True): + print("TILT: ", TILT) + print("lamda: ", lamda) + tiptilt = (np.multiply(TILT, lamda))/4 # translate the tip/tilt from lambda/D into RMS phase errors + proper.prop_zernikes(wfo, [2,3], tiptilt) # 2-->xtilt, 3-->ytilt + + return wfo diff --git a/heeps/adi/__init__.py b/heeps/adi/__init__.py new file mode 100644 index 0000000..53f82f5 --- /dev/null +++ b/heeps/adi/__init__.py @@ -0,0 +1,3 @@ +from __future__ import absolute_import + + diff --git a/heeps/config/__init__.py b/heeps/config/__init__.py new file mode 100644 index 0000000..9a33591 --- /dev/null +++ b/heeps/config/__init__.py @@ -0,0 +1,5 @@ +from __future__ import absolute_import + +from .download_from_gdrive import download_from_gdrive +from .read_config import * +from .pre_sim import pre_sim diff --git a/heeps/config/download_from_gdrive.py b/heeps/config/download_from_gdrive.py new file mode 100644 index 0000000..41b0619 --- /dev/null +++ b/heeps/config/download_from_gdrive.py @@ -0,0 +1,34 @@ +import requests +import os + +'''Downloads large mulit-cube phase screens from google drive''' + +def download_from_gdrive(id, destination, filename): + my_file = str(destination + '/'+ filename) + destination = destination + '/'+ filename + if (os.path.isfile(my_file) == False): + print('Downloading multi-cube phase screen from the google drive.....') + URL = "https://docs.google.com/uc?export=download" + session = requests.Session() + response = session.get(URL, params = { 'id' : id }, stream = True) + token = get_confirm_token(response) + if token: + params = { 'id' : id, 'confirm' : token } + response = session.get(URL, params = params, stream = True) + + save_response_content(response, destination) + +def get_confirm_token(response): + for key, value in response.cookies.items(): + if key.startswith('download_warning'): + return value + + return None + +def save_response_content(response, destination): + CHUNK_SIZE = 32768 + + with open(destination, "wb") as f: + for chunk in response.iter_content(CHUNK_SIZE): + if chunk: # filter out keep-alive new chunks + f.write(chunk) diff --git a/heeps/config/pre_sim.py b/heeps/config/pre_sim.py new file mode 100644 index 0000000..3f8ed31 --- /dev/null +++ b/heeps/config/pre_sim.py @@ -0,0 +1,40 @@ +import os +import proper +import zipfile + +from .download_from_gdrive import download_from_gdrive + +def pre_sim(conf): + # Creates required directories for simulation + os.makedirs(conf['OUT_DIR'] , exist_ok=True) + os.makedirs(conf['TMP_DIR'], exist_ok=True) + os.makedirs(conf['INPUT_DIR'], exist_ok=True) + proper.print_it = False # To suppress the printing of intermediate steps by PROPER routine + + # To check if input files are there + filename = 'ELT_2048_37m_11m_5mas_nospiders_cut.fits' + my_file = str(conf['INPUT_DIR'] + filename) + if (os.path.isfile(my_file) == False): + get_input_files() + + +def get_input_files(): + id2 = '1YR4G_8E7TpznTxumQA1zD4z_v4V5ZlF2' + dir1 = '.' + name = 'fits_files.zip' + download_from_gdrive(id2, dir1 , name) + + with zipfile.ZipFile(name,'r') as zip_ref: + zip_ref.extractall(conf['INPUT_DIR']) + os.remove(name) + + + + + + + + + + + diff --git a/heeps/config/read_config.py b/heeps/config/read_config.py new file mode 100644 index 0000000..6688c82 --- /dev/null +++ b/heeps/config/read_config.py @@ -0,0 +1,92 @@ +import collections +import os +import proper + +conf = collections.OrderedDict() + + + +# to create directories for inputing and storing fits files +conf['OUT_DIR'] = './output_files/' +conf['TMP_DIR'] = './temp_files/' +conf['INPUT_DIR'] = './input_files/' + + + +# ============================================================================= +# Define parameters for Telescope +# ============================================================================= +conf['GRIDSIZE'] = 1024 # (integar) grid size of the simulation array +conf['WAVELENGTH'] = 5.0*10**-6 # wavelength in micron +conf['DIAM'] = 37.0 # diameter of the telescope in meters +conf['R_OBSTR'] = 0.3 # secondary obstruction in percentage?? +conf['SPIDERS_WIDTH'] = 0.60 # width of ELT spiders in meters +conf['SPIDERS_ANGLE'] = [0, 60.0, 120.0] # angles of spiders +conf['MIS_SEGMENTS_NU'] = 0 # number of mission segments + +# to input pupil as a fits file, +# SCAO team currently uses circular pupil with secondary obstruction +conf['PUPIL_FILE'] = 'ELT_2048_37m_11m_5mas_nospiders_cut.fits' +conf['PREFIX'] = 'test' # prefix for saving fits_files +conf['PIXEL_SCALE'] = 5.0 # plate scale in milli_arc_sec/pixel +conf['F_LENS'] = 658.6 # float, meters, focal distance +conf['DEBUG'] = False # various fits files for coronagraphic propagation is saved in out_dir +conf['DEBUG_PRINT'] = False # prints various values + + +# ============================================================================= +# Parameters for Wavefront abberations +# ============================================================================= +conf['TILT_2D'] = [0.0, 0.0] +conf['TILT_CUBE'] = 10 # creates an array of (n,2) tip/tilt values +conf['ATM_SCREEN_NO'] = 0 # No phase screen +conf['ATM_SCREEN_2D'] = 'metis_370P_35L_HCI_Feb18_rwf8160_cut.fits' # Single phase screen +conf['ATM_SCREEN_CUBE'] = 'cube_atm_1000screens_Feb2018_RandomWind.fits' # 1000 phase screens +conf['GDRIVE_ID'] = '1AUtELRfn_xjnbsMM_SJG6W0c26zyzqNH' # google drive id linked to the fits file + +conf['ISLAND_PISTON'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +conf['STATIC_NCPA'] = True +conf['IMG_LM_SCAO'] = 'NCPA_IMG_LMPP1-SCAO_DET.fits' # NCPA phase screen @ 3.7 um +conf['IMG_NQ_SCAO'] = 'NCPA_IMG_NQPP1-SCAO_DET.fits' # NCPA phase screen @ 10 um + + + +# ============================================================================= +# Parameters for Coronagraphs +# ============================================================================= +""" Types of various modes include: + 1. Vortex coronagraph (VC), to use put conf['MODE'] = 'VC' + 2. Ring apodized vortex coronagraph (RAVC), to use put conf['MODE'] = 'RAVC' + 3. Apodizing phase plate (APP), to use put conf['MODE'] = 'APP' + 4. No coronagraph, just Lyot-stop, to use put conf['MODE'] = 'OFFAXIS' + 5. No coronagraph, but Ring apodizer and LS present, to use put conf['MODE'] = 'MASK' + 6. If conf['MODE'] = anything except above keywords ELT psf is generated +""" + +conf['MODE'] = 'VC' # By default vortex coronagraph +conf['CHARGE'] = 2 # For vortex coronagraph + +conf['PHASE_APODIZER_FILE'] = 'app_phase_cut.fits' +conf['APODIZER_MIS_ALIGN'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +conf['AMP_APODIZER'] = 0 + +conf['LYOT_STOP'] = True +conf['LS_PARA'] = [0.98, 0.03, 1.1] +conf['LS_MIS_ALIGN'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + + +# ============================================================================= +# Detector plane +# ============================================================================= +conf['N_D'] = 512 # final image size + + +collections.OrderedDict(sorted(conf.items())) + +# Creates required directories for simulation +os.makedirs(conf['OUT_DIR'] , exist_ok=True) +os.makedirs(conf['TMP_DIR'], exist_ok=True) +os.makedirs(conf['INPUT_DIR'], exist_ok=True) + + +proper.print_it = False # To suppress the printing of intermediate steps by PROPER routine diff --git a/heeps/coronagraphs/__init__.py b/heeps/coronagraphs/__init__.py new file mode 100644 index 0000000..c1091be --- /dev/null +++ b/heeps/coronagraphs/__init__.py @@ -0,0 +1,7 @@ +from __future__ import absolute_import + +from .circular_apodization import circular_apodization +from .apodization import apodization +from .vortex import vortex +from .lyotstop import lyotstop + diff --git a/heeps/coronagraphs/apodization.py b/heeps/coronagraphs/apodization.py new file mode 100644 index 0000000..b49d64b --- /dev/null +++ b/heeps/coronagraphs/apodization.py @@ -0,0 +1,65 @@ +import numpy as np +import proper +import math +import cv2 +from .circular_apodization import circular_apodization + + + + +def apodization(wf, conf, RAVC=False, phase_apodizer_file=0, amplitude_apodizer_file=0, apodizer_misalignment=0, Debug_print=False): + + r_obstr = conf['R_OBSTR'] + npupil = conf['NPUPIL'] + + apodizer_misalignment = np.array(conf['APODIZER_MIS_ALIGN']) + n = int(proper.prop_get_gridsize(wf)) + apodizer = 1 + if (RAVC == True): + t1_opt = 1. - 1./4*(r_obstr**2 + r_obstr*(math.sqrt(r_obstr**2 + 8.))) # define the apodizer transmission [Mawet2013] + R1_opt = (r_obstr/math.sqrt(1. - t1_opt)) # define the apodizer radius [Mawet2013] + if (Debug_print == True): + print ("r1_opt: ", R1_opt) + print ("t1_opt: ", t1_opt) + apodizer = circular_apodization(wf, R1_opt, 1., t1_opt, xc = apodizer_misalignment[0], yc = apodizer_misalignment[1],NORM=True) # define the apodizer + apodizer = proper.prop_shift_center(apodizer) + + if (isinstance(phase_apodizer_file, (list, tuple, np.ndarray)) == True): + xc_pixels = int(apodizer_misalignment[3]*npupil) + yc_pixels = int(apodizer_misalignment[4]*npupil) + apodizer_pixels = (phase_apodizer_file.shape)[0]## fits file size + scaling_factor = float(npupil)/float(apodizer_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + if (Debug_print==True): + print ("scaling_factor: ", scaling_factor) + apodizer_scale = cv2.resize(phase_apodizer_file.astype(np.float32), (0,0), fx=scaling_factor, fy=scaling_factor, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + if (Debug_print==True): + print ("apodizer_resample", apodizer_scale.shape) + apodizer_large = np.zeros((n,n)) # define an array of n-0s, where to insert the pupuil + if (Debug_print==True): + print("n: ", n) + print("npupil: ", npupil) + apodizer_large[int(n/2)+1-int(npupil/2)-1 + xc_pixels:int(n/2)+1+int(npupil/2)+ xc_pixels,int(n/2)+1-int(npupil/2)-1+ yc_pixels:int(n/2)+1+int(npupil/2)+ yc_pixels] =apodizer_scale # insert the scaled pupil into the 0s grid + apodizer = np.exp(1j*apodizer_large) + + + if (isinstance(amplitude_apodizer_file, (list, tuple, np.ndarray)) == True): + xc_pixels = int(apodizer_misalignment[0]*npupil) + yc_pixels = int(apodizer_misalignment[1]*npupil) + apodizer_pixels = (amplitude_apodizer_file.shape)[0]## fits file size + scaling_factor = float(npupil)/float(apodizer_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + if (Debug_print==True): + print ("scaling_factor: ", scaling_factor) + apodizer_scale = cv2.resize(amplitude_apodizer_file.astype(np.float32), (0,0), fx=scaling_factor, fy=scaling_factor, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + if (Debug_print==True): + print ("apodizer_resample", apodizer_scale.shape) + apodizer_large = np.zeros((n,n)) # define an array of n-0s, where to insert the pupuil + if (Debug_print==True): + print("n: ", n) + print("npupil: ", npupil) + apodizer_large[int(n/2)+1-int(npupil/2)-1 + xc_pixels:int(n/2)+1+int(npupil/2)+ xc_pixels,int(n/2)+1-int(npupil/2)-1+ yc_pixels:int(n/2)+1+int(npupil/2)+ yc_pixels] =apodizer_scale # insert the scaled pupil into the 0s grid + apodizer = apodizer_large + proper.prop_multiply(wf, apodizer) + + return wf + + diff --git a/heeps/coronagraphs/circular_apodization.py b/heeps/coronagraphs/circular_apodization.py new file mode 100644 index 0000000..36578e0 --- /dev/null +++ b/heeps/coronagraphs/circular_apodization.py @@ -0,0 +1,15 @@ +import proper + +def circular_apodization(wf, radius, t_in, t_out, xc = 0.0, yc = 0.0, **kwargs): + + if ("NORM" in kwargs and kwargs["NORM"]): + norm = True + else: + norm = False + + if (t_in > t_out): + apodizer = proper.prop_shift_center(proper.prop_ellipse(wf, radius, radius, xc, yc, NORM = norm))*(t_in-t_out)+t_out + else: + apodizer = proper.prop_shift_center(proper.prop_ellipse(wf, radius, radius, xc, yc, NORM = norm))*(t_out-t_in)+t_in + + return apodizer diff --git a/heeps/coronagraphs/lyotstop.py b/heeps/coronagraphs/lyotstop.py new file mode 100644 index 0000000..6b2426f --- /dev/null +++ b/heeps/coronagraphs/lyotstop.py @@ -0,0 +1,88 @@ +import numpy as np +import proper +import math +import cv2 +from astropy.io import fits + +def lyotstop(wf, conf, RAVC): + input_dir = conf['INPUT_DIR'] + diam = conf['DIAM'] + npupil = conf['NPUPIL'] + spiders_angle = conf['SPIDERS_ANGLE'] + r_obstr = conf['R_OBSTR'] + Debug = conf['DEBUG'] + Debug_print = conf['DEBUG_PRINT'] + LS_amplitude_apodizer_file = conf['AMP_APODIZER'] + LS_misalignment = np.array(conf['LS_MIS_ALIGN']) + if conf['PHASE_APODIZER_FILE'] == 0: + LS_phase_apodizer_file = 0 + else: + LS_phase_apodizer_file = fits.getdata(input_dir+'/'+conf['PHASE_APODIZER_FILE']) + LS = conf['LYOT_STOP'] + LS_parameters = np.array(conf['LS_PARA']) + n = proper.prop_get_gridsize(wf) + if (RAVC==True): # define the inner radius of the Lyot Stop + t1_opt = 1. - 1./4*(r_obstr**2 + r_obstr*(math.sqrt(r_obstr**2 + 8.))) # define the apodizer transmission [Mawet2013] + R1_opt = (r_obstr/math.sqrt(1. - t1_opt)) # define teh apodizer radius [Mawet2013] + r_LS = R1_opt + LS_parameters[1] # when a Ring apodizer is present, the inner LS has to have at least the value of the apodizer radius + else: + r_LS = r_obstr + LS_parameters[1] # when no apodizer, the LS has to have at least the radius of the pupil central obstruction + if LS==True: # apply the LS + if (Debug_print==True): + print("LS parameters: ", LS_parameters) + proper.prop_circular_aperture(wf, LS_parameters[0], LS_misalignment[0], LS_misalignment[1], NORM=True) + proper.prop_circular_obscuration(wf, r_LS, LS_misalignment[0], LS_misalignment[1], NORM=True) + if (LS_parameters[2]!=0): + for iter in range(0,len(spiders_angle)): + if (Debug_print==True): + print("LS_misalignment: ", LS_misalignment) + proper.prop_rectangular_obscuration(wf, LS_parameters[2], 2*diam,LS_misalignment[0], LS_misalignment[1], ROTATION=spiders_angle[iter]) # define the spiders + if (Debug==True): + out_dir = str('./output_files/') + fits.writeto(out_dir +'_Lyot_stop.fits', proper.prop_get_amplitude(wf)[int(n/2)-int(npupil/2 + 50):int(n/2)+int(npupil/2 + 50),int(n/2)-int(npupil/2 + 50):int(n/2)+int(npupil/2 + 50)], overwrite=True) + + if (isinstance(LS_phase_apodizer_file, (list, tuple, np.ndarray)) == True): + xc_pixels = int(LS_misalignment[3]*npupil) + yc_pixels = int(LS_misalignment[4]*npupil) + apodizer_pixels = (LS_phase_apodizer_file.shape)[0]## fits file size + scaling_factor = float(npupil)/float(apodizer_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + +# scaling_factor = float(npupil)/float(pupil_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + if (Debug_print==True): + print ("scaling_factor: ", scaling_factor) + apodizer_scale = cv2.resize(LS_phase_apodizer_file.astype(np.float32), (0,0), fx=scaling_factor, fy=scaling_factor, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + if (Debug_print==True): + print ("apodizer_resample", apodizer_scale.shape) + apodizer_large = np.zeros((n,n)) # define an array of n-0s, where to insert the pupuil + if (Debug_print==True): + print("npupil: ", npupil) + apodizer_large[int(n/2)+1-int(npupil/2)-1 + xc_pixels:int(n/2)+1+int(npupil/2)+ xc_pixels,int(n/2)+1-int(npupil/2)-1+ yc_pixels:int(n/2)+1+int(npupil/2)+ yc_pixels] =apodizer_scale # insert the scaled pupil into the 0s grid + phase_multiply = np.array(np.zeros((n,n)), dtype=complex) # create a complex array + phase_multiply.imag = apodizer_large # define the imaginary part of the complex array as the atm screen + apodizer = np.exp(phase_multiply) + proper.prop_multiply(wf, apodizer) + if (Debug == True): + fits.writeto('LS_apodizer.fits', proper.prop_get_phase(wf), overwrite=True) + + if (isinstance(LS_amplitude_apodizer_file, (list, tuple, np.ndarray)) == True): + print('4th') + xc_pixels = int(LS_misalignment[0]*npupil) + yc_pixels = int(LS_misalignment[1]*npupil) + apodizer_pixels = (LS_amplitude_apodizer_file.shape)[0]## fits file size + scaling_factor = float(npupil)/float(pupil_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + if (Debug_print==True): + print ("scaling_factor: ", scaling_factor) + apodizer_scale = cv2.resize(amplitude_apodizer_file.astype(np.float32), (0,0), fx=scaling_factor, fy=scaling_factor, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + if (Debug_print==True): + print ("apodizer_resample", apodizer_scale.shape) + apodizer_large = np.zeros((n,n)) # define an array of n-0s, where to insert the pupuil + if (Debug_print==True): + print("grid_size: ", n) + print("npupil: ", npupil) + apodizer_large[int(n/2)+1-int(npupil/2)-1 + xc_pixels:int(n/2)+1+int(npupil/2)+ xc_pixels,int(n/2)+1-int(npupil/2)-1+ yc_pixels:int(n/2)+1+int(npupil/2)+ yc_pixels] =apodizer_scale # insert the scaled pupil into the 0s grid + apodizer = apodizer_large + proper.prop_multiply(wf, apodizer) + if (Debug == True): + fits.writeto('LS_apodizer.fits', proper.prop_get_amplitude(wf), overwrite=True) + + return wf diff --git a/heeps/coronagraphs/vortex.py b/heeps/coronagraphs/vortex.py new file mode 100644 index 0000000..67fc49e --- /dev/null +++ b/heeps/coronagraphs/vortex.py @@ -0,0 +1,113 @@ +import numpy as np +import cv2 +import proper +import os +from ..fits import writefield +from ..fits import readfield + + + + + +def vortex(wfo, conf): + tmp_dir = conf['TMP_DIR'] + n = int(proper.prop_get_gridsize(wfo)) + ofst = 0 # no offset + ramp_sign = 1 #sign of charge is positive + ramp_oversamp = 11. # vortex is oversampled for a better discretization + + f_lens = conf['F_LENS'] + diam = conf['DIAM'] + charge = conf['CHARGE'] + pixelsize = conf['PIXEL_SCALE'] + Debug_print = conf['DEBUG_PRINT'] + + if charge!=0: + wavelength = proper.prop_get_wavelength(wfo) + gridsize = proper.prop_get_gridsize(wfo) + beam_ratio = pixelsize*4.85e-9/(wavelength/diam) + calib = str(charge)+str('_')+str(int(beam_ratio*100))+str('_')+str(gridsize) + my_file = str(tmp_dir+'zz_perf_'+calib+'_r.fits') + + proper.prop_propagate(wfo, f_lens, 'inizio') # propagate wavefront + proper.prop_lens(wfo, f_lens, 'focusing lens vortex') # propagate through a lens + proper.prop_propagate(wfo, f_lens, 'VC') # propagate wavefront + + if (os.path.isfile(my_file)==True): + if (Debug_print == True): + print ("Charge ", charge) + vvc = readfield(tmp_dir,'zz_vvc_'+calib) # read the theoretical vortex field + vvc = proper.prop_shift_center(vvc) + scale_psf = wfo._wfarr[0,0] + psf_num = readfield(tmp_dir,'zz_psf_'+calib) # read the pre-vortex field + psf0 = psf_num[0,0] + psf_num = psf_num/psf0*scale_psf + perf_num = readfield(tmp_dir,'zz_perf_'+calib) # read the perfect-result vortex field + perf_num = perf_num/psf0*scale_psf + wfo._wfarr = (wfo._wfarr - psf_num)*vvc + perf_num # the wavefront takes into account the real pupil with the perfect-result vortex field + + else: # CAL==1: # create the vortex for a perfectly circular pupil + if (Debug_print == True): + print ("Charge ", charge) + + wfo1 = proper.prop_begin(diam, wavelength, gridsize, beam_ratio) + proper.prop_circular_aperture(wfo1, diam/2) + proper.prop_define_entrance(wfo1) + proper.prop_propagate(wfo1, f_lens, 'inizio') # propagate wavefront + proper.prop_lens(wfo1, f_lens, 'focusing lens vortex') # propagate through a lens + proper.prop_propagate(wfo1, f_lens, 'VC') # propagate wavefront + + writefield(tmp_dir,'zz_psf_'+calib, wfo1.wfarr) # write the pre-vortex field + nramp = int(n*ramp_oversamp) #oversamp + # create the vortex by creating a matrix (theta) representing the ramp (created by atan 2 gradually varying matrix, x and y) + y1 = np.ones((nramp,), dtype=np.int) + y2 = np.arange(0, nramp, 1.) - (nramp/2) - int(ramp_oversamp)/2 + y = np.outer(y2, y1) + x = np.transpose(y) + theta = np.arctan2(y,x) + x = 0 + y = 0 + vvc_tmp = np.exp(1j*(ofst + ramp_sign*charge*theta)) + theta = 0 + vvc_real_resampled = cv2.resize(vvc_tmp.real, (0,0), fx=1/ramp_oversamp, fy=1/ramp_oversamp, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + vvc_imag_resampled = cv2.resize(vvc_tmp.imag, (0,0), fx=1/ramp_oversamp, fy=1/ramp_oversamp, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + vvc = np.array(vvc_real_resampled, dtype=complex) + vvc.imag = vvc_imag_resampled + vvcphase = np.arctan2(vvc.imag, vvc.real) # create the vortex phase + vvc_complex = np.array(np.zeros((n,n)), dtype=complex) + vvc_complex.imag = vvcphase + vvc = np.exp(vvc_complex) + vvc_tmp = 0. + writefield(tmp_dir,'zz_vvc_'+calib, vvc) # write the theoretical vortex field + + proper.prop_multiply(wfo1, vvc) + proper.prop_propagate(wfo1, f_lens, 'OAP2') + proper.prop_lens(wfo1, f_lens) + proper.prop_propagate(wfo1, f_lens, 'forward to Lyot Stop') + proper.prop_circular_obscuration(wfo1, 1., NORM=True) # null the amplitude iside the Lyot Stop + proper.prop_propagate(wfo1, -f_lens) # back-propagation + proper.prop_lens(wfo1, -f_lens) + proper.prop_propagate(wfo1, -f_lens) + writefield(tmp_dir,'zz_perf_'+calib, wfo1.wfarr) # write the perfect-result vortex field + + vvc = readfield(tmp_dir,'zz_vvc_'+calib) + vvc = proper.prop_shift_center(vvc) + scale_psf = wfo._wfarr[0,0] + psf_num = readfield(tmp_dir,'zz_psf_'+calib) # read the pre-vortex field + psf0 = psf_num[0,0] + psf_num = psf_num/psf0*scale_psf + perf_num = readfield(tmp_dir,'zz_perf_'+calib) # read the perfect-result vortex field + perf_num = perf_num/psf0*scale_psf + wfo._wfarr = (wfo._wfarr - psf_num)*vvc + perf_num # the wavefront takes into account the real pupil with the perfect-result vortex field + + proper.prop_propagate(wfo, f_lens, "propagate to pupil reimaging lens") + proper.prop_lens(wfo, f_lens, "apply pupil reimaging lens") + proper.prop_propagate(wfo, f_lens, "lyot stop") + + return wfo + + + + + + diff --git a/heeps/detector/__init__.py b/heeps/detector/__init__.py new file mode 100644 index 0000000..68278f2 --- /dev/null +++ b/heeps/detector/__init__.py @@ -0,0 +1,5 @@ +from __future__ import absolute_import + + +from .detector import detector + diff --git a/heeps/detector/detector.py b/heeps/detector/detector.py new file mode 100644 index 0000000..cc9039b --- /dev/null +++ b/heeps/detector/detector.py @@ -0,0 +1,24 @@ +import proper +from astropy.io import fits + +def detector(wfo, conf): + f_lens = conf['F_LENS'] + nd = conf['N_D'] + mode = conf['MODE'] + prefix = conf['PREFIX'] + Debug = conf['DEBUG'] + + n = proper.prop_get_gridsize(wfo) + if (n >= nd): + proper.prop_propagate(wfo, f_lens, "to reimaging lens") + proper.prop_lens(wfo, f_lens, "apply reimaging lens") + proper.prop_propagate(wfo, f_lens, "final focus") + (wfo, sampling) = proper.prop_end(wfo, NOABS = False) # conclude the simulation --> noabs= the wavefront array will be complex + else: + print('Error: final image is bigger than initial grid size') + psf = wfo[int(n/2-nd/2):int(n/2+nd/2),int(n/2-nd/2):int(n/2+nd/2)] + out_dir = str('./output_files/') + if (Debug==True): + fits.writeto(out_dir + prefix + '_' + mode+'_PSF'+'.fits', psf, overwrite=True) + return psf + diff --git a/heeps/fits/__init__.py b/heeps/fits/__init__.py new file mode 100644 index 0000000..78904f2 --- /dev/null +++ b/heeps/fits/__init__.py @@ -0,0 +1,6 @@ +from __future__ import absolute_import + + +from .writefield import writefield +from .readfield import readfield + diff --git a/heeps/fits/readfield.py b/heeps/fits/readfield.py new file mode 100644 index 0000000..2a6a61a --- /dev/null +++ b/heeps/fits/readfield.py @@ -0,0 +1,14 @@ +import numpy as np +from astropy.io.fits import getdata + + +def readfield(path, filename): + + data_r, hdr = getdata(path + filename + '_r.fits', header=True) + data_i = getdata(path + filename + '_i.fits') + + field = np.array(data_r, dtype=complex) + field.imag = data_i + + + return(field) diff --git a/heeps/fits/writefield.py b/heeps/fits/writefield.py new file mode 100644 index 0000000..7955cf6 --- /dev/null +++ b/heeps/fits/writefield.py @@ -0,0 +1,11 @@ +from astropy.io import fits + +def writefield(path, filename, field): + + fits.writeto(path + filename + '_r.fits', field.real, header=None, overwrite=True) + fits.writeto(path + filename + '_i.fits', field.imag, header=None, overwrite=True) + + return + + + diff --git a/heeps/metis_hci.py b/heeps/metis_hci.py new file mode 100644 index 0000000..b56de97 --- /dev/null +++ b/heeps/metis_hci.py @@ -0,0 +1,101 @@ +from .pupil import pupil +from .abberations import wavefront_abberations +from .coronagraphs import apodization, vortex,lyotstop +from .detector import detector + + +#from heeps import pupil, wavefront_abberations, detector, apodization, vortex,lyotstop # loads all HEEPS scripts required for simuation +from copy import deepcopy +import numpy as np +from astropy.io import fits +from multiprocessing import Pool +from functools import partial +from sys import platform + +def coronagraphs(wfo, conf): + mode = conf['MODE'] + if mode == 'APP': + RAVC = False + lyotstop(wfo, conf, RAVC) + conf['PHASE_APODIZER_FILE'] = 0 + if mode == 'RAVC': + RAVC = True + apodization(wfo, conf, RAVC=True) + vortex(wfo, conf) + lyotstop(wfo, conf, RAVC) + elif mode == 'VC': + RAVC = False + vortex(wfo, conf) + lyotstop(wfo, conf, RAVC) + elif mode == 'OFFAXIS': + print('No Coronagraph') + RAVC = False + lyotstop(wfo, conf, RAVC) + elif mode == 'MASK': + print('Ring apodizer and LS present') + RAVC = True + apodization(wfo, conf, RAVC=True) + lyotstop(wfo, conf, RAVC) + else: + print('ELT PSF') + return wfo + +def check_dim(input): + r = 0 + try: + r = input.ndim + except: AttributeError + return r + +def propagation_test(a,b): + if check_dim(a)==3 or check_dim(b)==2: + out = 'mutli' + else: + out = 'single' + return out + + +def multi_cube(atm_screen,TILT,conf,wfo1,iter): + wfo = deepcopy(wfo1) + if ((isinstance(atm_screen, (list, tuple, np.ndarray)) == True)): + if (atm_screen.ndim == 3): + atm_screen_iter = atm_screen[iter,:,:] + else: + atm_screen_iter = atm_screen + if (TILT.ndim == 2): + TILT_iter = TILT[iter,:] + else: + TILT_iter = TILT + wavefront_abberations(wfo, conf, atm_screen_iter, TILT_iter) + coronagraphs(wfo, conf) + psf = detector(wfo, conf) + return psf + + +def metis_hci(conf, atm_screen, TILT): + if (propagation_test(atm_screen,TILT)=='single'): + wfo = pupil(conf) + wavefront_abberations(wfo, conf, atm_screen, TILT) + coronagraphs(wfo, conf) + psf = detector(wfo, conf) + fits.writeto(conf['OUT_DIR'] + conf['PREFIX'] + '_PSF_'+ conf['MODE'] +'.fits', psf, overwrite=True) + else: + if (atm_screen.ndim == 3): + length_cube = atm_screen.shape[0] + if (TILT.ndim == 2): + length_cube = TILT.shape[0] + wfo1 = pupil(conf) + + if platform == "linux" or platform == "linux2": + p = Pool(4) + func = partial(multi_cube,atm_screen,TILT,conf,wfo1) + psf_cube = np.array(p.map(func, range(length_cube))) + else: + psf_cube = np.zeros((length_cube, conf['N_D'], conf['N_D'])) + for i in range(length_cube): + psf_cube[i,:,:] = multi_cube(atm_screen,TILT,conf,wfo1,i) + fits.writeto(conf['OUT_DIR'] + conf['PREFIX'] + '_PSF_cube_'+ conf['MODE'] +'.fits', psf_cube, overwrite=True) + psf = psf_cube[0] + return psf + + diff --git a/heeps/pupil/__init__.py b/heeps/pupil/__init__.py new file mode 100644 index 0000000..fae971f --- /dev/null +++ b/heeps/pupil/__init__.py @@ -0,0 +1,4 @@ +from __future__ import absolute_import + +from .pupil import pupil + diff --git a/heeps/pupil/pupil.py b/heeps/pupil/pupil.py new file mode 100644 index 0000000..4f3a686 --- /dev/null +++ b/heeps/pupil/pupil.py @@ -0,0 +1,101 @@ +import numpy as np +import cv2 +import proper +from astropy.io import fits + + +def pupil(conf): + + diam = conf['DIAM'] + gridsize = conf['GRIDSIZE'] + spiders_width = conf['SPIDERS_WIDTH'] + spiders_angle = conf['SPIDERS_ANGLE'] + pixelsize = conf['PIXEL_SCALE'] + + r_obstr = conf['R_OBSTR'] + wavelength = conf['WAVELENGTH'] + missing_segments_number = conf['MIS_SEGMENTS_NU'] + + Debug = conf['DEBUG'] + Debug_print = conf['DEBUG_PRINT'] + + prefix = conf['PREFIX'] + input_dir = conf['INPUT_DIR'] + out_dir = conf['OUT_DIR'] + + beam_ratio = pixelsize*4.85e-9/(wavelength/diam) + wfo = proper.prop_begin(diam, wavelength, gridsize, beam_ratio) + n = int(gridsize) + npupil = np.ceil(gridsize*beam_ratio) # compute the pupil size --> has to be ODD (proper puts the center in the up right pixel next to the grid center) + if npupil % 2 == 0: + npupil = npupil +1 + conf['NPUPIL'] = npupil + if (Debug_print == True): + print ("npupil: ", npupil) + print("lambda: ", wavelength) + pupil_file = fits.getdata(input_dir + conf['PUPIL_FILE']) + if (missing_segments_number == 0): + if (isinstance(pupil_file, (list, tuple, np.ndarray)) == True): + pupil = pupil_file + pupil_pixels = (pupil.shape)[0]## fits file size + scaling_factor = float(npupil)/float(pupil_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + if (Debug_print==True): + print ("scaling_factor: ", scaling_factor) + pupil_scale = cv2.resize(pupil.astype(np.float32), (0,0), fx=scaling_factor, fy=scaling_factor, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + if (Debug_print==True): + print ("pupil_resample", pupil_scale.shape) + pupil_large = np.zeros((n,n)) # define an array of n-0s, where to insert the pupuil + if (Debug_print==True): + print("n: ", n) + print("npupil: ", npupil) + pupil_large[int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2),int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2)] = pupil_scale # insert the scaled pupil into the 0s grid + + proper.prop_circular_aperture(wfo, diam/2) # create a wavefront with a circular pupil + + if (isinstance(pupil_file, (list, tuple, np.ndarray)) == True): + proper.prop_multiply(wfo, pupil_large) # multiply the saved pupil + else: + proper.prop_circular_obscuration(wfo, r_obstr, NORM=True) # create a wavefront with a circular central obscuration + if (spiders_width != 0): + for iter in range(0,len(spiders_angle)): + proper.prop_rectangular_obscuration(wfo, spiders_width, 2*diam, ROTATION = spiders_angle[iter]) # define the spiders + else: + if (missing_segments_number == 1): + pupil = fits.getdata(input_dir+'/ELT_2048_37m_11m_5mas_nospiders_1missing_cut.fits') + if (missing_segments_number == 2): + pupil = fits.getdata(input_dir+'/ELT_2048_37m_11m_5mas_nospiders_2missing_cut.fits') + if (missing_segments_number == 4): + pupil = fits.getdata(input_dir+'/ELT_2048_37m_11m_5mas_nospiders_4missing_cut.fits') + if (missing_segments_number == 7): + pupil = fits.getdata(input_dir+'/ELT_2048_37m_11m_5mas_nospiders_7missing_1_cut.fits') + + pupil_pixels = (pupil.shape)[0]## fits file size + scaling_factor = float(npupil)/float(pupil_pixels) ## scaling factor between the fits file size and the pupil size of the simulation + if (Debug_print==True): + print ("scaling_factor: ", scaling_factor) + pupil_scale = cv2.resize(pupil.astype(np.float32), (0,0), fx=scaling_factor, fy=scaling_factor, interpolation=cv2.INTER_LINEAR) # scale the pupil to the pupil size of the simualtions + if (Debug_print==True): + print ("pupil_resample", pupil_scale.shape) + pupil_large = np.zeros((n,n)) # define an array of n-0s, where to insert the pupuil + if (Debug_print==True): + print("n: ", n) + print("npupil: ", npupil) + pupil_large[int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2),int(n/2)+1-int(npupil/2)-1:int(n/2)+1+int(npupil/2)] =pupil_scale # insert the scaled pupil into the 0s grid + + + proper.prop_multiply(wfo, pupil_large) # multiply the saved pupil + if (spiders_width!=0): + for iter in range(0,len(spiders_angle)): + proper.prop_rectangular_obscuration(wfo, spiders_width, 2*diam,ROTATION=spiders_angle[iter]) # define the spiders + + if (Debug==True): + fits.writeto(out_dir + prefix +'_intial_pupil.fits', proper.prop_get_amplitude(wfo)[int(n/2)-int(npupil/2 + 50):int(n/2)+int(npupil/2 + 50),int(n/2)-int(npupil/2 + 50):int(n/2)+int(npupil/2 + 50)], overwrite=True) + + proper.prop_define_entrance(wfo) #define the entrance wavefront + wfo.wfarr *= 1./np.amax(wfo._wfarr) # max(amplitude)=1 + return wfo + + + + +