### Real layer Analysis

$$\textbf{S}\left(\boldsymbol\xi_{A},\boldsymbol r_{A};\nu\right) = \frac{1}{4\lambda^2}S_{o}\left(\boldsymbol\xi_{A};\nu\right) \left(\frac{z+\sqrt{z^2+|\boldsymbol r_{A}-\boldsymbol\xi_{A}|^2}} {z^2+|\boldsymbol r_{A}-\boldsymbol\xi_{A}|^2}\right)^2$$$$|\boldsymbol r_{A}-\boldsymbol\xi_{A}|^2 = (r_{A_{1}}-\xi_{A_{1}})^2 + (r_{A_{2}}-\xi_{A_{2}})^2$$


In [5]:

import numpy as np
import matplotlib.pyplot as plt

# COMPUTE REAL LAYER
intensity = 1
wavelength = 0.632
z = 1

sources1 = 4           # Point sources amount in xi_A1.
sources2 = 4           # Point sources amount in xi_A2.
pixel_number = 512     # sampling pixel amount in r_A.
alpha = np.sqrt(0.01)  # consider 99% of the free space diffraction envelope

AP_width1 = 10
AP_width2 = 10
OP_width1 = np.sqrt(np.power((1+np.sqrt(1+8*alpha))/(4*alpha),2)-1)*z + AP_width1
OP_width2 = np.sqrt(np.power((1+np.sqrt(1+8*alpha))/(4*alpha),2)-1)*z + AP_width2

xi_A1 = np.linspace(-0.5 * AP_width1, 0.5 * AP_width1, 2*sources1 - 1)  # AP Axis
xi_A2 = np.linspace(-0.5 * AP_width2, 0.5 * AP_width2, 2*sources2 - 1)  # AP Axis
r_A1 = np.linspace(-0.5 * OP_width1, 0.5 * OP_width1, pixel_number)    # OP Axis
r_A2 = np.linspace(-0.5 * OP_width2, 0.5 * OP_width2, pixel_number)    # OP Axis
#-----------------------------------------------------------------------------------------------

# BASE MATRIX
nf1 = int((np.size(xi_A1) + 1)/2)
nf2 = int((np.size(xi_A2) + 1)/2)
m1 = 2 * nf1 - 1
m2 = 2 * nf2 - 1
n1 = np.size(r_A1)
n2 = np.size(r_A2)
grid1, grid2, G2, G3 = np.meshgrid(np.arange(1, m1 + 1), np.arange(1, m2 + 1),
np.arange(1, n1 + 1), np.arange(1, n2 + 1))

# COORDINATE MATRICES OF THE APERTURE xi_A AND Observation r_A PLANE
xi_A1_matrix, xi_A2_matrix, r_A1_matrix, r_A2_matrix = np.meshgrid(xi_A1, xi_A2, r_A1, r_A2)

# WAVELENGTH, DISCRETIZATION OF APERTURE
k = (2 * np.pi) / wavelength
AP_width1 = 2 * xi_A1.max()   # AP width
AP_width2 = 2 * xi_A2.max()   # AP width
pixel_size1 = AP_width1 / (2 * nf1 - 2)
pixel_size2 = AP_width2 / (2 * nf2 - 2)

# COMPUTE REAL LAYER
radicand = np.power(z, 2) + np.power(r_A1_matrix - xi_A1_matrix, 2) + np.power(r_A2_matrix - xi_A2_matrix, 2)
factor = intensity / (4 * np.power(wavelength, 2))
real_layer = factor * lorentz_real * np.mod(grid1, 2)* np.mod(grid2, 2)
#-----------------------------------------------------------------------------------------------------------------




In [10]:

# COMPUTE VIRTUAL LAYER
virtual_layer = np.zeros((m1, m2, n1, n2))
structure = np.mod(grid1, 2)* np.mod(grid2, 2) <= 0
structure = 1*structure    # convert boolean array to int array

if nf >= 2:
for family in range(2, nf+1):
structure[:, 0:family - 1] = 0
structure[:, (2 * nf - 1) - (family - 2):] = 0

xi_d = 2 * (family - 1) * pixel_size

