Skip to content

GitLab

Commit 1e7805ce authored by Alexis Carlotti's avatar Alexis Carlotti
Browse files

small changes to MAIN.py

parent 93513be2
Hide whitespace changes
Inline Side-by-side
hrm-hc/utils_basic.py
View file @ 1e7805ce
......@@ -16,6 +16,62 @@ import warnings
warnings.simplefilter('ignore', category=AstropyWarning)
import imageio
WORKPATH = '/Users/carlotal/harmoni-hc/hrm-hc/COMMON_FILES/'
def ZELDA(P,PHI,lbd_vect,k_lbd,delta,D_win,D_tel,zelda_diam_mas,zelda_height,ref_index):
#P: pupil
#PHI: phase screen (m)
#lbd: wavelength vector (m)
#k_lbd: index within wavelength vector
#delta: dispersion from lbd_0 (mas)
#D_win: physical diameter of input array (m)
#D_tel: physical diameter of telescope (m)
#zelda_diam_mas: diameter of zelda mask (mas)
#zelda height: height of zelda mask (m)
#ref_index: refraction index at lbd
# ____________ _ _____
# |___ / ____| | | __ \ /\
# / /| |__ | | | | | | / \
# / / | __| | | | | | |/ /\ \
# / /__| |____| |____| |__| / ____ \
# /_____|______|______|_____/_/ \_\
N = len(P)
x,X,Y,R,T = VECT(N,D_win)
lbd = lbd_vect[k_lbd]
lbd_0 = np.mean(lbd_vect)
N_lbd = len(lbd_vect)
dispersion = np.exp(2*1j*np.pi*MAS2LD(delta,D_tel,lbd)*Y/D_tel) #phase due to dispersion
D_ratio = D_win/D_tel
lbd_ratio = lbd/lbd_0
FOV = 1.1*MAS2LD(zelda_diam_mas,D_tel,lbd)
# computation of electric fields in image plane
E_zelda_pre = ft(P*np.exp(2*1j*np.pi*PHI/lbd)*dispersion,1,D_ratio,D_ratio,FOV,FOV,lbd_ratio,N,N)
E_zelda_full = ft(P*np.exp(2*1j*np.pi*PHI/lbd)*dispersion,1,D_ratio,D_ratio,N/2,N/2,lbd_ratio,N,N)
E_zelda_noab = ft(P*dispersion,1,D_ratio,D_ratio,FOV,FOV,lbd_ratio,N,N) #no aberrations here; used to infer the diffraction due to the mask
# zelda mask model
u,U,V,R,T = VECT(N,FOV)
zelda_shape = (R<MAS2LD(zelda_diam_mas,D_tel,lbd)/2)
zelda_mask = np.exp(2*1j*np.pi*zelda_height*(ref_index-1)/lbd)*zelda_shape
E_zelda_post = E_zelda_pre*zelda_mask
E_zelda_diff = E_zelda_pre*zelda_shape
# computation of electric field in second pupil plane
P_zelda = ft(E_zelda_post,1,FOV,FOV,D_ratio,D_ratio,lbd_ratio,N,N)
P_zelda_diff = ft(E_zelda_diff,1,FOV,FOV,D_ratio,D_ratio,lbd_ratio,N,N)
P_zelda_full = ft(E_zelda_full,1,N/2,N/2,D_ratio,D_ratio,lbd_ratio,N,N)
P_camera = P_zelda_full + P_zelda - P_zelda_diff
# intensity in pupil plane
I = np.abs(P_camera)**2/N_lbd
# b factor (diffraction of zelda mask)
b = -np.real(ft(E_zelda_noab*zelda_shape,1,FOV,FOV,D_ratio,D_ratio,lbd_ratio,N,N))/N_lbd #for later calibration
return (I,b)
def MyInterp2(IN,x1,x2):
interp_f = interp2d(x1,x1,IN,kind='linear',fill_value=0)
OUT = interp_f(x2,x2)
return OUT
def BOUNDARIES1D(A):
A = (A.round()).astype(int)
# "Enclose" mask with sentients to catch shifts later on
......@@ -633,8 +689,8 @@ def AB2V_CONVFACT(WORKPATH,lambda_0, BW):
DATA_V=fits.getdata(WORKPATH + 'VEGA_Fnu.fits')
print(np.shape(DATA_V))
lbd=np.linspace(lambda_0-BW/2,lambda_0+BW/2,1000)
dlbd=lbd(2)-lbd(1)
F_nu_f = interp1d(DATA_V[:,1],DATA_V[:,2])
dlbd=lbd[1]-lbd[0]
F_nu_f = interp1d(DATA_V[:,0],DATA_V[:,1],fill_value=0)
F_nu = F_nu_f(lbd)
S1=np.sum(F_nu*dlbd*3.34e-19*(lbd*10)**2);
S2=np.sum(dlbd*np.ones(1000))
......@@ -654,7 +710,7 @@ def PhotonCount(T,S,lambda_0,BW,mag,mag_type,throughput):
if mag_type == 'AB':
AB_mag=mag
elif mag_type == 'VEGA':
AB_mag=mag+AB2V_CONVFACT(lambda_0, BW)
AB_mag=mag+AB2V_CONVFACT(WORKPATH,lambda_0*1e9, BW*1e9)
else:
AB_mag=mag
print('Error! Choose between AB mag or VEGA. Switching to AB mag')
......@@ -665,28 +721,28 @@ def PhotonCount(T,S,lambda_0,BW,mag,mag_type,throughput):
def MakeBinary(A):
N = len(A)
B = A
for j in range(N):
for i in range(N):
old_pixel = A[i,j]
new_pixel = np.round(A[i,j])
A[i,j] = new_pixel
old_pixel = B[i,j]
new_pixel = np.round(B[i,j])
B[i,j] = new_pixel
quant_error = old_pixel-new_pixel
if ((i+1)<N-1):
A[i+1,j]=A[i+1,j]+7/16*quant_error
B[i+1,j]=B[i+1,j]+7/16*quant_error
if (i>0 and j<N-1):
A[i-1,j+1]=A[i-1,j+1]+3/16*quant_error
B[i-1,j+1]=B[i-1,j+1]+3/16*quant_error
if ((j+1)<N-1):
A[i,j+1]=A[i,j+1]+5/16*quant_error
B[i,j+1]=B[i,j+1]+5/16*quant_error
if ((i+1)<N-1 and (j+1)<N-1):
A[i+1,j+1]=A[i+1,j+1]+1/16*quant_error
return A
B[i+1,j+1]=B[i+1,j+1]+1/16*quant_error
return B
def AZAV(I,OWA,PhotAp):
L = len(I)
......@@ -740,6 +796,8 @@ def MAKE_ELT(WORKPATH,N,D,n_miss,ref_err):
for k in range(0,n_miss):
MS *= (1-(color==MISS_LIST[k]))
#
hexa_nospider = hexa
hexa_nospider[R<15]=1
t=0.5
hexa=MS*hexa*(np.abs(X)>t/2)*((Y<(np.tan(np.pi/6)*X-t/2/np.cos(np.pi/6)))+(Y>(np.tan(np.pi/6)*X+t/2/np.cos(np.pi/6))))*((Y<(np.tan(-np.pi/6)*X-t/2/np.cos(-np.pi/6)))+(Y>(np.tan(-np.pi/6)*X+t/2/np.cos(-np.pi/6))))
......@@ -752,6 +810,27 @@ def MAKE_ELT(WORKPATH,N,D,n_miss,ref_err):
#
return (hexa,color,MS,referr)
def MAKE_ELT_CALUNIT(WORKPATH,N,D):
#Problem of coherence between ESO pupil and this one
#D = D *1.021
C = fits.getdata(WORKPATH + 'Coord_ELT.fits')
#
N=N+N%2
W=1.45*np.cos(np.pi/6.0)
[x,X,Y,R,T]=VECT(N,D)
#
hexa=np.zeros([N,N])
#
for k in range(0,len(C)):
Xt=X+C[k,0]
Yt=Y+C[k,1]
hexa=hexa+(k+1)*(Yt<0.5*W)*(Yt>=-0.5*W)*((Yt+np.tan(np.pi/3)*Xt)<W)*((Yt+np.tan(np.pi/3)*Xt)>=-W)*((Yt-np.tan(np.pi/3)*Xt)<W)*((Yt-np.tan(np.pi/3)*Xt)>=-W)
#
hexa=np.double(hexa>0)
hexa[R<15]=1
return hexa
def EELTNLS(N,a1,a2,a3):
D=1
d=0.3
......@@ -884,20 +963,25 @@ def MPS(N,power):
screen = np.real(screen*KER)
return screen
def MPS_DSP(N,DSP):
def MPS_DSP(M,DSP):
# generates a random phase screen of NxN points based on a 2D array representing the PSD
# The extent of the PSD must be 2*f_max where f_max is the maximum frequency is N/2
N = len(DSP)
f_max = N/2
f_vect,Xf,Yf,Rf,Tf = VECT(N,2*f_max)
screen = np.fft.fft2(np.exp(2*1j*pi*(np.random.rand(N,N)-0.5))*np.sqrt(DSP))
KER = [[-1,1],[1,-1]]
if N//2-N/2==0:
KER = np.pad(KER,[np.int(N/2-1),np.int(N/2-1)],'wrap')
else:
KER = np.pad(KER,[np.int(N/2),np.int(N/2)],'wrap')
KER = KER[0:N,0:N]
screen = np.real(screen*KER)
return screen
rand_screen = np.exp(2*1j*pi*(np.random.rand(N,N)-0.5))
#screen = np.fft.fftshift(np.fft.fft2(rand_screen*np.sqrt(DSP)))
screen = ft(rand_screen*np.sqrt(DSP),1,2*f_max,2*f_max,1,1,1,M,M)*(M/N)
#KER = [[-1,1],[1,-1]]
#if N//2-N/2==0:
# KER = np.pad(KER,[np.int(N/2-1),np.int(N/2-1)],'wrap')
#else:
# KER = np.pad(KER,[np.int(N/2),np.int(N/2)],'wrap')
# KER = KER[0:N,0:N]
#screen = np.real(screen*KER)
screen = np.real(screen)
#screen_ft = np.real(screen_ft)
return screen#,screen_ft)
def DSP(WF,gamma,P,alpha_min,alpha_max):
#WF: wavefront
......@@ -936,6 +1020,13 @@ def MAKE_WF_CUBE(N,D,WFE):
def DISPER(lbd_min,lbd_max,theta):
#lbd_min & lbd_max must be given in microns
#theta must be given in degrees
if lbd_min > lbd_max:
new_lbd_max = lbd_min
lbd_min = lbd_max
lbd_max = new_lbd_max
sign = -1
else:
sign = 1
TC = 12
T = TC + 273.16
RH = 15
......@@ -950,12 +1041,14 @@ def DISPER(lbd_min,lbd_max,theta):
S = 1.0/lbd
N0_1 = 1.0e-8*((2371.34+683939.7/(130-S0**2)+4547.3/(38.9-S0**2))*D1+(6487.31+58.058*S0**2-0.71150*S0**4+0.08851*S0**6)*D2)
N_1 = 1.0e-8*((2371.34+683939.7/(130-S**2)+4547.3/(38.9-S**2))*D1+(6487.31+58.058*S**2-0.71150*S**4+0.08851*S**6)*D2)
D = np.max(np.tan(theta*pi/180)*(N0_1-N_1)*206264.8)
D = sign*np.max(np.tan(theta*pi/180)*(N0_1-N_1)*206264.8)
return D
def DISPER_ADC(lbd_min,lbd_max,theta,theta_min,theta_max):
D = DISPER(lbd_min,lbd_max,theta)
D = np.abs(D-(DISPER(lbd_min,lbd_max,theta_min)+DISPER(lbd_min,lbd_max,theta_max))/2)
if lbd_min > lbd_max:
D = -D
return D
def ELEVATION(D,H):
......@@ -1040,6 +1133,17 @@ def ft(Ein, z, a, b, u, v, lbd, nu, nv):
Eout = Eout*dx*dx*np.exp(2*1j*pi*z/lbd)/(1j*lbd*z)
return Eout
def ft_cyl(Ein, z, a, b, u, v, lbd, nu, nv):
nx, ny = np.shape(Ein)
dx = a/nx
du = u/nu
x = (np.linspace(-nx/2,nx/2,nx,endpoint=False)+0.5)*dx
u = (np.linspace(-nu/2,nu/2,nu,endpoint=False)+0.5)*du
x_u = np.outer(x,u)
Eout = np.dot(np.transpose(np.exp(-2.*pi*1j*x_u/lbd)),Ein)
Eout = Eout*dx*np.exp(2*1j*pi*z/lbd)/np.sqrt(1j*lbd*z)
return Eout
def ftinv(Ein, z, a, b, u, v, lbd, nu, nv):
nx, ny = np.shape(Ein)
dx = a/nx
......@@ -1108,7 +1212,17 @@ def ftAMPL(Ein, a, nu, u, lbd):
Eout = np.dot(Eout,np.exp(-2.*pi*1j*x_u/lbd))
Eout = Eout/lbd*dx*dx
return Eout
def ft1D(Ein, z, a, u, lbd, nu):
nx = len(Ein)
dx = a/nx
du = u/nu
x = (np.linspace(-nx/2,nx/2,nx,endpoint=False)+0.5)*dx
u = (np.linspace(-nu/2,nu/2,nu,endpoint=False)+0.5)*du
x_u = np.outer(x,u)
Eout = np.exp(1j*np.pi*z/lbd)/(1j*np.sqrt(lbd*z))*np.dot(np.exp(-2*np.pi*1j*x_u/lbd/z),Ein)*dx
return Eout
def Fresnel(Ein, z, a, u, lbd, nu, trigger):
if trigger==False:
nx = len(Ein)
......@@ -1121,10 +1235,12 @@ def Fresnel(Ein, z, a, u, lbd, nu, trigger):
Eout = ft(Ein*np.exp(1j*pi/(lbd*z)*(X**2+Y**2)),z,a,a,u,u,lbd,nu,nu)*np.exp(1j*pi/(lbd*z)*(U**2+V**2))
else:
N = len(Ein)
four,x_int,dx_int = ft_BASIC(Ein,a,N/a,N,1)
four = ft_BASIC(Ein,a,N/a,N,1)
dx_int = 1/a
x_int = (np.linspace(-N/2,N/2,N,endpoint=False)+0.5)*dx_int
X_int,Y_int = np.meshgrid(x_int,x_int,sparse=False)
angular = 1j*lbd*z*np.exp(-1j*pi*z*lbd*(X_int**2+Y_int**2))
Eout,x,dx= ft_BASIC(angular*four,N/a,u,nu,-1)
Eout = ft_BASIC(angular*four,N/a,u,nu,-1)
Eout = Eout*np.exp(2*1j*pi*z/lbd)/(1j*lbd*z)
return Eout
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment