New docs

README.md

@@ -1,91 +1,64 @@
 # Axiprop
 A simple tool to compute optical propagation, based on the discrete Hankel and 
-Fourier transforms
+Fourier transforms.
 
-### Contents
+## Contents
 
 This library contains methods and convenience tools to model propagation of the 3D optical
-field. Computations can be done using a number of backends:
-- [NumPy](https://numpy.org) (**CPU**) optionally enhanced via [mkl_fft](https://github.com/IntelPython/mkl_fft) or
-[pyfftw](https://github.com/pyFFTW/pyFFTW)
-- [CuPy](https://cupy.dev) for **GPU** calculations via Nvidia CUDA API
-- [ArrayFire](https://arrayfire.com) for **GPU** calculations via CUDA or OpenCL APIs
-- [PyOpenCL](https://documen.tician.de/pyopencl) for **GPU** calculations via OpenCL API
+field. 
 
-Currently methods include:
+### Basic propagation methods include:
 - `PropagatorSymmetric`: cylindical axisymmetric propagator with the symmetric DHT proposed in 
 [[M. Guizar-Sicairos, J.C. Gutiérrez-Vega, JOSAA 21, 53 (2004)](https://doi.org/10.1364/JOSAA.21.000053)]
 - `PropagatorResampling`: cylindical axisymmetric propagator with a more generic DHT which allows for arbitrary
 sampling of radial axis [[K. Oubrerie, I.A. Andriyash et al, J. Opt. 24, 045503 (2022)](https://doi.org/10.1088/2040-8986/ac57d2)]
+- `PropagatorResamplingFresnel`: cylindical axisymmetric propagator based on the Fresnel approximation as given by `Eq. (4-17)` of [JW Goodman _Introduction to Fourier Optics_] and is suited for translations between far and near fields.
 - `PropagatorFFT2`: fully 3D FFT-based propagator
+- `PropagatorFFT2Fresnel`: fully 3D FFT-based propagator based on the Fresnel approximation 
 
-### Usage
-
-Consider a laser,
-```python 
-k0 = 2 * np.pi / 0.8e-6            # 800 nm wavelength
-tau = 35e-15/ (2*np.log(2))**0.5   # 35 fs FWHM duration
-R_las = 10e-3                      # 10 mm radius
-
-def fu_laser(kz, r):
-    """
-    Gaussian spot with the Gaussian temporal profile
-    """
-    profile_r = np.exp( -(r/R_las)**2 )
-    profile_kz = np.exp( -( (kz-k0) * c * tau / 2 )**2 )
-    return profile_r * profile_kz
-```
-
-and some focusing optics,
-```python
-f_N = 40                      # f-number f/40 
-f0 = 2 * R_las * f_N          # focal length
-```
+### Propagation with plasma models:
+- `Simulation` class includes a few algorithms (`'Euler'` `'Ralston'`, `'MP'`, `'RK4'`), method to provide adjustive steps, diagnostics
+- Plasma models: `PlasmaSimple` (linear current, $`n_{pe}(z)`$ ), `PlasmaSimpleNonuniform` (linear current, $`n_{pe}(z, r)`$), `PlasmaRelativistic` (non-linear current, $`n_{pe}(z, r)`$), `PlasmaIonization` (non-linear current and ionization, $`n_{g}(z, r)`$)
 
-Define the propagator,
-```python
-prop = PropagatorSymmetric((Rmax, Nr), (k0, L_kz, Nkz), Nr_end)
-```
-and setup the laser reflected from the focusing mirror
-```python
-A0 = laser_from_fu( fu_laser, prop.kz, prop.r )
-A0 = A0 * mirror_parabolic( f0, prop.kz, prop.r )
-```
+### Convenience tools include:
+- Container classes to create and handle time-frequency transformations for the **carrier-frequency-resolved** and **enveloped** fields.
+- Methods to convert fields between temporal and spatial representations
 
-Use `AXIPROP` to compute the field after propagation of `dz` distance 
-(e.g. `dz=f0` for field at focus):
-```python
-A0 = prop.step(A0, f0)
-```
-or evaluate it at `Nsteps` along some `Distance` around the focus,
-```python
-dz =  Distance / Nsteps
-zsteps = Nsteps * [dz,]
-zsteps[0] = f0 - Distance/2
-A_multi = prop.steps(A0, zsteps)
-```
+Computations can be done using a number of backends:
+- [NumPy](https://numpy.org) (**CPU**) optionally enhanced via [mkl_fft](https://github.com/IntelPython/mkl_fft) or
+[pyfftw](https://github.com/pyFFTW/pyFFTW)
+- [CuPy](https://cupy.dev) for **GPU** calculations via Nvidia CUDA to be discontinuedI
+- [PyOpenCL](https://documen.tician.de/pyopencl) for **GPU** calculations/OpenCL
+- [ArrayFire](https://arrayfire.com) for **GPU** calculations via CUDA/OpenCL (to be discontinued)
 
-Plot the results using your favorite tools 
+## Documentation
 
-![example_image](https://github.com/hightower8083/axiprop/blob/main/examples/example_figure.png)
+Documentation is organized in a few sections:
 
-For more info checkout the example notebooks for [radial](https://github.com/hightower8083/axiprop/blob/main/examples/example.ipynb) and [cartesian](https://github.com/hightower8083/axiprop/blob/main/examples/test2d.ipynb) cases, and also look for methods documentation.
+ - [Main equations for optical field](https://github.com/hightower8083/axiprop/blob/new-docs/docs/main_equations.md)
+ - [Representation of optical field](https://github.com/hightower8083/axiprop/blob/new-docs/docs/field_representation.md)
+ - [Propagation of optical field](https://github.com/hightower8083/axiprop/blob/new-docs/docs/propagators_vacuum.md)
 
-### Installation
+## Installation
+### From PyPI
+Major versions are released at [PyPI](https://pypi.org) and can be installed with `pip`:
+```bash
+pip install axiprop
+```
 
-Install `axiprop` by cloning the source 
+### From source
+You can build and install `axiprop` by cloning the source 
 ```bash
 git clone https://github.com/hightower8083/axiprop.git
 cd axiprop
 python setup.py install
 ```
-or directly via PiPy
+or directly from Github  with `pip`:
 ```bash
 pip install git+https://github.com/hightower8083/axiprop.git
 ```
 
-#### Additional requirements
-
+### Additional dependencies
 Note that, while base backend `NP` requires only NumPy and SciPy, other 
 backends have specific dependencies:
 - `NP_MKL`:  [mkl_fft](https://github.com/IntelPython/mkl_fft)

docs/field_representation.md

@@ -0,0 +1,25 @@
+### Field representation
+
+#### Temporal representation
+
+We assumed the field propagation along z-axis is  and the temporal representation. The later means that the full field $E(t,x,y,z)$ is considered at a fixed $z=z_0$, and is measured on a finite $(x,y)$-plane over a time interval $t\in[t_{min},t_{max}]$, so the input as $E_0 = E(t,x,y,z=z_0)$. 
+
+The field equation can be written as following:
+
+$$ \partial_z^2 E =  \frac{1}{c^2}  \partial_t^2 E - \nabla_{\perp}^2 E + \mu_0  \partial_t J   $$
+
+#### Geometries
+
+Two main geometries are usually considered
+- cartesian $(x, y, t)$, where transverse Laplace operator is $\nabla_\perp^2 = \partial_x^2 + \partial_y^2$
+- cylindrical $(r, \theta, t)$, where $\nabla_{\perp}^2 = r^{-1} \partial_r (r \partial_r) +  r^{-2} \partial_\theta^2$. 
+
+#### Envelope
+
+If one considers the field around its central frequency $\omega_0$ as:
+
+$$ E(t) = \mathrm{Re}[\hat{E}(t) \exp(- i \omega_0 t) ], $$
+
+in many practical cases the complex function $\hat{E}(t)$ can be assumed to be much slower than the actual $E(t)$. This presentation called envelope is often preferrable for analysis, and is also used in `Lasy`
+
+## [Return to README](https://github.com/hightower8083/axiprop/blob/new-docs/README.md#documentation)

docs/main_equations.md

@@ -0,0 +1,152 @@
+# Main considerations
+
+In the optical propagation calculations we describe evolution of electromagnetic field from a state defined at some time and space region to another state, that this field takes after some time and in a different space region.
+
+This problem can be described by the _Maxwell equations_:
+
+```math
+\begin{aligned}
+ & \nabla \times \mathbf{E} = - \partial_t \mathbf{B}\,,\\
+ & \nabla \times \mathbf{H} = \partial_t \mathbf{D} + \mathbf{J}\,,\\
+ & \nabla \mathbf{D} = \rho\,,\\
+ & \nabla  \mathbf{B} =0\,,
+\end{aligned}
+```
+
+accompanied by the _constitutive equations_
+
+```math
+\begin{aligned}
+ & \mathbf{D}  = \epsilon_0 \mathbf{E} + \mathbf{P} \,,\\
+ & \mathbf{H}  = \cfrac{1}{\mu_0} \mathbf{B} - \mathbf{M} \,.
+\end{aligned}
+```
+
+In the following we will assume that all charges of the system are free, i.e. $\epsilon_0$ and $\mu_0$ are the electric permittivity and magnetic permeability of vacuum. Moreover, at the current stage of development we will only consider **plasma** processes and discard the polarization and magentization of the media, thus setting $\mathbf{P}=0$ and $\mathbf{M}=0$.
+
+# Full equation
+
+In the optical field the electric and magnetic components are consistent with each other, and in most partical situations it is enough to follow the electric component, fot which the equation can be derived as:
+
+```math
+\begin{aligned}
+ & \nabla^2 \mathbf{E} - \frac{1}{c^2} \; \partial_t^2 \mathbf{E} = \mu_0\; \partial_t\mathbf{J}  + \frac{1}{\epsilon_0}\nabla \rho\,.
+\end{aligned}
+```
+
+Here the left hand side is the wave equation which describes propagation of the initial field in vacuum (the _optical beam_ whose dynamics we study). The terms in the right hand side are the currents and density modulations, which act as the _sources_ for the field generated by the media as a response to the optical field.
+
+#### Source terms
+
+The two source terms are essentially different in nature. The first one is excited _instantaneously_ as electrons generate electric currents being moved by the field. The second term constitutes the charge density modulations, and it presents a slower response and is typically considered to be much weaker. Let us demonstrate this by assuming the charge continuity in plasma and re-writing the source terms in the intergral form:
+
+```math
+\begin{aligned}
+ & \mu_0 \; \partial_t\mathbf{J}  + \frac{1}{\epsilon_0}\nabla \rho =  \frac{1}{\epsilon_0 } \int_{-\infty}^{t} dt' \left(\frac{1}{c^2} \; \partial_t^2 \mathbf{J} -   \nabla (\nabla \mathbf{J}) \right)\,.
+\end{aligned}
+```
+
+In the Fourier space, the terms under the integral appear as $k^2 \hat{\mathbf{J}}$ and $\mathbf{k}\cdot (\mathbf{k} \cdot \hat{\mathbf{J}})$, where $k=\omega/c$, and for the polarized field, it is enough to consider a single transverse component (e.g. $\mathbf{J} \parallel \mathbf{E} \parallel \mathbf{e}_x$) and discard others. The ratio between source terms therefore reads, $k^2/k_x^2$, and for the field wavelength $\lambda_0$, and transverse size $R$, it scales as $\propto (R/\lambda_0)^2$. It is easy to see, that in most practical situations this indeed allows to assume $\mu_0 \; \partial_t\mathbf{J} \gg \frac{1}{\epsilon_0}\nabla \rho$, and discard contribution of the charge density modulations compared to the induced currents.
+
+# Working equation
+
+To conclude this part, in the calculations we will consider the equation for the electric field in the following form:
+
+```math
+\begin{aligned}
+ & \nabla^2 \mathbf{E} - \frac{1}{c^2} \partial_t^2 \mathbf{E} = \mu_0  \partial_t\mathbf{J} \,, & 
+\qquad (1)
+\end{aligned}
+```
+
+and moreover we will use it as a scalar equation, for the field and current components along the laser field polarisation.
+
+#### Vacuum case
+
+In a case when no sources are present in the system, i.e. field propagates in vacuum, the equation simplifies down to the uniform wave equation
+
+```math
+\begin{aligned}
+& \left(\nabla^2 - \frac{1}{c^2} \; \partial_t^2\right) E = 0
+\end{aligned}
+```
+
+#### Plasma dispersion
+
+A simplest linear plasma response can be described with the help of non-relativistic motion equation:
+
+```math
+\begin{aligned}
+& \mu_0 \partial_t J = -e \; \mu_0 \;n_{pe} \; \partial_t v = \omega_{pe}^2 / c^2 E
+\end{aligned}
+```
+
+where $n_{pe}$ is a number density of plasma electrons, and $\omega_{pe} = \sqrt{\frac{e^2 n_{pe} }{\epsilon_0  m_e}}$ is a plasma frequency.
+
+Introducing this term in the optical equation leads to the following equation:
+
+```math
+\left(c^2 \nabla^2 - \partial_t^2 \; - \;  \omega_{pe}^2 \right) \; E =  0\,,
+```
+
+which if considered in the spatio-temporal Fourier domain leads to the well-known dispersion relation
+
+```math
+\omega^2 = k^2 c^2 + \omega_{pe}^2
+```
+
+# Field representations
+
+When approaching the computational problem of optical propagation we typically need to solve Eq. (1) with an initial condition that prepresents the initial beam state. There are two standard ways to define and consider the field.
+
+#### Spatial representation
+
+In so-called spatial form, the field is defined at a time $t_0$ in a full space as $E_0(x,y,z)=E(x,y,z, t_0)$, and we look for the state after some time pass, so we calculate, $E_1(x,y,z) = E(x,y,z, t_0+\Delta t)$. The spatial domain where $E_1$ is localized is usually adjusted to the propagated disctance, typically $z\to z+c\Delta t$ (we will always consider propagation along $z$-axis). This appropach is typical for the particle-in-cell (PIC) codes, where electromagnetic solver is couples with the charges particles motion, that can lead to the complex charge and current density structures in the spatial domain.
+
+#### Temporal representation
+
+Here the field is defined in the spatial plane located at $z_0$ as a function of time $E_0(x, y, t)=E(x, y, z_0, t)$, and we look at its state at a different location $E_1(x,y,t) = E(x,y,z_0+\Delta z, t)$ at a later time when the field arrives there, e.g. $t \to t+\Delta z/c$. In this approach, media is presented by an oscillating plane, and it is more convenient be used with the temporal spectral decompositions of the field, while resolving the longitudinal dynamics of the plasma on a scale of the pulse length $c\tau$, may be more challenging.
+
+While both approaches are mathematically equavalent, they would lead to different implementations, and usually serve different purposes. The spatial representation is used in PIC simulation as it is giving better description of the excited plasma structures and can be easilty extended to the more rigorous account of the source terms. **Temporal representation** is more common for the optical propagation problems, and in the following we will always be considering it as a **default** one. In vacuum without source terms, the conversion between representations can be easily done using propagation techniques.
+
+# Propagation equation with temporal representation
+
+Let us calculate the field $E_1(x,y,t) = E(x,y,z_0+\Delta z, t)$ for a known $E_0(x, y, t)=E(x, y, z_0, t)$, i.e. considering $z$ the propagation variable, and the field being defined in the temporal domain. For this we divide the Laplace operator in Eq(1) into longitudinal and transverse parts:
+
+```math
+\begin{aligned}
+& \partial_z^2 \; \mathbf{E} =  \frac{1}{c^2} \;\partial_t^2 \;\mathbf{E} \;-\; \nabla_\perp^2 \;\mathbf{E} + \mu_0 \; \partial_t \; \mathbf{J} \,, \qquad (2)
+\end{aligned}
+```
+
+where $\nabla_\perp^2 = \partial_x^2 + \partial_y^2$ for the problem defined in the cartesian geometry $(x,y,t)$, and $\nabla_{\perp}^2 = \cfrac{1}{r}\; \partial_r (r \; \partial_r) +  \cfrac{1}{r^2} \; \partial_\theta^2$ in the cylindrical coordiantes $(r, \theta, t)$.
+
+# Spectral methods
+
+All our computations will adopt the so-called *spectral* methods. In the optical problems this commonly involves presentation of the field and the source terms as the sums of the waves both temporally and in the transverse plane.
+
+For the **cartesian** geometry such decomposition is well-known as the spatio-temporal Fourier series:
+
+```math
+\begin{aligned}
+& A (x,y,t) = Re \left [ \sum_{ k_x, k_y, \omega} \; \hat{A}(k_x, k_y, \omega) \; \exp(ik_x x + ik_y y-i \omega t) \right ]\,,
+\end{aligned}
+```
+
+and for the **cylindrical** coordiantes we define is as the Fourier-Bessel (or Fourier-Hankel) series:
+
+```math
+\begin{aligned}
+& A (r,\theta,t) = Re \left [ \sum_{k_r, m, \omega} \; \hat{A}(k_r, m, \omega) \; J_m(k_r r) \; \exp(-i m \theta-i \omega t) \right ]
+\end{aligned}
+```
+
+The properties of such representations are defined by the sampling that we choose in the real and spectral spaces. For the case of spatio-temporal Fourier series, it is well known that uniform sampling in space
+
+Both these representations transform the first two terms in the right hand side of Eq(2) into:
+
+```math
+\frac{1}{c^2} \;\partial_t^2 \;\mathbf{E} \;-\; \nabla_\perp^2 \;\mathbf{E} \to -(\omega^2 - k_r^2) \hat{E}
+```
+
+## [Return to README](https://github.com/hightower8083/axiprop/blob/new-docs/README.md#documentation)

docs/propagators_vacuum.md

@@ -0,0 +1,19 @@
+## Propagation
+
+Propagation calculations in `axiprop` are realized with spectral transformations in time and transverse spatial domains.
+
+### Non-paraxial propagator
+
+#### 3D
+
+#### RZ
+
+### Fresnel propagator
+
+For a mode $k=2\pi / \lambda$
+
+$$ E(x, y, k, z) = \cfrac{ \exp\left(i k z \left(1 + \frac{x^2+y^2}{2 z^2}\right) \right) }{i \lambda z} \int  \int_{-\infty}^{\infty} 
+dx' dy' \left[E(x', y', k, z=0) \exp\left(\frac{i k}{2 z} \left(x'^2+y'^2\right) \right) \right] \exp\left(-\frac{i k}{z} 
+\left(xx'+yy'\right)\right) $$
+
+## [Return to README](https://github.com/hightower8083/axiprop/blob/new-docs/README.md#documentation)