utils Package¶

This package contains functions used by other methods in the toolbox. Most of these functions are from version 1 of the toolbox with some minor modifications and bug fixes.

Warning

These functions are likely to move in future releases and are not very well documented.

Special functions
Coordinate transformations
Translations and rotations
Helper functions
Geometry functions
Unclassified
Polarizability calculation

Special functions ¶

ott.utils.sbesselh1(n, kr, varargin)¶

SBESSELH1 spherical hankel function hn(kr) of the first kind, hn(kr) = sqrt(pi/2kr) Hn+0.5(kr).

hn = SBESSELH1(n,z) calculates spherical hankel function of first kind.

[hn,dzhn] = SBESSELH1(n,z) additionally, calculates the derivatives of the appropriate Ricatti-Bessel function divided by z.

Coordinate transformations ¶

ott.utils.rtp2xyz(r, theta, phi)¶

RTP2XYZ coordinate transformation from spherical to cartesian: r radial distance [0, Inf) theta polar angle, measured from +z axis [0, pi] phi azimuthal angle, measured from +x towards +y axes [0, 2*pi)

[x,y,z] = RTP2XYZ(r,theta,phi) takes vectors or scalars, outputs the spherical coordinates as vectors/scalars of the same size.

[x,y,z] = RTP2XYZ(r) same as above but with the coordinate packed into the vector/matrix r = [ r theta phi ].

x = RTP2XYZ(…) same as above with the result packed into the vector/matrix x = [ x y z ].

ott.utils.rtpv2xyzv(rv, thetav, phiv, r, theta, phi)¶

RTPV2XYZV spherical to cartiesn vector field conversion

[xv,yv,zv,x,y,z] = RTPV2XYZV(rv,thetav,phiv,r,theta,phi)

[vec_cart,pos_cart] = rtpv2xyzv(vec,pos)

Inputs must be column vectors or Nx3 matrices.

Translations and rotations ¶

ott.utils.translate_z(nmax, z, varargin)¶

Calculates translation matrices for translation of VSWFs along z axis.

Usage

[A,B] = translate_z(nmax,z) calculates translation matrices. The matrices are use as:

\[ \begin{align}\begin{aligned}M' = A M + B N\\N' = B M + A N\end{aligned}\end{align} \]

[A,B,C] = ott.utils.translate_z(nmax,z) additionally, calculates C, the scalar SWF translation coefficients in 3d packed form.

Parameters

nmax (int) – Determines the number of multipole terms to include in the translation matrices (multipole order). Can be a single integer or two integers for the [row, column] nmax. If the row, column indices don’t match, A and B will not be square.
z (numeric) – Translation distance.

A and B are sparse matrices, since only m’ = m VSWFs couple

If z is a vector/matrix only A’s and B’s will be outputted. A and B will be cells of matrices the length of the number of elements in z. To save time only use unique values of z.

Time may be saved by taking the conjugate transpose instead of calculating translations in the positive or negative direction.

Optional named parameters

‘type’ (enum) – Type of translation matrix to generate.
‘method’ (enum) – Method to calculate translation matrix.

Translation matrix types

‘sbesselj’ regular to regular. Default. Used for most particle or beam translations.
‘sbesselh1’ outgoing to regular. Can be useful for doing multiple scattering calculations.
‘sbesselh2’ incoming to regular
‘sbesselh1farfield’ outgoing to regular far-field limit
‘sbesselh2farfield’ incoming to regular far-field limit

Methods

‘gumerov’ – Use Gumerov method (default, recommended)
‘videen’ – Use Videen method. (not recommended for large translations of the beam, unstable)

Example usage

The following example calculates the A and B translation matrices and applies them to a Gaussian beam. A procedure similar to this is done when calling Bsc.translateZ with a distance:

beam = ott.BscPmGauss('NA', 0.9, 'index_medium', 1.0, ...
    'polarisation', [1, 0], 'wavelength0', 1064e-9);

z = 1.0e-6 ./ beam.wavelength;
[A, B] = translate_z([beam.Nmax, beam.Nmax], z);

% Apply translation matrices
new_beam = beam.translate(A, B);

ott.utils.rotx(angle_deg, varargin)¶

Create a 3x3 rotation matrix for rotation about x axis

R = rotx(angle_deg) calculate the rotation matrix for rotations from the +z towards +y axis.

