# !/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
This module generates Scalar_source_XY class for defining sources.
Its parent is Scalar_field_XY.
The main atributes are:
* self.x - x positions of the field
* self.y - y positions of the field
* self.u - field XY
* self.wavelength - wavelength of the incident field. The field is monocromatic
The magnitude is related to microns: `micron = 1.`
*Class for unidimensional scalar masks*
*Functions*
* plane_wave
* gauss_beam
* spherical_wave
* vortex_beam
* laguerre_beam
* hermite_gauss_beam
* zernike_beam
* bessel_beam
* plane_waves_dict
* plane_waves_several_inclined
* gauss_beams_several_parallel
* gauss_beams_several_inclined
*Also*
* laguerre_polynomial_nk
* fZernike
* delta_kronecker
"""
from math import factorial
from numpy import arctan2, cos, exp, pi, sign, sin, sqrt, zeros
from scipy.special import eval_hermite, j0, j1, jv
from . import degrees, np, um
from .scalar_fields_XY import Scalar_field_XY
from .utils_math import fZernike, laguerre_polynomial_nk
# from scipy.special.orthogonal import hermite
[docs]
class Scalar_source_XY(Scalar_field_XY):
"""Class for XY scalar sources.
Parameters:
x (numpy.array): linear array with equidistant positions. The number of data is preferibly :math:`2^n` .
y (numpy.array): linear array wit equidistant positions for y values
wavelength (float): wavelength of the incident field
info (str): String with info about the simulation
Attributes:
self.x (numpy.array): linear array with equidistant positions. The number of data is preferibly :math:`2^n` .
self.y (numpy.array): linear array wit equidistant positions for y values
self.wavelength (float): wavelength of the incident field.
self.u (numpy.array): (x,z) complex field
self.info (str): String with info about the simulation
"""
def __init__(self, x=None, y=None, wavelength=None, info=""):
super().__init__(x, y, wavelength, info)
self.type = 'Scalar_source_XY'
[docs]
def plane_wave(self, A=1, theta=0 * degrees, phi=0 * degrees, z0=0 * um):
"""Plane wave. self.u = A * exp(1.j * k *
(self.X * sin(theta) * cos(phi) +
self.Y * sin(theta) * sin(phi) + z0 * cos(theta)))
According to https://en.wikipedia.org/wiki/Spherical_coordinate_system: physics (ISO 80000-2:2019 convention)
Parameters:
A (float): maximum amplitude
theta (float): angle in radians
phi (float): angle in radians
z0 (float): constant value for phase shift
"""
k = 2 * pi / self.wavelength
self.u = A * exp(1.j * k *
(self.X * sin(theta) * cos(phi) +
self.Y * sin(theta) * sin(phi) + z0 * cos(theta)))
[docs]
def gauss_beam(self,
r0,
w0,
z0,
alpha=0 * degrees,
beta=0 * degrees,
A=1,
theta=0. * degrees,
phi=0 * degrees):
"""Gauss Beam.
Parameters:
r0 (float, float): (x,y) position of center
w0 (float, float): (wx,wy) minimum beam width
z0 (float): (z0x, z0y) position of beam width for each axis (could be different)
alpha (float): rotation angle of the axis of the elliptical gaussian amplitude
beta (float): rotation angle of the axis of the main directions of the wavefront (it can be different from alpha)
A (float): maximum amplitude
theta (float): angle in radians (angle of k with respect to z))
References:
"""
if isinstance(w0, (float, int, complex)):
w0 = (w0, w0)
if isinstance(z0, (float, int, complex)):
z0 = (z0, z0)
w0x, w0y = w0
# w0 = sqrt(w0x * w0y)
x0, y0 = r0
z0x, z0y = z0
k = 2 * np.pi / self.wavelength
# only for x axis.
z_rayleigh_x = k * w0x**2 / 2
z_rayleigh_y = k * w0y**2 / 2
phaseGouy_x = np.arctan2(z0x, z_rayleigh_x)
phaseGouy_y = np.arctan2(z0y, z_rayleigh_y)
wx = w0x * np.sqrt(1 + (z0x / z_rayleigh_x)**2)
wy = w0y * np.sqrt(1 + (z0y / z_rayleigh_y)**2)
# w = sqrt(wx * wy)
if z0x == 0:
R_x = 1e10
else:
R_x = -z0x * (1 + (z_rayleigh_x / z0x)**2)
if z0y == 0:
R_y = 1e10
else:
R_y = -z0y * (1 + (z_rayleigh_y / z0y)**2)
amplitude = (A * (w0x / wx) * (w0y / wy) * np.exp(
-(self.X * np.cos(alpha) + self.Y * np.sin(alpha) - x0)**2 /
(wx**2)) * np.exp(
-(-self.X * np.sin(alpha) + self.Y * np.cos(alpha) - y0)**2 /
(wy**2)))
phase1 = np.exp(1.j * k * (self.X * np.sin(theta) * np.cos(phi) +
self.Y * np.sin(theta) * np.sin(phi)))
phase2 = np.exp(1j * (k * z0x - phaseGouy_x + k * ((self.X * np.cos(beta) + self.Y * np.sin(beta))**2) / (2 * R_x))) * \
np.exp(1j * (k * z0y - phaseGouy_y + k * ((-self.X * np.sin(beta) + self.Y * np.cos(beta))**2) / (2 * R_y)))
self.u = amplitude * phase1 * phase2
[docs]
def spherical_wave(self, r0, z0, A=1, radius=0, normalize=False):
"""Spherical wave.
Parameters:
A (float): maximum amplitude
r0 (float, float): (x,y) position of source
z0 (float) or (float, float): z0 or (zx0, zy0) position of source. When two parameters, the source is stigmatic
radius (float): radius for mask
normalize (bool): If True, maximum of field is 1
"""
k = 2 * np.pi / self.wavelength
x0, y0 = r0
if isinstance(z0, (int, float)):
z0 = (z0, z0)
z0x, z0y = z0
R2 = (self.X - x0)**2 + (self.Y - y0)**2
# Rz = sqrt((self.X - x0)**2 + (self.Y - y0)**2 + z0**2)
if z0x == 0:
R_x = 1e10
else:
R_x = np.sqrt((self.X - x0)**2 + z0x**2)
if z0y == 0:
R_y = 1e10
else:
R_y = np.sqrt((self.Y - y0)**2 + z0y**2)
if radius > 0:
amplitude = (R2 <= radius**2)
else:
amplitude = 1
self.u = amplitude * A * np.exp(-1.j * np.sign(z0x) * k * R_x)*np.exp(-1.j * np.sign(z0y) * k * R_y) / np.sqrt(R_x*R_y)
if normalize is True:
self.u = self.u / np.abs(self.u.max() + 1.012034e-12)
[docs]
def vortex_beam(self, A, r0, w0, m):
"""Vortex beam.
Parameters:
A (float): Amplitude
r0 (float, float): (x,y) position of source
w0 (float): width of the vortex beam
m (int): order of the vortex beam
Example:
vortex_beam(r0=(0 * um, 0 * um), w0=100 * um, m=1)
"""
if isinstance(w0, (float, int, complex)):
w0x, w0y = w0, w0
else:
w0x, w0y = w0
x0, y0 = r0
amplitude = ((self.X - x0) + 1.j * sign(m) *
(self.Y - y0))**np.abs(m) * np.exp(-(
(self.X - x0)**2 / w0x**2 +
(self.Y - y0)**2 / w0y**2))
self.u = A * amplitude / np.abs(amplitude).max()
[docs]
def hermite_gauss_beam(self, r0, A, w0, n, m, z, z0):
"""Hermite Gauss beam.
Parameters:
A (float): amplitude of the Hermite Gauss beam.
r0 (float, float): (x,y) position of the beam center.
w0 (float, float): Gaussian waist.
n (int): order in x.
m (int): order in y.
z (float): Propagation distance.
z0 (float, float): Beam waist position at each dimension
Example:
hermite_gauss_beam(A=1, r0=(0,0), w0=(100*um, 50*um), n=2, m=3, z=0)
"""
# Prepare space
X = self.X - r0[0]
Y = self.Y - r0[1]
r2 = sqrt(2)
if isinstance(w0, (float, int, complex)):
w0x, w0y = w0, w0
else:
w0x, w0y = w0
if isinstance(z0, (float, int, complex)):
z0x, z0y = z0, z0
else:
z0x, z0y = z0
k = 2 * pi / self.wavelength
# Calculate propagation
zx = z - z0x
zRx = k * w0x**2 / 2
wx = w0x * sqrt(1 + zx**2 / zRx**2)
if zx == 0:
Rx = np.inf
else:
Rx = zx + zRx**2 / zx
zy = z - z0y
zRy = k * w0y**2 / 2
wy = w0y * sqrt(1 + zy**2 / zRy**2)
if zy == 0:
Ry = np.inf
else:
Ry = zy + zRy**2 / zy
# Calculate amplitude
A = A * sqrt(2**(1 - n - m) /
(pi * factorial(n) * factorial(m))) * sqrt(w0x * w0y /
(wx * wy))
Ex = eval_hermite(n, r2 * X / wx) * exp(-X**2 / wx**2)
Ey = eval_hermite(m, r2 * Y / wy) * exp(-Y**2 / wy**2)
# Calculate phase
Ef = exp(1j * k * (X**2 / Rx + Y**2 / Ry)) * exp(
-1j * (0.5 + n) * np.arctan(zx / zRx)) * exp(
-1j * (0.5 + m) * np.arctan(zy / zRy)) * exp(1j * k *
(zx + zy) / 2)
self.u = A * Ex * Ey * Ef
[docs]
def laguerre_beam(self, r0, A, w0, n, l, z, z0):
"""Laguerre beam.
Parameters:
A (float): amplitude of the Hermite Gauss beam.
r0 (float, float): (x,y) position of the beam center.
w0 (float): Gaussian waist.
n (int): radial order.
l (int): angular order.
z (float): Propagation distance.
z0 (float): Beam waist position.
Example:
laguerre_beam(A=1, r0=(0 * um, 0 * um), w0=1 * um, p=0, l=0, z=0)
"""
# Prepare space
X = self.X - r0[0]
Y = self.Y - r0[1]
Ro2 = X**2 + Y**2
Ro = np.sqrt(Ro2)
Th = np.arctan2(Y, X)
# Parameters
r2 = sqrt(2)
z = z - z0
k = 2 * pi / self.wavelength
# Calculate propagation
zR = k * w0**2 / 2
w = w0 * sqrt(1 + z**2 / zR**2)
if z == 0:
R = np.inf
else:
R = z + zR**2 / z
# Calculate amplitude
A = A * w0 / w
Er = laguerre_polynomial_nk(2 * Ro2 / w**2, n, l) * exp(
-Ro2 / w**2) * (r2 * Ro / w)**l
# Calculate phase
Ef = exp(1j * (k * Ro2 / R + l * Th)) * \
exp(-1j * (1 + n) * np.arctan(z / zR))
self.u = A * Er * Ef
[docs]
def zernike_beam(self, A, r0, radius, n, m, c_nm):
"""Zernike beam.
Parameters:
A (float): amplitude of the Zernike beam beam
r0 (float, float): (x,y) position of source
radius (float): width of the beam
n (list): list of integers with orders
m (list): list of integers with orders
c_nm (list): list of coefficients, in radians
Example:
zernike_beam(A=1, r0=(0,0), radius=5 * mm, n=[1, 3, 3, 5, 5, 5], m=[1, 1, 3, 1, 3, 5], c_nm=[.25, 1, 1, 1, 1, 1])
"""
# normalizing to radius 1
x0, y0 = r0
R = sqrt((self.X - x0)**2 + (self.Y - y0)**2) / radius
# phase as sum of Zernike functions
phase = zeros(self.X.shape, dtype=float)
for s in range(len(n)):
phase = phase + c_nm[s] * fZernike(self.X - x0, self.Y - y0, n[s],
m[s], radius)
self.u = A * exp(1.j * np.real(phase))
[docs]
def bessel_beam(self,
A,
r0,
alpha,
n,
theta=0 * degrees,
phi=0 * degrees,
z0=0):
"""Bessel beam produced by an axicon. Bessel-beams are generated using 2D axicons.
Parameters:
A (float): amplitude of the Bessel beam
r0 (float, float): (x,y) position of source
alpha (float): angle of the beam generator
n (int): order of the beam
theta (float): angle in radians
phi (float): angle in radians
z0 (float): constant value for phase shift
References:
J. Durnin, J. Miceli, and J. H. Eberly, Phys. Rev. Lett. 58, 1499 (1987).
F. Courvoisier, et al. "Surface nanoprocessing with nondiffracting femtosecond Bessel beams" Optics Letters Vol. 34, No. 20 3163 (2009)
"""
k = 2 * np.pi / self.wavelength
x0, y0 = r0
R = sqrt((self.X - x0)**2 + (self.Y - y0)**2)
beta = k * np.cos(alpha)
if n == 0:
jbessel = j0(k * np.sin(alpha) * R)
elif n == 1:
jbessel = j1(k * np.sin(alpha) * R)
else:
jbessel = jv(n, k * np.sin(alpha) * R)
self.u = A * jbessel * np.exp(1j * beta * z0) * np.exp(
1.j * k *
(self.X * sin(theta) * cos(phi) + self.Y * sin(theta) * sin(phi)) +
z0 * cos(theta))
[docs]
def plane_waves_dict(self, params):
"""Several plane waves with parameters defined in dictionary
Parameters:
params: list with a dictionary:
A (float): maximum amplitude
theta (float): angle in radians
phi (float): angle in radians
z0 (float): constant value for phase shift
"""
# Definicion del vector de onda
k = 2 * pi / self.wavelength
self.u = np.zeros_like(self.u, dtype=complex)
for p in params:
self.u = self.u + p['A'] * exp(
1.j * k *
(self.X * sin(p['theta']) * cos(p['phi']) + self.Y *
sin(p['theta']) * sin(p['phi']) + p['z0'] * cos(p['theta'])))
[docs]
def plane_waves_several_inclined(self, A, num_beams, max_angle, z0=0):
"""Several paralel plane waves
Parameters:
A (float): maximum amplitude
num_beams (int, int): number of beams in the x and y directions
max_angle (float, float): maximum angle of the beams
z0 (float): position of the beams
"""
num_beams_x, num_beams_y = num_beams
max_angle_x, max_angle_y = max_angle
t = np.zeros_like(self.u, dtype=complex)
anglex = max_angle_x / num_beams_x
angley = max_angle_y / num_beams_y
for i in range(num_beams_x):
for j in range(num_beams_y):
theta = np.pi / 2 - max_angle_x / 2 + anglex * (i + 0.5)
phi = np.pi / 2 - max_angle_y / 2 + angley * (j + 0.5)
self.plane_wave(A, theta, phi, z0)
t = t + self.u
self.u = t
[docs]
def gauss_beams_several_parallel(self,
r0,
A,
num_beams,
w0,
z0,
r_range,
theta=0 * degrees,
phi=0 * degrees):
"""Several parallel gauss beams
Parameters:
A (float): maximum amplitude
num_beams (int, int): number of gaussian beams (equidistintan) in x and y direction.
w0 (float): beam width of the bemas
z0 (float): constant value for phase shift
r0 (float, float): central position of rays (x_c, y_c)
r_range (float, float): range of rays x, y
theta (float): angle
phi (float): angle
"""
x_range, y_range = r_range
x_central, y_central = r0
num_beams_x, num_beams_y = num_beams
t = np.zeros_like(self.u, dtype=complex)
dist_x = x_range / num_beams_x
dist_y = y_range / num_beams_y
for i in range(num_beams_x):
xi = x_central - x_range / 2 + dist_x * (i + 0.5)
for j in range(num_beams_y):
yi = y_central - y_range / 2 + dist_y * (j + 0.5)
self.gauss_beam(r0=(xi, yi),
w0=w0,
z0=z0,
A=A,
theta=theta,
phi=phi)
t = t + self.u
self.u = t
[docs]
def gauss_beams_several_inclined(self, A, num_beams, w0, r0, z0,
max_angle):
"""Several inclined gauss beams
Parameters:
A (float): maximum amplitude
num_beams (int, int): number of gaussian beams (equidistintan) in x and y direction.
w0 (float): beam width
r0 (float, float): central position of rays (x_c, y_c)
z0 (float): constant value for phase shift
max_angle (float, float): maximum angles
"""
num_beams_x, num_beams_y = num_beams
max_angle_x, max_angle_y = max_angle
t = np.zeros_like(self.u, dtype=complex)
angle_x = max_angle_x / num_beams_x
angle_y = max_angle_y / num_beams_y
for i in range(num_beams_x):
thetai = np.pi / 2 - max_angle_x / 2 + angle_x * (i + 0.5)
for j in range(num_beams_y):
phii = np.pi / 2 - max_angle_y / 2 + angle_y * (j + 0.5)
self.gauss_beam(r0=r0,
w0=w0,
z0=z0,
A=A,
theta=thetai,
phi=phii)
t = t + self.u