"""Collection of tools for computing different maps."""
import numpy as np
from scipy.special import erf
[docs]def make_coordinates(N, a):
"""Makes the 1D and square 2D grid of coordinates for map calculations.
Parameters
----------
N_samples : int
Number of samples in each dimension (`N` x `N`)
a : float
Dimension of square grid (`a` x `a`) to generate
Returns
-------
x,y : array_like[float]
1D arrays for the x and y coordinates
X,Y,R : array_like[float]
2D Arrays of the X, Y, and radial R coordinates
"""
x = y = np.linspace(-a, a, N)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X ** 2 + Y ** 2)
return x, y, X, Y, R
[docs]def circular_aperture(x, y, R_ap, x_offset=0, y_offset=0):
"""Circular aperture map.
Parameters
----------
x, y : array
1D arrays describing uniform 2D grid to compute map
over in meters
R : float
Radius of aperture in meters
x_offset, y_offset : float, optional
Offset of aperture from origin
"""
X, Y = np.meshgrid(x, y)
R = np.sqrt((X - x_offset) ** 2 + (Y - y_offset) ** 2)
ap_map = np.ones_like(X)
ap_map[R > R_ap] = 0
return ap_map
[docs]def surface_point_absorber(
xs, ys, w, h, power_absorbed, alpha=0.55e-6, kappa=1.38, zero_min=False
):
"""Models the surface deformation from a small point absorber in a coating of a
mirror. It calcaultes the thermo-elastic deformation due to excess heat being
deposited in the mirror.
Parameters
----------
xs, ys : array
1D array for the x and y axis to calculate the distortion over
w : double
Area of the absorber size
h : double
Thickness of mirror
power_absorbed : double
Amount of power absorbed over the area w
alpha : double, optional
Thermo-elastic coefficient of material, default value for fused silica
kappa : double, optional
Thermal conductivity of the material, default value for fused silica
Returns
-------
Height map in meters
Notes
-----
Equation from:
A. Brooks, et.al "Point absorbers in Advanced LIGO," Appl. Opt. (2021)
"""
dx = xs[1] - xs[0]
dy = ys[1] - ys[0]
r = np.sqrt(np.add.outer(ys ** 2, xs ** 2))
c0 = -1 / 2 - np.log((h ** 2 + np.sqrt(w ** 2 + h ** 2)) / w)
out = np.zeros_like(r)
# most of the time the absorber is very small, much smaller than
# the discretisation of the map, so we can simplify this calculation
if w > dx or w > dy:
b1 = np.abs(r) <= w
b2 = np.abs(r) > w
out[b1] = -1 / 2 * (r[b1] / w) ** 2
out[b2] = c0 + np.log((h ** 2 + np.sqrt(r[b2] ** 2 + h ** 2)) / r[b2])
else:
out = c0 + np.log((h ** 2 + np.sqrt(r ** 2 + h ** 2)) / r)
if zero_min:
out -= np.min(out)
return alpha * power_absorbed / (2 * np.pi * kappa) * out
[docs]def overlap_piston_coefficient(x, y, z, weight_spot: float):
"""Computes the amount of weighted piston term there is in some 2D data like a
surface map. This is computed by evaluating a weighted Hermite polynomial overlap
integral in an efficient manner.
Parameters
----------
x, y : array_like
1D Array of x and y describing the 2D plane of `z`
z : array_like
2D optical path depth [metres]
weight_spot : float
Beam spot size to weight over
Returns
-------
piston : float
amount of piston term
"""
dx = (x[1] - x[0]) / weight_spot
dy = (y[1] - y[0]) / weight_spot
Wx = np.exp(-((x / weight_spot) ** 2))
Wy = np.exp(-((y / weight_spot) ** 2))
k00 = np.dot((z @ Wx), Wy) * dx / (np.sqrt(np.pi) * 2) * dy / (np.sqrt(np.pi) * 2)
return 4 * k00
[docs]def overlap_tilt_coefficients(x, y, z, weight_spot: float):
"""Computes the amount of yaw and pitch terms present in a map's displacement data.
This is computed by evaluating a weighted Hermite polynomial overlap integral in an
efficient manner.
Parameters
----------
x, y : array_like
1D Array of x and y describing the 2D plane of `z`
z : array_like
2D optical path depth [metres]
weight_spot : float
Beam spot size to weight over
Returns
-------
k10, k01 : complex
Complex-valued overlap tilt for the HG10 and HG01 mode
"""
weight_spot /= np.sqrt(2)
dx = (x[1] - x[0]) / weight_spot
dy = (y[1] - y[0]) / weight_spot
Wx = np.exp(-((x / weight_spot) ** 2))
Wy = np.exp(-((y / weight_spot) ** 2))
Hn = 2 * x / weight_spot * Wx
Hm = 2 * y / weight_spot * Wy
k10 = np.dot((z @ Hn), Wy) * dx / (np.sqrt(np.pi) * 2) * dy / (np.sqrt(np.pi) * 2)
k01 = np.dot((z @ Wx), Hm) * dx / (np.sqrt(np.pi) * 2) * dy / (np.sqrt(np.pi) * 2)
return (
# convert back to normal x units in displacement map
4 * k10 / weight_spot,
4 * k01 / weight_spot,
)
[docs]def overlap_curvature_coefficients(x, y, z, weight_spot: float):
"""Computes the amount of x and y curvature terms present in a map's displacement
data.
This is computed by evaluating a weighted Hermite polynomial overlap integral in an
efficient manner.
Parameters
----------
x, y : array_like
1D Array of x and y describing the 2D plane of `z`
z : array_like
2D optical path depth [metres]
weight_spot : float
Beam spot size to weight over
Returns
-------
k20, k02 : complex
Complex-valued overlap coefficients for the HG20 and HG02 mode
"""
# Normalisation constant
weight_spot /= np.sqrt(2)
dx = (x[1] - x[0]) / weight_spot
dy = (y[1] - y[0]) / weight_spot
# Spot size weightings
Wx = np.exp(-((x / weight_spot) ** 2))
Wy = np.exp(-((y / weight_spot) ** 2))
# 2nd order HG mode
Hn = ((2 * x / weight_spot) ** 2 - 2) * Wx
Hm = ((2 * y / weight_spot) ** 2 - 2) * Wy
# Compute the overlap integrals
k20 = np.dot((z @ Hn), Wy) * dx / (np.sqrt(np.pi) * 2 ** 3) * dy / (np.sqrt(np.pi))
k02 = np.dot((z @ Wx), Hm) * dx / (np.sqrt(np.pi) * 2 ** 3) * dy / (np.sqrt(np.pi))
return (
# convert back to normal x units in displacement map
4 * k20 / weight_spot / weight_spot,
4 * k02 / weight_spot / weight_spot,
)
[docs]def overlap_1D_curvature_coefficients(x, z, weight_spot: float):
"""Computes the amount of spot size weighted quadratic `x**2` term there is present
in some 1D data.
Parameters
----------
x : ndarray(dtype=float)
Sample points, ideally should be symmetric about 0.
z : ndarray(dtype=float)
function at sample points x
weight_spot : float
Spot size to use as weighting
Returns
-------
quadratic_term : double
Weighted quadratic term of z
Notes
-----
This function is essentially evaluating a weighted Hermite polynomial
overlap integral with `z(x)` to determine the linear term.
.. math::
\\int_{\\min(x)}^{\\max(x)} H_{2}(x) z(x) W(x) dx
Where the weighting function is :math:`W(x) = exp(-x**2)`.
"""
weight_spot /= np.sqrt(2)
dx = (x[1] - x[0]) / weight_spot
# exponential weighting of spot size
W = np.exp(-((x / weight_spot) ** 2))
# Second order Hermite, H_2(x)
Hn = ((2 * x / weight_spot) ** 2 - 2) * W
# normalisation constant as integral of H_0(x) over
# domain is not 1
norm = (
0.5 * np.sqrt(np.pi) * (erf(x.max() / weight_spot) - erf(x.min() / weight_spot))
)
k20 = (z @ Hn).sum() * dx / norm
return (
# convert back to normal x units in displacement map
k20
/ weight_spot
/ weight_spot
/ 2
)
[docs]def overlap_1D_tilt_coefficients(x, z, weight_spot: float):
"""Computes the amount of spot size weighted linear `x` term there is present in
some 1D data.
Parameters
----------
x : ndarray(dtype=float)
Sample points, ideally should be symmetric about 0.
z : ndarray(dtype=float)
function at sample points x
weight_spot : float
Spot size to use as weighting
Returns
-------
linear_term : double
Weighted linear term of z
Notes
-----
This function is essentially evaluating a weighted Hermite polynomial
overlap integral with `z(x)` to determine the linear term.
.. math::
\\int_{\\min(x)}^{\\max(x)} H_{1}(x) z(x) W(x) dx
Where the weighting function is :math:`W(x) = exp(-x**2)`.
"""
weight_spot /= np.sqrt(2)
dx = (x[1] - x[0]) / weight_spot
# exponential weighting of spot size
W = np.exp(-((x / weight_spot) ** 2))
# Second order Hermite, H_1(x)
Hn = 2 * x / weight_spot * W
# normalisation constant as integral of H_0(x) over
# domain is not 1
norm = np.sqrt(np.pi) * (erf(x.max() / weight_spot) - erf(x.min() / weight_spot))
k10 = (z @ Hn).sum() * dx / norm
return (
# convert back to normal x units in displacement map
2
* k10
/ weight_spot
)
[docs]def rms(x, y, z, weight_spot: float, xo: float = 0, yo: float = 0):
"""Computes the spot weight RMS over some 2D data, such as optical path depth.
Parameters
----------
x, y : array_like
1D Array of x and y describing the 2D plane of `z`
z : array_like
2D optical path depth [metres]
weight_spot : float
Beam spot size to weight over
xo, yo : float
Origin of the beam position
Returns
-------
rms : float
Root mean squared in units of whatever `z` is
Notes
-----
Based on Equation 4 in:
A. Brooks, et.al
Overview of Advanced LIGO adaptive optics
Appl. Opt. 55, 8256-8265 (2016)
"""
dx = x[1] - x[0]
dy = y[1] - y[0]
Wx = np.exp(-(((x - xo) / weight_spot) ** 2))
Wy = np.exp(-(((y - yo) / weight_spot) ** 2))
scaling = np.sqrt(2 / np.pi / weight_spot ** 2)
deltaW = z
deltaWbar = ((deltaW @ Wx) @ Wy) * dx * dy * scaling
return np.sqrt((((deltaW - deltaWbar) ** 2 @ Wx) @ Wy) * dx * dy * scaling)
[docs]class BinaryReaderEOFException(Exception):
def __init__(self):
pass
def __str__(self):
return "Not enough bytes in file to satisfy read request"
[docs]class BinaryReader:
"""MetroPro binary data format reader."""
# Map well-known type names into struct format characters.
typeNames = {
"int8": "b",
"uint8": "B",
"int16": "h",
"uint16": "H",
"int32": "i",
"uint32": "I",
"int64": "q",
"uint64": "Q",
"float": "f",
"double": "d",
"char": "s",
}
[docs] def __init__(self, fileName):
self.file = open(fileName, "rb")
[docs] def read(self, typeName, size=None):
"""Read a datatype from the binary file.
Parameters
----------
typename : str
See `BinaryReader.typeNames`.
size : int
Number of bytes to read
"""
import struct
typeFormat = BinaryReader.typeNames[typeName.lower()]
typeFormat = ">" + typeFormat
typeSize = struct.calcsize(typeFormat)
if size is None:
value = self.file.read(typeSize)
if typeSize != len(value):
raise BinaryReaderEOFException
unpacked = struct.unpack(typeFormat, value)[0]
else:
value = self.file.read(size * typeSize)
if size * typeSize != len(value):
raise BinaryReaderEOFException
unpacked = np.zeros(size)
for k in range(size):
i = k * typeSize
unpacked[k] = struct.unpack(typeFormat, value[i : i + typeSize])[0]
return unpacked
[docs] def seek(self, offset, refPos=0):
"""offset in bytes and refPos gives reference position, where 0 means origin of
the file, 1 uses current position and 2 uses the end of the file."""
self.file.seek(offset, refPos)
def __del__(self):
self.file.close()
[docs]def read_metropro_file(filename):
"""Reading the metroPro binary data files. Translated from Hiro Yamamoto's
'LoadMetroProData.m'.
Parameters
----------
filename
Name of metropro data file.
"""
f = BinaryReader(filename)
# Read header
hData = read_metropro_header(f)
if hData["format"] < 0:
print(
"Error: Format unknown to readMetroProData()\nfilename: {:s}".format(
filename
)
)
return 0
# Read phase map data
# Skipping header and intensity data
f.seek(hData["size"] + hData["intNBytes"])
# Reading data
dat = f.read("int32", size=hData["Nx"] * hData["Ny"])
# Marking unmeasured data as NaN
dat[dat >= hData["invalid"]] = np.nan
# Scale data to meters
dat = dat * hData["convFactor"]
# Reshaping into Nx * Ny matrix
dat = dat.reshape(hData["Ny"], hData["Nx"])
# Flipping up/down, i.e., change direction of y-axis.
dat = dat[::-1, :]
# Auxiliary data to return
# dxy = hData['cameraRes']
# Unnecessary
# x1 = dxy * np.arange( -(len(dat[0,:])-1)/2, (len(dat[0,:])-1)/2 + 1)
# y1 = dxy * np.arange( -(len(dat[:,0])-1)/2, (len(dat[:,1])-1)/2 + 1)
return dat, hData # , x1, y1