R = rotx([a1, a2, a3, …]) returns a 3xN matrix of rotation matrices for each angle in the input.

Optional named arguments:

usecell bool True to output as cell array instead of 3xN matrix.: Default: false. The cell array has the same shape as angle_deg.

Replacement/extension to Matlab rotx function provided in the Phased Array System Toolbox.

ott.utils.roty(angle_deg, varargin)¶

Create a 3x3 rotation matrix for rotation about y axis

R = roty(angle_deg) calculate the rotation matrix for rotations from the +z towards +x axis.

R = roty([a1, a2, a3, …]) returns a 3xN matrix of rotation matrices for each angle in the input.

Optional named arguments:

usecell bool True to output as cell array instead of 3xN matrix.: Default: false. The cell array has the same shape as angle_deg.

Replacement/extension to Matlab roty function provided in the Phased Array System Toolbox.

ott.utils.rotz(angle_deg, varargin)¶

Create a 3x3 rotation matrix for rotation about z axis

R = rotz(angle_deg) calculate the rotation matrix for rotations from the +x towards +y axis.

R = rotz([a1, a2, a3, …]) returns a 3xN matrix of rotation matrices for each angle in the input.

Optional named arguments:

usecell bool True to output as cell array instead of 3xN matrix.: Default: false. The cell array has the same shape as angle_deg.

Replacement/extension to Matlab rotz function provided in the Phased Array System Toolbox.

ott.utils.rotation_matrix(rot_axis, rot_angle)¶

ROTATION_MATRIX calculates rotation matrix using Euler-Rodrigues formula.

R = rotation_matrix( axis, angle ) calculates the rotation about axis by angle (in radians).

R = rotation_matrix( axis_angle ) calculates the rotation about vector axis_angle, the angle is specified as the length of the vector (in radians).

ott.utils.wigner_rotation_matrix(nmax, R)¶

WIGNER_ROTATION_MATRIX rotation matrix for rotation of spherical harmonics or T-matrices.

D = WIGNER_ROTATION_MATRIX(nmax,R) calculates the rotation matrix for the VSH given a 3x3 coordinate rotation matrix R. Usage: a’ = D a.

This method from Choi et al., J. Chem. Phys. 111: 8825-8831 (1999) Note change in notation - here, use x’ = Rx (x is row vector), a’ = Da (a is column vector) etc.

Helper functions ¶

ott.utils.matchsize(varargin)¶

Checks that all vector inputs have the same number of rows.

Usage

[A,B,…] = matchsize(A,B,…) checks inputs have same number of rows, and expands single-row inputs by repetition to match the input row number.

Parameters

A,B,… (numeric) – Numeric arrays whose number of rows are to be matched.

Example

The following example shows has two inputs, a scalar and a row vector. The scalar is duplicated to match the length of the row vector:

A = 5;
B = [1; 2; 3];

[A,B] = matchsize(A, B);
disp(A)  % -> [5; 5; 5]
disp(B)  % -> [1; 2; 3]

ott.utils.threewide(a)¶

THREEWIDE creates colum vector with input repeated in 3 columns: the function can take a column of row vector input, the output will be a matrix with three columns.

You might find this useful for multiplying a vector of scalars with a column vector of 3-vectors.

ott.utils.iseven(input)¶

ISEVEN determines if an integer is even Outputs a matrix of the same size as input with 1 for even and 0 for odd entries.

Warning: Plays up if the the integer is of the order 10^16

ott.utils.isodd(input)¶

ISODD determines if an integer is odd Outputs a matrix the same size as input with 1 for odd and 0 for even entries.

Warning: Plays up if the the integer is of the order 10^16

ott.utils.rotate_3x3tensor(ualpha, varargin)¶

Apply a set of rotations to a 3x3 tensor

alpha = rotate_3x3tensor(ualpha, R, …) applies the operation alpha = R*ualpha*inv(R) if R is a 3x3N matrix of rotation matrices.

alpha = rotate_3x3tensor(ualpha, ‘direction’, dir, …) computes the appropriate rotation rotation matrix to rotate from the z-axis to the 3xN matrix of directions. This is useful for uniaxial materials. ‘direction’ is the [x; y; z] Cartesian coordinate. ‘sphdirection’ is the [phi; theta] Spherical coordinate.

Optional named parameters:

