Commit 80ba72fc authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

fix hankel propagation

parent d5d34d92
Pipeline #174480 passed with stages
in 1 minute and 44 seconds
%% 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 0x7f3746568670>
<matplotlib.colorbar.Colorbar at 0x7fefc1ca36d0>
%% 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 0x7f37464c0f10>]
[<matplotlib.lines.Line2D at 0x7fefc1bfbf70>]
%% 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 0x7f37460c5310>]
[<matplotlib.lines.Line2D at 0x7fefc17fd2e0>]
%% 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 0x7f374602b520>
<matplotlib.colorbar.Colorbar at 0x7fefc17634f0>
%% Cell type:code id: tags:
``` python
propTF2 = vac.FresnelTFPropagator(u0.shape, 2.5 * Fpix)
propTF2 = vac.FresnelTFPropagator(u0.shape, 2*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)
ax.set_xlim(0,0.25)
```
%%%% Output: display_data
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7f3746062910>]
(0.0, 0.25)
%% Cell type:markdown id: tags:
## Hankel Propagatoin
%% Cell type:code id: tags:
``` python
xx_hankel = fresnel.hankel.hankelSamples(n)
n_hankel = n//2
xx_hankel = fresnel.hankel.hankelSamples(n_hankel, xmax=0.25)
maskcs = (xx_hankel**2 + yy**2 < radius**2)
maskcs = (xx_hankel**2 < radius**2)
u02r = u02[n//2,:]
u02r = maskcs * 1+0j
propCS = vac.FresnelPropagatorCS(u02r.shape[0], 2.5 * Fpix)
propCS = vac.FresnelPropagatorCS(u02r.shape[0], 2*Fpix)
u2rprop = propCS(u02r)
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots()
ax.plot(xx_hankel)
ax.plot(xx_hankel, np.abs(u02r)**2)
ax.plot(xx_hankel, np.abs(u2rprop)**2)
ax.set_xlim(0,0.25)
```
%%%% Output: display_data
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7f3745da7dc0>]
(0.0, 0.25)
%% 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 0x7f3746020a00>]
......
......@@ -47,7 +47,7 @@ def hankelMatrix(N, order=0):
return Y
def hankelSamples(N, kmax=0.5, order=0):
def hankelSamples(N, xmax=0.5, order=0):
"""
Returns the real space sampling grid for the DHT
......@@ -57,19 +57,23 @@ def hankelSamples(N, kmax=0.5, order=0):
number of pixels
order : int, optional
order of DHT
kmax : float, optional
maximum sampling frequency in dimensionless units
xmax : float, optional
width of real space grid
Returns
-------
array
real space sampling grid
References
----------
.. [2] A. W. Norfolk and E. J. Grace, “Reconstruction of optical fields with the Quasi-discrete Hankel transform,” Optics Express, vol. 18, no. 10, p. 10551, May 2010. https://doi.org/10.1364%2Foe.18.010551
"""
jn = np.array(scipy.special.jn_zeros(order, N + 1))
return jn[:-1] * kmax / jn[N]
return jn[:-1] * xmax / jn[N]
def hankelFreq(N, kmax=0.5, order=0):
......@@ -91,9 +95,9 @@ def hankelFreq(N, kmax=0.5, order=0):
frequency sampling grid
"""
jn = np.array(scipy.special.jn_zeros(order, N))
jn = np.array(scipy.special.jn_zeros(order, N + 1))
return jn / (kmax * 2 * np.pi)
return jn[:-1] * kmax / jn[N]
class DiscreteHankelTransform:
......@@ -104,35 +108,3 @@ class DiscreteHankelTransform:
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
......@@ -9,7 +9,7 @@ from .misc import fftfreqn, gridn
def fresnelKernelCS(N, fresnelNumber):
hFreq = hankel.hankelFreq(N)
kern = np.exp(-1j * np.pi * hFreq ** 2 / fresnelNumber)
kern = np.exp(-1j * np.pi / fresnelNumber * hFreq ** 2)
return kern
......@@ -42,7 +42,7 @@ def fresnelTFKernel(shape, fresnelNumbers):
f = fftfreqn(shape)
kernel = np.ones(shape, dtype=np.complex128)
for dim in range(ndim):
kernel *= np.exp((-1j * np.pi / (fresnelNumbers[dim])) * f[dim] ** 2)
kernel *= np.exp(-1j * np.pi / fresnelNumbers[dim] * f[dim] ** 2)
return kernel
......
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