Any light field can be expressed as a sum of plane waves. Such an ensemble of plane waves is called the plane wave spectrum of the light field. The plane wave spectrum is the Fourier transform of the light field in the real space representation.
Since, this is such a basic technique, I wonder whether there is a standard tool, preferably open-source, with a simple interface via Python, Matlab or similar, to do the following (numerically):
- Start with a light field in real space, say a monochromatic *-Gaussian beam.
- compute the Fourier transform
- apply some operation that changes the plane wave spectrum, say a boundary between two media
- compute the inverse Fourier transfrom
- produce some output, e.g. a plot
Computationally, all it takes is a library to perform the FFT. I'm looking for a framework that would wrap this along with some phyiscal concepts.
Answer
The reason that these sorts of libraries don't exist is because the particular algorithm that you use to do the calculation will depend upon the exact details of the light field and the input and output planes you are trying to compute.
For example, let me outline the the simplest case for this sort of calcualtion:
- The input light field $g_0$ has a slowly varying phase profile
- You want to propagate a short distance relative to the size of the beam (or, in a more pure-mathematical approach, your input field has infinite extent)
- The input and output planes are parallel and differ only by a translation along the z-axis
- The input and output computation grids have the same spacing
- This grid spacing is sufficiently small to nyquist sample the (spatially band-limited) input field.
In this case, you can easily use the angular spectrum propagation (ASP) method with FFTs to quickly compute your propagated wavefront. You seem to understand basic scalar diffraction, but to clarify for other readers, this means you compute:
$$\begin{eqnarray} G_0(\xi,\eta) = \mathcal{F}\left[ g_0(x,y) \right] \\ G_z(\xi,\eta) = A_z(\xi,\eta)G_0(\xi,\eta) \\ g_z(x,y) = \mathcal{F}^{-1}\left[ G_z(\xi,\eta) \right] \end{eqnarray}$$
Where $\mathcal{F}$ is the fourier transform from spatial coordinates $(x,y)$ to spatial frequency coordinates $(\xi,\eta)$, and $A_z(\xi,\eta)$ is the propagation kernel for distance $z$. In the future I will simply let capital letters denote Fourier transform of corresponding lower case letters. We typically define for ASP:
$$ A_z(\xi,\eta) = \exp \left[ i 2 \pi \frac{z}{\lambda} \sqrt{1- (\lambda \xi)^2 - (\lambda \eta)^2} \right]$$
where $\lambda$ is the wavelength.
However, this quickly falls apart if any of the conditions I outlined above are not satisfied. If the input field varies too rapidly (or conversely if the computation grid isn't fine enough to Nyquist sample those variations) then you will have severe aliasing in the Fourier transforms, and your result will be corrupted. If the input and output planes aren't parallel, you simply cannot use this method as it is formulated. If you need the output grid to be sampled differently than the input grid, then an FFT algorithm will no longer work, because the essential trade-off that you make for the speed of an FFT is that you cannot select arbitrary points to compute in Fourier space.
There are other propagation techniques of course. Lets say you have a converging beam, such as that immediately after a focusing lens. Barring an extremely fine computation grid spacing, or a very long focal length lens, ASP will fail because $G_0$ will be badly aliased. This is the textbook example of the applicability of Fresnel diffraction:
$$ g_z(x,y) = \frac{e^{i k z}}{i \lambda z} e^{i \frac{\pi}{\lambda z} (x^2 + y^2)} \mathcal{F} \left[ g_0(u,v) e^{i \frac{\pi}{\lambda z} (u^2 + v^2)} \right] $$
where $(u,v)$ are spatial coordinates in the input plane.
If we consider the input field $g_0$ to be composed of the product of a quadratic phase front $q(x,y) = \exp \left[-i \frac{\pi}{\lambda \mathcal{f}} (x^2 + y^2) \right]$ (representing the phase of the converging lens of focal length $\mathcal{f}$) and a residual flat component $s_0(x,y)$, we can see that over some range of propagation distances, the quadratic phase term inside the Fourier transform will act to "flatten" the transformed function such that there is no aliasing. Indeed, when $\mathcal{f} = z$ the transformed field will equal to the roughly flat $s_0$! In this case the resulting computation is basically identical to yet another propagation technique, Fraunhofer Diffraction, neglecting the minor detail of the quadratic phase term outside of the Fourier transform in Fresnel diffraction, and a scaling from angular coordinates to spatial coordinates.
Even with Fresnel/Fraunhoffer diffraction, we still run into problems if the propagation distance is too small. Note that the term $e^{i \frac{\pi}{\lambda z} (u^2 + v^2)}$ will have very high curvature when $z$ is small, once again producing aliasing problems in the Fourier transform. Typically, the solution to this is to propagate the field twice; first forwards by a distance $z_1$, then backwards by $-z_2$ such that the overall propagation distance is $z=z_1-z_2$, but realizing that this is necessary, and that this (or any other) particular solution is applicable is non-trivial.
Now what if your input field is significantly tilted (i.e. it isn't propagating parallel to the z-axis)? You will want your computation grid to include the output field, but most of the energy may be located far from the origin of the $(x,y)$ coordinate system. You could make your grid large enough to accommodate this, but ideally you'd compute an off-axis chunk of the output plane. Once again you need a new algorithm!
Consider that in general some people may have a need to propagate complicated optical fields to non-parallel planes or curved surfaces, to change the grid spacing in the output plane, or whatever other complication you can image, and you will begin to see how much decision making and complexity can be involved in selecting a propagation algorithm. Doing this in an automated way isn't necessarily impossible, but it is very difficult, and there is simply nobody who has undertaken the task of writing an open source library for this. I suspect such code exists in proprietary software packages, like CodeV or other lens design software, but remember that for a library to be developed there has to be somebody who is simultaneously interested (or paid sufficiently), qualified, and has the time to produce such a thing. Also remember that the audience for it would be limited to optical engineers and researchers, so there is little reward.
No comments:
Post a Comment