</h4>

</h4>

OPTICAL COHERENCE WAVELETS

This Jupyther notebook was designed to illustrate the complete python implementation of the optical coherence wavelets modeling of optical fields under the paraxial and nonparaxial approximations. Some of the central ideas involve intensity distributions, description of the degree of spatial coherence, and marginal power spectrum. You will find details cell by cell in python code will while the theoretical fundamentals are developed in a modular way. Cell execution is done after "shift + enter".

'SPATIAL COHERENCE WAVELETS' is a project under construction with a strong focus on unidimensional diffraction processes under the nonparaxial modeling of optical fields [1], all simulations are based on the theory of spatial coherence wavelets which purpose is the. Further developments are welcome and should point towards the more general 2D case. Enjoy your modeling of optical fields in any state of spatial coherence.

REFERENCES:

[1] Castañeda R, Sucerquia J. Non-approximated numerical modeling of propagation of light in any state of spatial coherence

[2] R. Castañeda,* D. Vargas, E. Franco. Discreteness of the set of radiant point sources: a physical feature of the second-order wave -fronts

[3] R. Castañeda,* D. Vargas, E. Franco. Spatial coherence of light and a fundamental discontinuity of classical second-order wave fronts


In [ ]:
#*********************************************************************
# AUTHOR: David Alejandro Vargas O                                   *
#         david.vargas@geophysik.uni-muenchen.de                     *
# AFFILIATION: LMU Ludwig Maximilian University of Munich            *
# DATE: 01.04.2015                                                   *
#*********************************************************************
#     COPYRIGHT                                                      *
# Copyright (C) 2015  David Alejandro Vargas Otalora                 *
#                                                                    *
# This program is free software: you can redistribute it and/or modif*
# y it under the terms of the GNU General Public License as published*
# by the Free Software Foundation, either version 3 of the License,  *
# or(at your option) any later version.                              *
#                                                                    *
# This program is distributed in the hope that it will be useful,    *
# but WITHOUT ANY WARRANTY; without even the implied warranty of     *
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
# GNU General Public License for more details.                       *
#                                                                    *
# You should have received a copy of the GNU General Public License  *
# along with this program.  If not, see http://www.gnu.org/licenses/ *
#*********************************************************************

outline

0. Import Modules

1. Intensity Distributions.

1.1. Single slide distribution
1.2. Double slide distribution 
1.3. Single pinhole distribution 
1.4. Double pinhole distribution
1.5. Sine distribution 
1.6. Ronchi slide distribution 
1.7. DC Component distribution

2. spatial coherence degree Distributions

2.1. Gaussian distribution    
2.2. Lorentzian distribution

3. Marginal Power Spectrum

3.1. Setup of parameter    
3.2. paraxial Marginal Power Spectrum    
3.3. nonparaxial Marginal Power Spectrum    
3.4. Graphics

4. 2D Modelling of Marginal Power Spectrum

0. Import Modules


In [ ]:
import Intensity
import CoherenceDegree
import MarPowSpectrum
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Show the plots in the Notebook.
plt.switch_backend("nbagg")

1. Intensity Distributions.

Intensity.py is a class which methods will be used as functions to represent the different intensity distributions at the entrance plane.

METHODS:

slide             : Single slide distribution  
double_slide      : Double slide distribution
pinhole           : Pin hole distrubution
double_pinhole    : Double pinhole distribution
sine              : Sine distribution 
ronchi_slide      : Ronchi slide distribution
straight_line     : DC level distribution with slope variation

The whole range of intensity distributions considered here, will be exposed in the following section. The python implementation can be apreciated for each case.

1.1. Single slide distribution


In [ ]:
# SLIDE GRAPHICS
so = 100            # Maximum intensity.
wide = 100          # Slide wide.
xia = np.linspace(-wide, wide, 500)
slide = lambda x: Intensity.slide(x, so, wide)

plt.plot(xia, slide(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ SLIDE $')
plt.axis([-wide, wide, -0.2*so, 1.3*so])
plt.grid(True)
plt.show()

1.2. Double slide distribution


In [ ]:
# DOUBLE SLIDE GRAPHICS
so_1 = 100          # Maximum intensity left slide.
so_2 = 100          # Maximum intensity right slide.
wide_1 = 40         # Left slide wide.
wide_2 = 40         # Right slide wide.
slide_space = 60    # Space between slides.
x = wide_1 + wide_2 + slide_space
xia = np.linspace(-x, x, 500)
double_slide = lambda x: Intensity.double_slide(x, so_1, so_2, wide_1, wide_2, slide_space)

plt.plot(xia, double_slide(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ DOUBLE \ SLIDE $')
plt.axis([-x, x, -0.2*so_1, 1.3*so_1])
plt.grid(True)
plt.show()

1.3. Single pinhole distribution


In [ ]:
# PINHOLE GRAPHICS
so = 100            # Maximum intensity.
wide = 200          # Slide wide.
xia = np.linspace(-wide, wide, 500)
pinhole = lambda x: Intensity.pinhole(x, so, wide)

plt.plot(xia, pinhole(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ PINHOLE $')
plt.axis([-wide, wide, -0.2*so, 1.3*so])
plt.grid(True)
plt.show()

1.4. Double pinhole distribution


In [ ]:
# DOUBLE PINHOLE GRAPHICS
so_1 = 100          # Maximum intensity left slide.
so_2 = 100          # Maximum intensity right slide.
wide_1 = 40         # Left slide wide.
wide_2 = 40         # Right slide wide.
slide_space = 100    # Space between slides.
x = wide_1 + wide_2 + slide_space
xia = np.linspace(-x, x, 500)
double_pinhole = lambda x: Intensity.double_pinhole(x, so_1, so_2, wide_1, wide_2, slide_space)

plt.plot(xia, double_pinhole(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ DOUBLE PINHOLE $')
plt.axis([-x, x, -0.2*so, 1.3*so_1])
plt.grid(True)
plt.show()

1.5. Sine distribution


In [ ]:
# SINE GRAPHICS
so = 200            # Maximum intensity.
wide = 200          # Slide wide.
a = 100             # Amplitude.
w = 4               # Spatial frequency.
xia = np.linspace(-0.5*wide, 0.5*wide, 500)
sine = lambda x: Intensity.sine(x, so, wide, a, w) 

plt.plot(xia, sine(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ SINE $')
plt.axis([-0.5*wide, 0.5*wide, -0.2*so, 2*so])
plt.grid(True)
plt.show()

1.6. Ronchi slide distribution


In [ ]:
# RONCHI SLIDE PINHOLE GRAPHICS
so = 100            # Maximum intensity.
wide = 200          # Slide wide.
w = 4               # Spatial frequency.
xia = np.linspace(-0.5*wide, 0.5*wide, 500)
ronchi_slide = lambda x: Intensity.ronchi_slide(x, so, wide, w)

plt.plot(xia, ronchi_slide(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ RONCHI \ SLIDE $')
plt.axis([-0.5*wide, 0.5*wide, -0.2*so, 1.3*so])
plt.grid(True)
plt.show()

1.7. DC Component distribution


In [ ]:
# STRAIGHT LINE GRAPHICS
so_1 = 50          # Maximum intensity left slide.
so_2 = 200          # Maximum intensity right slide.
wide = 200          # Slide wide.
xia = np.linspace(-wide, wide, 500)
straight_line = lambda x: Intensity.straight_line(x, so_1, so_2, wide)

plt.plot(xia, straight_line(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ STRAIGHT \ LINE $')
plt.axis([-wide, wide, -0.25*so, 2*so_2])
plt.grid(True)
plt.show()

2. spatial coherence degree Distributions

CoherenceDegree.py is a class which methods will be used as functi ons to represent different degrees of spatial coherence distributi ons at the entrance plane.

METHODS:

gauss_coherence     : Gaussian spacial coherence degree 
lorentz_coherence   : Lorentzian spacial coherence degree 


The whole range of coherence degree distributions considered here, will be exposed in the following section. The python implementation can be apreciated for each case.

2.1. Gaussian distribution


In [ ]:
# GAUSSIAN GRAPHICS
sigma = 100     # Gaussian standard deviation.
wide = 400          # Slide wide.
xia = np.linspace(-wide, wide, 100)
gauss_coherence = lambda x: CoherenceDegree.gauss_coherence(x, sigma)

plt.plot(xia, gauss_coherence(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ GAUSSIAN $')
plt.axis([-wide, wide, -0.2, 1.3])
plt.grid(True)
plt.show()

2.2. Lorentzian distribution


In [ ]:
# LORENTZIAN GRAPHICS
gamma = 100     # Lorentzian standard deviation.
wide = 400          # Slide wide.
xia = np.linspace(-wide, wide, 100)
lorentz_coherence = lambda x: CoherenceDegree.lorentz_coherence(x, gamma)

plt.plot(xia, gauss_coherence(xia))
plt.xlabel(r'$\xi_A \ [\mu m]$')
plt.ylabel(r'$Intensity \ I(\xi_A)$')
plt.title(r'$ LORENTZIAN $')
plt.axis([-wide, wide, -0.2, 1.3])
plt.grid(True)
plt.show()

3. Marginal Power Spectrum

3.1. Setup of parameter


In [ ]:
# PARAMETERS
so = 1            # Maximum intensity.
wide = 10          # Slide wide.
sources = 50  # Point sources amount.
pixel_number = 512  # sampling pixel amount in xa 512
wavelength = 0.632  # Optical field wavelength.
z = 5  # Aperture plane - Exit plane distance.
aperture_size = wide  # Aperture size
exit_size = 1.5*wide  # Exit size

# INTENSITY 
slide = lambda x: Intensity.slide(x, so, wide)

# SPACIAL COHERENCE DEGREE 
sigma = 100  # Gaussian standard deviation
gamma = 100
gauss_coherence = lambda x: CoherenceDegree.gauss_coherence(x, sigma)
lorentz_coherence = lambda x: CoherenceDegree.lorentz_coherence(x, gamma)

# COORDINATES
xia = 0.5 * np.linspace(-aperture_size, aperture_size, 2 * sources - 1)  # Entrance coordinate XiA.
xa = 0.5 * np.linspace(-exit_size, exit_size, pixel_number)  # Exit coordinate XA.

# COORDINATE MATRICES OF THE ENTRANCE XiA AND EXIT XA WINDOW
xia_matrix, xa_matrix = np.meshgrid(xia, xa)

# BASE MATRIX
n = np.size(xa)  # number of pixels at the exit window XA
m = 2 * sources - 1  # number of pixels at the entrance window XiA
grid, malla = np.meshgrid(np.arange(1, m + 1), np.arange(1, n + 1))

# WAVELENGTH, DISCRETIZATION OF SLIDE
k = (2 * np.pi) / wavelength
pixel_size = aperture_size / (2 * sources - 2)

intensity = slide
coherence_degree = gauss_coherence

3.2. Paraxial Marginal Power Spectrum


In [ ]:
# COMPUTE REAL LAYER PARAXIAL
real_layer = intensity(xia_matrix) * np.mod(grid, 2)

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

if sources >= 2:
    for family in range(2, sources):
        structure[:, 0:family - 1] = 0  # ojo, verify differences between python and matlab, array initialization
        structure[:, (2 * sources - 1) - (family - 2):] = 0
        
        xid = 2 * (family - 1) * pixel_size
        phase = 0
                         
        rad_1 = np.sqrt(intensity(xia_matrix + 0.5*xid))
        rad_2 = np.sqrt(intensity(xia_matrix - 0.5*xid))
        cosine = np.cos((k/z) * xid * xia_matrix + (k/z) * xid * xa_matrix + phase)
                        
        virtual_source = 2 * coherence_degree(xid) * rad_1 * rad_2 * cosine
        virtual_source = 2 * coherence_degree(xid) * rad_1 * rad_2 * cosine
        
        virtual = structure * virtual_source
        virtual_layer = virtual_layer + virtual  # ojo, do not trust

        structure = structure == 0  # invert zeros for ones.       
        structure = 1*structure    # convert boolean array to int array
else:
    print('Error: np, Point sources amount must be mayor or equal 2')

3.3. Nonparaxial Marginal Power Spectrum

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 [ ]:
# COMPUTE REAL LAYER NONPARAXIAL
radicand = np.power(z, 2) + np.power((xa_matrix - xia_matrix), 2)
lorentz_real = np.power((z + np.sqrt(radicand)) / radicand, 2)
real_layer = intensity(xia_matrix) / (4 * np.power(wavelength, 2)) * lorentz_real * np.mod(grid, 2)

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

if sources >= 2:
    for family in range(2, sources):
        structure[:, 0:family - 1] = 0
        structure[:, (2 * sources - 1) - (family - 2):] = 0
        
        xid = 2 * (family - 1) * pixel_size

        term_1 = np.power(z, 2) + np.power(xa_matrix - xia_matrix, 2) + (np.power(xid, 2) / 4)
        
        radicand_1 = term_1 - (xa_matrix - xia_matrix) * xid
        radicand_2 = term_1 + (xa_matrix - xia_matrix) * xid

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

        rad_1 = np.sqrt(intensity(xia_matrix + 0.5*xid))
        rad_2 = np.sqrt(intensity(xia_matrix - 0.5*xid))
        #print(xa_matrix + 0.5*xid)

        lorentz_virtual = (1 / (2 * np.power(wavelength, 2))) * rad_1 * rad_2 * lor_1 * lor_2
        cosine = np.cos(k * (np.sqrt(radicand_1) - np.sqrt(radicand_2)))
        
        virtual_source = coherence_degree(xid) * lorentz_virtual * cosine
        
        virtual = structure * virtual_source
        virtual_layer = virtual_layer + virtual  # ojo, do not trust

        structure = structure == 0  # invert zeros for ones.
        structure = 1*structure    # convert boolean array to int array
else:
    print('Error: np, Point sources amount must be mayor or equal 2')

3.4. Graphics


In [ ]:
# MARGINAL POWER SPECTRUM (MPS)
Mar_Pow_Spectrum = real_layer + virtual_layer

s_xaReal = real_layer.sum(axis=1)
s_xa = Mar_Pow_Spectrum.sum(axis=1)
s_xia = Mar_Pow_Spectrum.sum(axis=0)

#--------------------------------------------------------------
# MARGINAL POWER SPECTRUM, RADIANT AND VIRTUAL PART
#--------------------------------------------------------------
plt.figure(figsize=(12,4), dpi=80) 
plt.subplot(1, 3, 1)
plt.imshow(real_layer, cmap='magma', aspect='auto',
           extent=[-0.5*aperture_size, 0.5*aperture_size,
                   -0.5*exit_size, 0.5*exit_size])
plt.title('MPS - Real Part')
plt.xlabel(r'$ \xi_A \ (\mu m)$')
plt.ylabel(r'$ r_A \ (\mu m)$')

plt.subplot(1, 3, 2)
plt.imshow(virtual_layer, cmap='inferno', aspect='auto',
           extent=[-0.5*aperture_size, 0.5*aperture_size,
                   -0.5*exit_size, 0.5*exit_size])
plt.title('MPS - Virtual Part')
plt.xlabel(r'$ \xi_A \ (\mu m)$')
plt.ylabel(r'$ r_A \ (\mu m)$')

plt.subplot(1, 3, 3)
plt.imshow(Mar_Pow_Spectrum, cmap='jet', aspect='auto',
           extent=[-0.5*aperture_size, 0.5*aperture_size,
                   -0.5*exit_size, 0.5*exit_size])
plt.title('Marginal Power Spectrum')
plt.xlabel(r'$ \xi_A \ (\mu m)$')
plt.ylabel(r'$ r_A \ (\mu m)$')

plt.show()

In [ ]:
#--------------------------------------------------------------
# MARGINAL POWER SPECTRUM, ENTRANCE AND EXIT POWER SPECTRUM
#--------------------------------------------------------------
plt.figure(figsize=(12,4), dpi=80) 
plt.subplot(1, 3, 1)
plt.imshow(Mar_Pow_Spectrum, cmap='gray', aspect='auto',
           extent=[-0.5*aperture_size, 0.5*aperture_size,
                   -0.5*exit_size, 0.5*exit_size])
plt.title('Marginal Power Spectrum')
plt.xlabel(r'$ \xi_A \ (\mu m)$')
plt.ylabel(r'$ r_A \ (\mu m)$')

plt.subplot(1, 3, 2)
plt.fill(xia, s_xia, 'c', alpha=0.9)
plt.title('Entrance Power Spectrum')
plt.xlabel(r'$ \xi_A \ (\mu m)$')
plt.ylabel(r'$ S(\xi_A) $')

plt.subplot(1, 3, 3)
plt.fill(xa, s_xa, 'r', alpha=0.7)
plt.title('Exit Power Spectrum')
plt.xlabel(r'$ r_A \ (\mu m) $')
plt.ylabel(r'$ S(r_A) $')

plt.show()

4. 2D Modelling of Marginal Power Spectrum (UNFINISHED WORK)

4.1. 3D - Structure


In [ ]:
# PARAMETERS
sources = 2  # Point sources amount.

# BASE MATRIX
n = 2 * sources - 1  # number of pixels at the entrance window XiA
grid, malla = np.meshgrid(np.arange(1, n + 1), np.arange(1, n + 1))

Structure = np.mod(grid, 2) * np.transpose(np.mod(grid, 2))

print Structure

In [ ]:
import numpy as np

# PARAMETERS
sources = 2  # Point sources amount.
pixel_number = 3  # sampling pixel amount in xa
aperture_size = 4  # Aperture size
exit_size = 10  # Exit size

# COORDINATES
xia1 = 0.5 * np.linspace(-aperture_size, aperture_size, 2 * sources - 1)  # Entrance coordinate XiA.
xia2 = 0.5 * np.linspace(-aperture_size, aperture_size, 2 * sources - 1)  # Entrance coordinate XiA.
xa = 0.5 * np.linspace(-exit_size, exit_size, pixel_number)  # Exit coordinate XA.

# COORDINATE MATRICES OF THE ENTRANCE XiA AND EXIT XA WINDOW
xia1_matrix, xa_matrix , xia2_matrix = np.meshgrid(xia1, xa, xia2)
xia_matrix = np.sqrt(np.power(xia1_matrix,2) + np.power(xia2_matrix,2))  # Lo logré!
#print xia_matrix



# BASE MATRIX
#n = np.size(xa)  # number of pixels at the exit window XA
m = 2 * sources - 1  # number of pixels at the entrance window XiA
n = m
grid1, grid2, grid3 = np.meshgrid(np.arange(1, m + 1), np.arange(1, n + 3), np.arange(1, m + 1))
Structure = np.mod(grid1, 2) * np.mod(grid3, 2)  # Lo logré!
#print Structure