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)
lorentz_real = np.power((z + np.sqrt(radicand)) / radicand, 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

        lor_1 = (z + np.sqrt(radicand_1)) / radicand_1
        lor_2 = (z + np.sqrt(radicand_2)) / radicand_2

        rad_1 = np.sqrt(intensity(xi_A_matrix + 0.5*xi_d))
        rad_2 = np.sqrt(intensity(xi_A_matrix - 0.5*xi_d))
            
        factor2 = (1 / (2 * np.power(wavelength, 2)))
        lorentz_virtual = factor2 * rad_1 * rad_2 * lor_1 * lor_2
        phase = k * (np.sqrt(radicand_1) - np.sqrt(radicand_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)
    lorentz_real = np.power((z + np.sqrt(radicand)) / radicand, 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

            lor_1 = (z + np.sqrt(radicand_1)) / radicand_1
            lor_2 = (z + np.sqrt(radicand_2)) / radicand_2

            rad_1 = np.sqrt(intensity(xi_A_matrix + 0.5*xi_d))
            rad_2 = np.sqrt(intensity(xi_A_matrix - 0.5*xi_d))
            
            factor2 = (1 / (2 * np.power(wavelength, 2)))
            lorentz_virtual = factor2 * rad_1 * rad_2 * lor_1 * lor_2
            phase = k * (np.sqrt(radicand_1) - np.sqrt(radicand_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 [ ]: