Skip to content

GitLab

Commit 86f7c517 authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

fix resampling

parent 9410e501
Pipeline #174548 failed with stages
in 56 seconds
Hide whitespace changes
Inline Side-by-side
examples/Vacuum.ipynb
View file @ 86f7c517
%% 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 0x7fefc1ca36d0>
<matplotlib.colorbar.Colorbar at 0x7f63670f0760>
%% 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 0x7fefc1bfbf70>]
[<matplotlib.lines.Line2D at 0x7f63670b6040>]
%% 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 0x7fefc17fd2e0>]
[<matplotlib.lines.Line2D at 0x7f6366c49190>]
%% 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 0x7fefc17634f0>
<matplotlib.colorbar.Colorbar at 0x7f6366baf5b0>
%% Cell type:code id: tags:
``` python
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)
#ax.set_xlim(0,0.25)
```
%%%% Output: display_data
%%%% Output: execute_result
(0.0, 0.25)
[<matplotlib.lines.Line2D at 0x7f6366be86d0>]
%% Cell type:markdown id: tags:
## Hankel Propagatoin
%% Cell type:code id: tags:
``` python
n_hankel = n//2
xx_hankel = fresnel.hankel.hankelSamples(n_hankel, xmax=0.25)
maskcs = (xx_hankel**2 < radius**2)
u02r = maskcs * 1+0j
propCS = vac.FresnelPropagatorCS(u02r.shape[0], 2*Fpix)
propCS = vac.FresnelPropagatorCS(n_hankel, 10*Fpix)
u2rprop = propCS(u02r)
xx_new = np.linspace(-0.25, 0.25, n)
resample = fresnel.hankel.ResampleTransform(n_hankel, xx_new, 0.25)
u2rprop_new = resample(u2rprop)
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots()
ax.plot(xx_hankel, np.abs(u02r)**2)
ax.plot(xx_hankel, np.abs(u2rprop)**2)
ax.set_xlim(0,0.25)
ax.plot(xx_new, np.abs(u2rprop_new)**2)
```
%%%% Output: display_data
%%%% Output: execute_result
(0.0, 0.25)
[<matplotlib.lines.Line2D at 0x7f636679ea00>]
%% Cell type:code id: tags:
``` python
```
......
fresnel/hankel.py
View file @ 86f7c517
......@@ -106,9 +106,9 @@ class DiscreteHankelTransform:
return self._matrix @ x
def hankel_resample_matrix(Nin, samplesout, xmax=1, order=0, eps=1e-10):
class ResampleTransform:
"""
Create a matrix to resample a function sampled on the sampling grid of the DHT as in [2]_.
Create a callable to resample a function sampled on the sampling grid of the DHT as in [2]_.
Parameters
----------
......@@ -137,33 +137,32 @@ def hankel_resample_matrix(Nin, samplesout, xmax=1, order=0, eps=1e-10):
.. [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
"""
Nout = len(samplesout)
def __init__(self, Nin, samplesout, xmax=1, order=0, eps=1e-10):
Nout = len(samplesout)
jn = np.array(sci.special.jn_zeros(order, Nin + 1))
jN = jn[Nin]
jn = np.array(scipy.special.jn_zeros(order, Nin + 1))
jN = jn[Nin]
K = jN / xmax
K = jN / xmax
samplesin = jn[:-1] / K
samplesin = jn[:-1] * xmax / jN
samplesin = np.expand_dims(samplesin, axis=0)
samplesout = np.expand_dims(samplesout, axis=1)
samplesin = np.expand_dims(samplesin, axis=0)
samplesout = np.expand_dims(samplesout, axis=1)
k = np.arange(Nin)
m = np.arange(Nout)
diff = samplesin ** 2 - samplesout ** 2
singular_points = np.abs(diff) < eps
diff = samplesin[k] ** 2 - samplesout[m] ** 2
singular_points = diff < eps
samen = singular_points.sum(axis=1) > 0
same = np.expand_dims(samen, axis=1)
samen = (singular_points).sum(axis=1) > 0
same = samen[m]
tmp = (
2
* samplesin
* scipy.special.jn(order, K * samplesout)
/ (K * (diff + singular_points) * scipy.special.jn(order + 1, K * samplesin))
)
self.matrix = (~same) * tmp + same * singular_points
tmp = (
2
* samplesin[k]
* sci.special.jn(order, K * samplesout[m])
/ (K * (diff + samenm) * sci.special.jn(order + 1, K * samplesin[k]))
)
matrix = (~same) * tmp + same * singular_points
return matrix
def __call__(self, f):
return self.matrix @ f
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