Theoretical background¶
In very simple terms, harmoni-pm provides a framework to simulate how objects in the sky (described as intensity functions in some arbitrary system of coordinates) are detected by a CCD with a given set of parameters and uncertainties describing the behavior of the instrument.
Transform theory¶
As light traverses the optical system, the relative position of the photons with respect to the optical axis changes not only due to changes in the focal length, but also due to aberrations, reflections, etc. For the sake of simplicity, we can assume geometrical optics and understand the effects of the optical elements as \(T: \mathbb R^2\to\mathbb R^2\) functions between pairs of coordinates.
Transforms are bi-directional, with the forward direction being that that transforms object plane coordinates into image plane coordinates, and the backward direction that that transforms image plane coordinates into object plane coordinates. This implies that transforms must be invertible. Although theoretically \(T_b=T_f^{-1}\) and therefore \(T_b\circ T_f=\mathbb I\), this is rarely achievable in the real world. In practical terms, it is required that the composition of the same transform in both backward and forward directions equals the identity transform up to certain (hopefully small) numerical error:
Integration and energy transport¶
We refer by integration to the calculation of the number of photons that arrive to each pixel in a given period of time, assuming an intensity pattern at the other end of the optical system. This calculation is implemented the ImageSampler component and is based on radiative transfer. Let \(dE_t\) the energy differential corresponding to the amount of radiative energy coming from a region in the sky of radiative intensity \(I\) and angular size \(d\Omega\) entering an aperture of area \(A\) in a time \(dt\) and a small range of frequencues \(d\nu\). This quantity can be written as:
\[dE_t = Id\Omega A\text{cos }\theta d\nu dt\]
Since the field is relatively small, \(\text{cos }\theta\approx1\) for all the rays that fall in the detector. On the other hand, the amount of energy \(dE_d\) arriving to a point in the detector at the instrument’s focal plane will be a fraction of the total energy that entered the telescope due to opacities, obstructions and aperture stops. This dependence can be written in terms of a dimensionless efficiency factor \(\beta\):
\[dE_d = \beta dE_t\]
TODO: consider \(e^{-\tau}\)!
Let \(F\) the flux corresponding to the region \(d\Omega\) that is focused in a small area \(dS\) in the focal plane. The energy differential in the detector will be:
\[dE_d = FdSd\nu dt\]
which allows us to deduce the following expression for the flux in the focal plane:
\[F = \beta IA \frac{d\Omega}{dS} = \beta\frac{IA}{f^2}\]
where \(f\) is the total focal length, if we assume aberration-free optics. If the power spectral density of the light is known, we can calculate the photon flux by:
\[Q = \int F \frac{d\nu}{h\nu}\]
which, multiplied by the pixel area \(A_p\) and assuming that the CCD Nyquist-samples the image plane, provides the number of photons per unit time striking a pixel. The true number of photons that arrived to a specific pixel in an integration time \(\Delta t\) can be modelled as a Poisson-distributed random variable with \(\lambda=QA_p\Delta t\).
Focal distance and aberrations¶
In the previous point, we stated that \(\frac{d\Omega}{dS}=\frac{1}{f^2}\) if we assume aberration-free optics. This will rarely be the case in the real world, and therefore a proper understanding of the relationship between aberrations and focal distance is necessary.
Let \(T_b(\vec{x})=[\phi(\vec{x}), \psi(\vec{x})]\) be the compound backward transform of all the optical elements in the optical path, where \(\vec{x}=x\hat{e}_x+y\hat{e}_y\) is the position vector of a point over the focal plane and \(\phi(\vec{x})\) and \(\psi(\vec{x})\) are two angular coordinates of the corresponding point in the sky with respect to the telescope pointing and centered in the equator. In this system \(\phi\) is the longitude and \(\psi\) its latitude. The orthonormal basis \(\hat{e}_x,\hat{e}_y\) can be assumed to be conveniently aligned to the rows and columns of pixels of the detector.
The solid angle differential \(d\Omega=\text{cos }\psi d\phi d\psi\) is related to the focal plane area differential \(dS=dxdy\):
Where \(J_b\) is the Jacobian matrix of the backward transform:
In the paraxial approximation, \(\text{cos }\psi\approx1\). \(J_b\) allows a polar decomposition of the form \(AR\), where \(A\) is a symmetric matrix and \(R\) a rotation matrix that can be connected to field rotation and the specific pick-off arm configuration. In the absence of rotations \(J_b=A\) and in the absence of aberrations we can choose \(\phi\) so that it only depends on \(x\) and \(\psi\) only depends on \(y\), which implies that \(\frac{\partial\phi}{\partial y} = \frac{\partial\psi}{\partial x}=0\) and \(A\) will be simplified down to a diagonal matrix:
which is the ideal case assumed earlier. As radial aberrations are introduced, \(A\) will lose diagonality and the determinant will change. This relationship suggests a way to estimate the focal distances from the coordinate transform and measure the degree of aberration in the optics:
with \(E\) a dimensionless quality figure that measures how much the aberrations are affecting the image in a given point \(\vec{x}\) in the focal plane, \(f\) the theoretical focal distance and \(\tilde{J}_b(\vec{x})\) a numerical calculation of the Jacobian of the transform in the surroundings of \(\vec{x}\).
Oversampling¶
Coordinate transforms must be calculated at least once for every pixel in the detector. If the detector Nyquist-samples the image plane, it should be enough to measure the flux at a fixed relative offset in every pixel. However, sub-Nyquist structure may exist in the image plane and oversampling would be needed to prevent aliasing / Moiré.
For a given oversampling value \(M\), ImageSampler will measure the intensity of the image plane in \(M^2\) points inside the pixel surface and calculate the average intensity of all of them. The offset of the sampling point \(0 \leq i, j < M\) from the bottom left corner of a pixel of width \(\Delta x\) and height \(\Delta y\) will be:
Implementation¶
From the implementation perspective, the generation of these subpixel coordinates is vectorized and performed in several steps:
Generation of the matrix of row indices. A square matrix \(r_{ij}=i\) is created as the tensor product of a
numpy.linspace(0, M - 1, M)by a column vector of ones of the same size.Generation of the x offsets. The vector \(x_k=[f(r_{ij})_k + \frac12]\delta x\) is created, with \(f\) the
flattenoperator (which returns a vector of \(M^2\) components such that \(f(r_{ij})_k=r_{k\text{ mod }M,\lfloor k/M\rfloor}\)) and \(\delta x = \Delta x/M\).Generation of the y offsets. The vector \(y_k=[f(r^T_{ij})_k + \frac12]\delta y\) is created, with \(\delta y = \Delta y/M\).
Generation of the coordinate offset vector. The matrix \(\mathtt{xy}_{k2}\) is created, such that \(\mathtt{xy}_{k0}=x_k\) and \(\mathtt{xy}_{k1}=y_k\).
The resulting \(\mathtt{xy}\) vector is then added to the coordinate of the bottom-left corner of each pixel to derive the full oversampling points.
Slicing and evaluation¶
The integration of the light arriving at \(h\times w\) pixels with an oversampling of \(M\) requires the calculation of \(h\times w\times M^2\) coordinate transforms. In addition to the computatational cost of this operation, the memory allocation required by the vectorization of these operations may be too big to fit in the computer’s memory. These two problems are addressed by slicing and parallelization.
Slicing refers to the process of selecting subsets of pixels in the detector from which integration operations are vectorized in numpy arrays. The maximum size of a slice is that of a contiguous square of side HARMONI_IMAGE_SAMPLER_SLICE_SIZE pixels, currently 128. The actual size of the slice may be smaller near the right and upper edges of the detector. The total number of coordinate transforms to perform can be calculated by multiplying the number of pixels in a slice \(K\) by the square of the oversampling \(M^2\).
Details on coordinate generation¶
The generation of coordinates is vectorized as well, and is performed in two stages. In a first stage (implemented in ImageSampler._integrate_slice), matrices of row indices \(r_{ij}\) and column indices \(c_{ij}\) are created using the same stratregy as described in the oversampling implementation, and flattened in order to obtain a \(\mathtt{ij}_{k2}\) matrix of \(K\) rows in which each row represents the row and column indices of the pixels to integrate.
In a second stage (implemented in ImageSampler._integrate_pixels, called by ImageSampler._integrate_slice), the full coordinate list including subsampling is generated from an initial pixel list \(\mathtt{ij}_{k2}\). This stage is divided again in the following steps:
Tiling of the oversampling offset matrix. A matrix \(\mathtt{p\_{xy}}\) consisting of concatenating \(\mathtt{xy}\) \(K\) times is generated, producing an oversampling offset matrix of \(KM^2\) rows.
Repetition of the pixel index matrix. Each individual row of the index matrix \(\mathtt{ij}_{k2}\) is repeated contiguosly \(M^2\) times, growing the matrix up to \(KM^2\) rows.
Generation of the sampling point matrix. The sampling coordinates are generated by adding the following matrix to \(\mathtt{p\_{xy}}\):
where \(\vec{x}_0\) is a position vector that encodes the physical displacement of the detector.
Evaluation¶
The evaluation of the intensity field involves a data reduction due to oversampling, as the resulting intensity vector must match the original pixel index matrix passed to ImageSampler._integrate_pixels. The intensity vector returned is evaluated by averaging the intensities belonging to the same pixel:
where \(T_b\) is the backward transform, \(\bar{I}_k\) the reduced intensity vector and \(I\) the intensity function of the image plane.
Parallelization¶
TODO: include details on thread parallelization
Time integration¶
TODO: include details on Poisson sampling, etc.