Skip to content

GitLab

Commit d5d34d92 authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

hankel

parent 82a21c92
Pipeline #174470 failed with stages
in 29 seconds
Hide whitespace changes
Inline Side-by-side
examples/Vacuum.ipynb
View file @ d5d34d92
%% Cell type:code id: tags:
``` python
import fresnel.vacuum as vac
import fresnel.hankel
import numpy as np
import matplotlib.pyplot as plt
```
%% Cell type:code id: tags:
``` python
%matplotlib widget
```
%% Cell type:code id: tags:
``` python
n = 250
xx = np.linspace(-0.25, 0.25, n, endpoint=True).reshape(1,-1)
yy = np.linspace(-0.25, 0.25, n, endpoint=True).reshape(-1,1)
structuresize = 0.1
radius = structuresize/2
mask = (xx**2 < radius**2) * (yy**2 < radius**2)
u0 = mask * 1+0j
Fpix = 1/n
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots()
im = ax.imshow(np.abs(u0)**2)
fig.colorbar(im)
```
%%%% Output: display_data
%%%% Output: execute_result
<matplotlib.colorbar.Colorbar at 0x7fc8a219aa60>
<matplotlib.colorbar.Colorbar at 0x7f3746568670>
%% Cell type:code id: tags:
``` python
# Voelz sampling criterion: wl * z / dx / L = wl * z / dx**2 / N
# pixel Fresnel Number = dx**2 / z / wl
# Fpix is chosen such that the Voelz criterion is satisfied!
# For other fresnel numbers, the sampling must be adapted (padding)
```
%% Cell type:code id: tags:
``` python
# Transfer function propagator
propTF = vac.FresnelTFPropagator(u0.shape, Fpix)
upropTF = propTF(u0)
fig, ax = plt.subplots()
ax.plot(xx.flatten(), np.abs(u0[n//2,:])**2)
ax.plot(xx.flatten(), np.abs(upropTF[n//2,:])**2)
```
%%%% Output: display_data
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7fc8a2101340>]
[<matplotlib.lines.Line2D at 0x7f37464c0f10>]
%% Cell type:code id: tags:
``` python
propIR = vac.FresnelIRPropagator(u0.shape, Fpix)
upropIR = propIR(u0)
fig, ax = plt.subplots()
ax.plot(xx.flatten(), np.abs(u0[n//2,:])**2)
ax.plot(xx.flatten(), np.abs(upropIR[n//2,:])**2)
```
%%%% Output: display_data
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7fc8a1cf65e0>]
[<matplotlib.lines.Line2D at 0x7f37460c5310>]
%% Cell type:code id: tags:
``` python
mask2 = (xx**2 + yy**2 < radius**2)
u02 = mask2 * 1+0j
fig, ax = plt.subplots()
im = ax.imshow(np.abs(u02)**2)
fig.colorbar(im)
```
%%%% Output: display_data
%%%% Output: execute_result
<matplotlib.colorbar.Colorbar at 0x7fc8a1cda910>
<matplotlib.colorbar.Colorbar at 0x7f374602b520>
%% Cell type:code id: tags:
``` python
propTF2 = vac.FresnelTFPropagator(u0.shape, 2.5 * Fpix)
u2propTF = propTF2(u02)
fig, ax = plt.subplots()
ax.plot(xx.flatten(), np.abs(u02[n//2,:])**2)
ax.plot(xx.flatten(), np.abs(u2propTF[n//2,:])**2)
```
%%%% Output: stream
<ipython-input-49-f19e038b5538>:5: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
fig, ax = plt.subplots()
%%%% Output: display_data
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7fc8a1196a90>]
[<matplotlib.lines.Line2D at 0x7f3746062910>]
%% Cell type:markdown id: tags:
## Hankel Propagatoin
%% Cell type:code id: tags:
``` python
xx_hankel = fresnel.hankel.hankelSamples(n)
maskcs = (xx_hankel**2 + yy**2 < radius**2)
u02r = u02[n//2,:]
propCS = vac.FresnelPropagatorCS(u02r.shape[0], 2.5 * Fpix)
u2rprop = propCS(u02r)
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots()
ax.plot(xx.flatten(), np.abs(u02r)**2)
ax.plot(xx.flatten(), np.abs(u2rprop)**2)
ax.plot(xx_hankel)
```
%%%% Output: stream
%%%% Output: display_data
<ipython-input-51-29fa783d9e8a>:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
fig, ax = plt.subplots()
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7f3745da7dc0>]
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots()
ax.plot(xx.flatten(), np.abs(u02r)**2)
ax.plot(xx.flatten(), np.abs(u2rprop)**2)
```
%%%% Output: display_data
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7fc8a10f64c0>]
[<matplotlib.lines.Line2D at 0x7f3746020a00>]
......
fresnel/hankel.py
View file @ d5d34d92
......@@ -7,15 +7,15 @@ import numpy as np
import scipy.special
def hankelMatrix(N, n=0):
def hankelMatrix(N, order=0):
"""
Create a N x N matrix for discrete Hankel transfrom (DHT) [1]_ of n-th order.
Create a N x N matrix for discrete Hankel transfrom (DHT) [1]_ of order `order`.
Parameters
----------
N : int
number of pixels
n : int, optional
order : int, optional
order of DHT
Returns
......@@ -34,73 +34,105 @@ def hankelMatrix(N, n=0):
.. [1] N. Baddour and U. Chouinard, “Theory and operational rules for the discrete Hankel transform,” Journal of the Optical Society of America A, vol. 32, no. 4, p. 611, Mar. 2015. https://doi.org/10.1364/JOSAA.32.000611
"""
jn = np.array(scipy.special.jn_zeros(n, N + 1))
jn = np.array(scipy.special.jn_zeros(order, N + 1))
k = np.expand_dims(np.arange(N), axis=0)
m = np.expand_dims(np.arange(N), axis=1)
jN = jn[-1]
Y = scipy.special.jn(n, jn[m] * jn[k] / jN) # matrix
Y *= 2 / (jN * scipy.special.jn(n + 1, jn[k]) ** 2) # prefactor
Y = scipy.special.jn(order, jn[m] * jn[k] / jN) # matrix
Y *= 2 / (jN * scipy.special.jn(order + 1, jn[k]) ** 2) # prefactor
return Y
def hankelFreq(N, n=0, kmax=0.5):
def hankelSamples(N, kmax=0.5, order=0):
"""
Returns the Hankel space (frequency) sampling grid for the inverse DHT
Returns the real space sampling grid for the DHT
Parameters
----------
N : int
number of pixels
n : int, optional
order : int, optional
order of DHT
kmax : float, optional
maximum sampling frequency in dimensionless units
Returns
-------
array
frequency sampling grid
real space sampling grid
"""
jn = np.array(scipy.special.jn_zeros(n, N + 1))
jn = np.array(scipy.special.jn_zeros(order, N + 1))
return jn[:-1] * kmax / jn[N]
def hankelSamples(N, n=0, kmax=0.5):
def hankelFreq(N, kmax=0.5, order=0):
"""
Returns the real space sampling grid for the DHT
Returns the Hankel space (frequency) sampling grid for the inverse DHT
Parameters
----------
N : int
number of pixels
n : int, optional
order : int, optional
order of DHT
kmax : float, optional
maximum sampling frequency in dimensionless units
Returns
-------
array
real space sampling grid
frequency sampling grid
"""
jn = np.array(scipy.special.jn_zeros(n, N))
jn = np.array(scipy.special.jn_zeros(order, N))
return jn / (kmax * 2 * np.pi)
class DiscreteHankelTransform:
def __init__(self, N, n=0, kmax=0.5):
def __init__(self, N, kmax=0.5, order=0):
self._matrix = hankelMatrix(N, n)
self._matrix = hankelMatrix(N, order)
def __call__(self, x):
return self._matrix @ x
def hankel_resample_matrix(N1, new, kmax=0.5, n=0):
"""
As in: Reconstruction of optical fields with the Quasi-discrete Hankel transform
By: Andrew W. Norfolk and Edward J. Grace
https://www.osapublishing.org/DirectPDFAccess/31AA535C-9410-DDAC-67361BE3117DB160_199359/oe-18-10-10551.pdf?da=1&id=199359&seq=0&mobile=no
"""
jn = np.array(sci.special.jn_zeros(n, N1 + 1))
jN = jn[N1]
samples = jn[:-1] / kmax
K = kmax
k = np.expand_dims(np.arange(N), axis=0)
m = np.expand_dims(np.arange(N), axis=1)
samenm = samples[k] == abs(new[m])
samen = (samenm).sum(axis=1) > 0
same = samen[m]
S = (
2
* samples[k]
* sci.special.jn(n, K * new[m])
/ (K * (samples[k] ** 2 - new[m] ** 2 + samenm) * sci.special.jn(n + 1, K * samples[k]))
)
S = S * (~same) + same * samenm
return S
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment