new architecture

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")
+

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()
+

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
+
+

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
+

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

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

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

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

heeps/adi/__init__.py

@@ -0,0 +1,3 @@
+from __future__ import absolute_import
+
+

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

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)