‘inverse’ bool When true, returns the inverse polarisability.: Default: false.

ott.utils.ka2nmax(ka)¶

Finds a reasonable Nmax to truncate at for given size parameter

Usage: Nmax = ka2nmax(ka) calculates reasonable maximum order, Nmax, to truncate beam beam coefficients/T-matrix at for a given size parameter.

Returns \(Nmax = |ka| + 3 (|ka|)^(1/3)\)

ott.utils.nmax2ka(Nmax)¶

NMAX2KA finds size parameter ka corresponding to Nmax

ka = NMAX2KA(Nmax) finds size parameter for maximum order, Nmax, which spherical expansions are truncated.

Truncation order is given by Nmax = ka + 3 (ka)^(1/3)

Geometry functions ¶

ott.utils.angulargrid(ntheta, nphi, behaviour)¶

ANGULARGRID makes a angular grid of points over a sphere

[theta,phi] = ANGULARGRID(N) generates two column N^2-by-1 matrices with theta (polar) and phi (azimuthal) angle pairs for N discrete evenly spaced polar and azimuthal angles.

[theta,phi] = ANGULARGRID(ntheta, nphi) specifies the number of evenly spaced points to use in the theta and phi direction.

[theta,phi] = ANGULARGRID(…, behaviour) uses behaviour to control the output type:

behaviour | output¶

0 | column vectors of all points 1 | vectors of all theta and phi values 2 | ntheta x nphi matrix of all points

Note that the output data values are the same for behaviours 0 and 2; they’re just arranged differently. To convert from one format to another: 2 -> 0: theta = theta(:); phi = phi(:); 0 -> 2: theta = reshape(theta,ntheta,nphi);

phi = reshape(phi,ntheta,nphi);

ott.utils.perpcomponent(A, n)¶

PERPCOMPONENT finds perpendicular (and optionally) parallel components of a vector relative to a reference vector.

perp_component = PERPCOMPONENT(A,n) calculates perpendicular component of row vector A or Nx3 matrix A of vectors. n is the reference vector.

[perp_component,parallel_component] = PERPCOMPONENT(A,n) calculates perpendicular and parallel components.

ott.utils.inpolyhedron(varargin)¶

INPOLYHEDRON Tests if points are inside a 3D triangulated (faces/vertices) surface

IN = INPOLYHEDRON(FV,QPTS) tests if the query points (QPTS) are inside the patch/surface/polyhedron defined by FV (a structure with fields ‘vertices’ and ‘faces’). QPTS is an N-by-3 set of XYZ coordinates. IN is an N-by-1 logical vector which will be TRUE for each query point inside the surface. By convention, surface normals point OUT from the object (see FLIPNORMALS option below if to reverse this convention).

INPOLYHEDRON(FACES,VERTICES,…) takes faces/vertices separately, rather than in an FV structure.

IN = INPOLYHEDRON(…, X, Y, Z) voxelises a mask of 3D gridded query points rather than an N-by-3 array of points. X, Y, and Z coordinates of the grid supplied in XVEC, YVEC, and ZVEC respectively. IN will return as a 3D logical volume with SIZE(IN) = [LENGTH(YVEC) LENGTH(XVEC) LENGTH(ZVEC)], equivalent to syntax used by MESHGRID. INPOLYHEDRON handles this input faster and with a lower memory footprint than using MESHGRID to make full X, Y, Z query points matrices.

INPOLYHEDRON(…,’PropertyName’,VALUE,’PropertyName’,VALUE,…) tests query points using the following optional property values:

TOL - Tolerance on the tests for “inside” the surface. You can think of tol as the distance a point may possibly lie above/below the surface, and still be perceived as on the surface. Due to numerical rounding nothing can ever be done exactly here. Defaults to ZERO. Note that in the current implementation TOL only affects points lying above/below a surface triangle (in the Z-direction). Points coincident with a vertex in the XY plane are considered INside the surface. More formal rules can be implemented with input/feedback from users.

GRIDSIZE - Internally, INPOLYHEDRON uses a divide-and-conquer algorithm to split all faces into a chessboard-like grid of GRIDSIZE-by-GRIDSIZE regions. Performance will be a tradeoff between a small GRIDSIZE (few iterations, more data per iteration) and a large GRIDSIZE (many iterations of small data calculations). The sweet-spot has been experimentally determined (on a win64 system) to be correlated with the number of faces/vertices. You can overwrite this automatically computed choice by specifying a GRIDSIZE parameter.

FACENORMALS - By default, the normals to the FACE triangles are computed as the cross-product of the first two triangle edges. You may optionally specify face normals here if they have been pre-computed.

FLIPNORMALS - (Defaults FALSE). To match a wider convention, triangle face normals are presumed to point OUT from the object’s surface. If your surface normals are defined pointing IN, then you should set the FLIPNORMALS option to TRUE to use the reverse of this convention.

Example:
tmpvol = zeros(20,20,20); % Empty voxel volume tmpvol(5:15,8:12,8:12) = 1; % Turn some voxels on tmpvol(8:12,5:15,8:12) = 1; tmpvol(8:12,8:12,5:15) = 1; fv = isosurface(tmpvol, 0.99); % Create the patch object fv.faces = fliplr(fv.faces); % Ensure normals point OUT % Test SCATTERED query points pts = rand(200,3)*12 + 4; % Make some query points in = inpolyhedron(fv, pts); % Test which are inside the patch figure, hold on, view(3) % Display the result patch(fv,’FaceColor’,’g’,’FaceAlpha’,0.2) plot3(pts(in,1),pts(in,2),pts(in,3),’bo’,’MarkerFaceColor’,’b’) plot3(pts(~in,1),pts(~in,2),pts(~in,3),’ro’), axis image % Test STRUCTURED GRID of query points gridLocs = 3:2.1:19; [x,y,z] = meshgrid(gridLocs,gridLocs,gridLocs); in = inpolyhedron(fv, gridLocs,gridLocs,gridLocs); figure, hold on, view(3) % Display the result patch(fv,’FaceColor’,’g’,’FaceAlpha’,0.2) plot3(x(in), y(in), z(in),’bo’,’MarkerFaceColor’,’b’) plot3(x(~in),y(~in),z(~in),’ro’), axis image

See also: UNIFYMESHNORMALS (on the <a href=”http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=43013”>file exchange</a>)

This is a third party function, not part of OTTv2, see tplicenses/stl_Sven.txt for information about licensing.

Unclassified ¶

Todo

These functions should be moved to other categories

ott.utils.paraxial_transformation_matrix(paraxial_order, basis_in, basis_out, normal_mode)¶

PARAXIAL_TRANSFORMATION_MATRIX produces paraxial beam mode conversion in a particular order.

[modeweights, col_modes, row_modes] = …: PARAXIAL_TRANSFORMATION_MATRIX(degree, basis_in, basis_out) or
[modeweights, col_modes, row_modes] = …: PARAXIAL_TRANSFORMATION_MATRIX( degree, basis_in, basis_out, normal_mode)

inputs:

degree : paraxial degree of modes e.g. gaussian is 0. basis_in : 0 vortex LG, 1 vortex HG, [2,xi] vortex IG. basis_out : 0 LG, 1 HG, [2,xi] IG. normal_mode : 0 is default vortex->non-vortex. (because of toolbox modes)

1 makes the conversion non-vortex->non-vortex.

outputs:

modeweightsweights of the conversion basis_in->basis_out. Format: each: row is the corresponding mode of the output basis, each column for the input basis, such that conj. tranpose(basis_1 –> basis_2) = basis_2 –> basis_1. holds for (vortex->non-vortex)’ == non-vortex->vortex.
col_modesoutputs the LG, HG or IG indices corresponding to each: COLUMN of the matrix. [p,l], [m,n], [o,m,p].
row_modesoutputs the LG, HG or IG indices corresponding to each ROW: of the matrix. [p,l], [m,n], [o,m,p].

ott.utils.paraxial_beam_waist(paraxial_order)¶

ott.utils.lgmode(p, l, r, phi, z, theta)¶

LGMODE calculates LG mode amplitude at z = 0

A = LGMODE(p,l,r,phi) calculates the LG mode amplitude for mode [p,l] at locations given in polar coordinates [r, phi]. r is in units of the beam width; r and phi can be matrices of equal size.

A = LGMODE(p,l,r,phi,z) computes the modes with z as well.

A = LGMODE(p,l,r,phi,z,theta) scales the beam waist according to the beam convergence angle theta (in degrees): w0=1/tan(theta).