term_1 = np.power(z, 2)
term_2 = np.power(r_A_matrix - xi_A_matrix, 2)
term_3 = (np.power(xi_d, 2) / 4)
term = term_1 + term_2 + term_3

radicand_1 = term - (r_A_matrix - xi_A_matrix) * xi_d
radicand_2 = term + (r_A_matrix - xi_A_matrix) * xi_d

factor2 = (1 / (2 * np.power(wavelength, 2)))
cosine = np.cos(phase)

source = coherence_degree(xi_d) * lorentz_virtual * cosine

virtual = structure * source
virtual_layer = virtual_layer + virtual

structure = structure == 0  # invert zeros for ones.
# convert boolean array to int array
structure = 1*structure

else:
raise ValueError('np must be mayor or equal 2')
return real_layer, virtual_layer




In [6]:

# POWER SPECTRUM
Sr_AReal = real_layer.sum(axis=1)
Sr_AReal = Sr_AReal.sum(axis=0)

# GRAPHICS
plt.imshow(Sr_AReal, cmap='gray', aspect='auto',
extent=[-0.5*OP_width1, 0.5*OP_width1,
-0.5*OP_width2, 0.5*OP_width2])
plt.title('Observation plane - Power Spectrum')
plt.xlabel(r'$r_{A1} \ (\mu m)$')
plt.ylabel(r'$r_{A2} \ (\mu m)$')

plt.show()




In [73]:

def nonparaxial(xi_A, r_A, z, wavelength, intensity, coherence_degree):
"""
# NOTE: All lenght units must be in micrometers.
:param xi_A:               Entrance coordinate xi_A.
:param r_A:                Observation coordinate r_A.
:param wavelength:         Optical field wavelength.
:param z:                  Aperture plane - Observation plane distance.
:param intensity:          Intensity distribution
:param coherence_degree:   Coherence degree
:return: real_layer:       Marginal Power Spectrum's real layer
:return: virtual_layer:    Marginal Power Spectrum's virtual layer
"""

# BASE MATRIX
nf = int((np.size(xi_A) + 1)/2)   # Point sources amount.
n = np.size(r_A)  # number of pixels at the Observation plane OP
m = 2 * nf - 1  # number of pixels at the aperture plane AP
grid, nonused = np.meshgrid(np.arange(1, m + 1), np.arange(1, n + 1))

# COORDINATE MATRICES OF THE APERTURE xi_A AND Observation r_A PLANE
xi_A_matrix, r_A_matrix = np.meshgrid(xi_A, r_A)

# WAVELENGTH, DISCRETIZATION OF SLIDE
k = (2 * np.pi) / wavelength
AP_width = 2 * xi_A.max()   # AP width
pixel_size = AP_width / (2 * nf - 2)

# COMPUTE REAL LAYER
radicand = np.power(z, 2) + np.power((r_A_matrix - xi_A_matrix), 2)
factor1 = intensity(xi_A_matrix) / (4 * np.power(wavelength, 2))
real_layer = factor1 * lorentz_real * np.mod(grid, 2)

# COMPUTE VIRTUAL LAYER
virtual_layer = np.zeros((n, m))
structure = np.mod(grid, 2) <= 0
structure = 1*structure    # convert boolean array to int array

if nf >= 2:
for family in range(2, nf+1):
structure[:, 0:family - 1] = 0
structure[:, (2 * nf - 1) - (family - 2):] = 0

xi_d = 2 * (family - 1) * pixel_size

term_1 = np.power(z, 2)
term_2 = np.power(r_A_matrix - xi_A_matrix, 2)
term_3 = (np.power(xi_d, 2) / 4)
term = term_1 + term_2 + term_3

radicand_1 = term - (r_A_matrix - xi_A_matrix) * xi_d
radicand_2 = term + (r_A_matrix - xi_A_matrix) * xi_d

factor2 = (1 / (2 * np.power(wavelength, 2)))
cosine = np.cos(phase)

source = coherence_degree(xi_d) * lorentz_virtual * cosine

virtual = structure * source
virtual_layer = virtual_layer + virtual

structure = structure == 0  # invert zeros for ones.
# convert boolean array to int array
structure = 1*structure

else:
raise ValueError('np must be mayor or equal 2')
return real_layer, virtual_layer




In [ ]: