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.
[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/ *
#*********************************************************************
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.1. Gaussian distribution
2.2. Lorentzian distribution
3.1. Setup of parameter
3.2. paraxial Marginal Power Spectrum
3.3. nonparaxial Marginal Power Spectrum
3.4. Graphics
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")
Intensity.py is a class which methods will be used as functions to represent the different intensity distributions at the entrance plane.
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.
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()
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()
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()
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()
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()
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()
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()
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.
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.
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()
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()
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
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')
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')
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()
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