ott.utils.legendrerow(n, theta)¶

LEGENDREROW gives the spherical coordinate recursion in m

pnm = LEGENDREROW(n, theta) gives the spherical recursion for a given n, theta.

This provides approximately no benefit over the MATLAB implimentation. It may provide a benefit in Octave. Inspiration from [Holmes and Featherstone, 2002] and [Jekeli et al., 2007].

ott.utils.laguerre(p, l, X)¶

LAGUERRE associated Laguerre function

L = LAGUERRE(p,l,X) evaluate associated Laguerre function. p and l must be integer scalars greater than zero

ott.utils.incecoefficients(p, xi)¶

ott.utils.hgmode(m, n, x, y, z, theta)¶

ott.utils.hermite(n, X)¶

ott.utils.emField(krtp, type, nm, ab, varargin)¶

EMFIELD calculates field from the vector spherical wave functions

[E, H, data] = emField(krtp, type, nm, ab, …) calculates the E and H field for unit-less spherical coordinates krtp. krtp is a Nx3 matrix of [radial, polar, azimuthal] coordinates. type must be ‘incoming’, ‘regular’, or ‘outgoing’.

Optional named parameters:

‘saveData’ bool saves data that can be used for repeated: calculations of the fields at these locations (default: nargout==3).

‘data’ data data to use from previous calculation ‘calcE’ bool calculate E field (default: true) ‘calcH’ bool calculate H field (default: true)

If internal fields are calculated only the theta and phi components of E are continuous at the boundary. Conversely, only the kr component of D is continuous at the boundary.

ott.utils.combined_index(in1, in2)¶

COMBINED_INDEX translates between (n,m) and combined index Mode indices and combined index are related by: ci = n * (n+1) + m.

[n,m] = COMBINED_INDEX(ci) calculates (n,m) from the combined index.

ci = COMBINED_INDEX(n,m) calculates the combined index from mode indices.

length = COMBINED_INDEX(Nmax, Nmax) calculates length of the beam vectors.

Nmax = COMBINED_INDEX(length) calculates Nmax from length of beam vectors.

Polarizability calculation ¶

ott.utils.polarizability.FCD(spacing, index, varargin)¶

Filtered coupled dipole polarizability

Usage

alpha = FCD(spacing, index) Calculates a Nx1 element vector containing the isotropic polarizabilities for N dipoles.

Parameters

spacing (numeric scalar) – lattice spacing parameter
index (Nx1 numeric) – Relative refractive indices for N dipoles.

Optional named arguments

k0 (numeric) – Wavenumber to scale spacing by. Default: 2*pi.

ott.utils.polarizability.LDR(spacing, index, varargin)¶

Lattice dispersion relation polarizablity

Polarizability calculation based on

Draine & Goodman, Beyond Clausius-Mossoti: wave propagation on a polarizable point lattice and the discrete dipole approximation, The Astrophysical Journal, 405:685-697, 1993 March 10

Usage

alpha = LDR(spacing, index, …) Calculates a Nx1 element vector containing the isotropic polarisabilities for N dipoles.

alpha = LDR(spacing, index, kvec, E0, …) As above but specifies the polarisability information for use with plane wave illumination.

Parameters

spacing (numeric scalar) – lattice spacing parameter
index (Nx1 numeric) – Relative refractive indices for N dipoles.
kvec (1x3 numeric) – Wave vector [kx, ky, kz]
E0 (1x3 numeric) – E-field polarisation [Ex, Ey, Ez]

Optional named arguments

k0 (numeric) – Wavenumber to scale spacing by. Default: 2*pi.

ott.utils.polarizability.CM(spacing, index)¶

Clausius-Mossoti Polarizability

Usage

alpha = CM(spacing, index) Calculates a Nx1 element vector containing the isotropic polarisabilities for N dipoles.

Parameters

spacing (numeric scalar) – lattice spacing parameter
index (Nx1 numeric) – Relative refractive indices for N dipoles.

utils Package¶

Special functions¶

Coordinate transformations¶

Translations and rotations¶

Helper functions¶

Geometry functions¶

behaviour | output¶

Unclassified¶

Polarizability calculation¶

Special functions ¶

Coordinate transformations ¶

Translations and rotations ¶

Helper functions ¶

Geometry functions ¶

Unclassified ¶

Polarizability calculation ¶