ESBO-ETC/esbo_etc/classes/psf/Zemax.py
2020-04-29 17:37:07 +02:00

119 lines
5.6 KiB
Python

import numpy as np
import re
from logging import warning
import astropy.units as u
from scipy.interpolate import RegularGridInterpolator
from scipy.integrate import nquad
from scipy.optimize import bisect
from .IPSF import IPSF
from typing import Union
from ...lib.helpers import error
class Zemax(IPSF):
"""
A class for modelling the PSF from a Zemax output file
"""
@u.quantity_input(wl="length", d_aperture="length")
def __init__(self, file: str, f_number: float, wl: u.Quantity, d_aperture: u.Quantity):
"""
Initialize a new PSF from a Zemax file.
Parameters
----------
file : str
Path to the Zemax-file. The origin of the coordinate system is in the lower left corner of the matrix
f_number : float
The working focal number of the optical system
wl : Quantity
The central wavelength which is used for calculating the PSF
d_aperture : Quantity
The diameter of the telescope's aperture.
"""
# Store parameters
self.__f_number = f_number
self.__wl = wl
self.__d_aperture = d_aperture
# Read PSF from file
with open(file, encoding="utf16") as fp:
self.__psf = np.genfromtxt((x.replace(",", ".") for x in fp), delimiter='\t', skip_header=21)
# Read header parameters from the file
with open(file, encoding="utf16") as fp:
head = [next(fp) for _ in range(21)]
# Parse shape of the grid and check the read PSF-array
shape = [int(x) for x in re.findall("[0-9]+", list(filter(re.compile("Image grid size: ").match, head))[0])]
if shape != list(self.__psf.shape):
warning("Not all PSF entries read.")
# Parse and calculate the grid width
grid_delta = [float(x.replace(",", ".")) for x in
re.findall("[0-9]+,*[0-9]*", list(filter(re.compile("Data area is ").match, head))[0])]
unit = re.findall(".+(?=\\.$)", re.sub("Data area is [0-9]+,*[0-9]* by [0-9]+,*[0-9]* ", "",
list(filter(re.compile("Data area is ").match, head))[0]))[0]
self.__grid_delta = np.array(grid_delta) / np.array(shape) << u.Unit(unit)
# Parse the center point of the PSF in the grid
self.__center_point = [int(x) for x in
re.findall("[0-9]+", list(filter(re.compile("Center point is: ").match, head))[0])]
def calcReducedObservationAngle(self, contained_energy: Union[str, int, float, u.Quantity]) -> u.Quantity:
"""
Calculate the reduced observation angle in lambda / d_ap for the given contained energy.
Parameters
----------
contained_energy : Union[str, int, float, u.Quantity]
The percentage of energy to be contained within a circle with the diameter reduced observation angle.
Returns
-------
reduced_observation_angle: Quantity
The reduced observation angle in lambda / d_ap
"""
if type(contained_energy) == str:
try:
contained_energy = float(contained_energy) / 100.0 * u.dimensionless_unscaled
except ValueError:
error("Could not convert encircled energy to float.")
elif type(contained_energy) in [int, float]:
contained_energy = contained_energy * u.dimensionless_unscaled
# Create an linear interpolation function for the PSF and the corresponding grid coordinates
x_range = np.arange(-(self.__center_point[0] - 1), self.__psf.shape[0] - self.__center_point[0] + 1)
y_range = np.arange(-(self.__center_point[1] - 1), self.__psf.shape[1] - self.__center_point[1] + 1)
interp = RegularGridInterpolator((y_range, x_range), np.flip(self.__psf, axis=0))
# Calculate the maximum possible radius as the smallest distance from the center of the PSF to the borders of
# the grid.
r_max = min(self.__center_point[0] - 1, self.__center_point[1] - 1,
self.__psf.shape[0] - self.__center_point[0], self.__psf.shape[1] - self.__center_point[1])
# Calculate the overall integral of the PSF
total = np.sum(self.__psf)
# Find the radius of the circle containing the given percentage of energy. Therefore, the interpolation
# function is numerically integrated within the radius. The Integration radius is optimized using bisection.
try:
r = bisect(lambda r_c: contained_energy.value -
nquad(lambda x, y: interp(np.array([y, x])),
[lambda y: [-1 * np.sqrt(r_c ** 2 - y ** 2), np.sqrt(r_c ** 2 - y ** 2)],
[-r_c, r_c]], opts={"epsrel": 1e-1})[0] / total, 0, r_max, xtol=0.1)
except ValueError:
r = r_max
# Calculate the reduced observation angle for the radius of the circle. Therefore, first convert the radius in
# grid elements to plate size, then calculate the corresponding observation angle and reduce it.
reduced_observation_angle = r * self.__grid_delta[0] / (self.__f_number * self.__d_aperture) * \
self.__d_aperture / self.__wl
return reduced_observation_angle * u.dimensionless_unscaled
def mapToGrid(self, grid: np.ndarray) -> np.ndarray:
"""
Map the integrated PSF values to a sensor grid.
Parameters
----------
grid : ndarray
The grid to map the values to. The values will only be mapped onto entries with the value 1.
Returns
-------
grid : ndarray
The grid with the mapped values.
"""
pass