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

add hankel interpolation without storing matrix

parent 072dcf4e
Pipeline #175724 passed with stages
in 56 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
%matplotlib inline
```
%% 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 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 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 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 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)
```
%%%% Output: display_data
%%%% Output: execute_result
[<matplotlib.lines.Line2D at 0x7f6366be86d0>]
%% Cell type:markdown id: tags:
## Hankel Propagatoin
## Hankel Propagation
%% 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(n_hankel, 10*Fpix)
propCS = vac.FresnelPropagatorCS(n_hankel, 2*Fpix)
u2rprop = propCS(u02r)
interp = fresnel.hankel.interp1d(u2rprop, n_hankel, 0.25)
xx_new = np.linspace(-0.25, 0.25, n)
resample = fresnel.hankel.ResampleTransform(n_hankel, xx_new, 0.25)
u2rprop_new = resample(u2rprop)
u2rprop_new = interp(xx_new)
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots()
ax.plot(xx_new, np.abs(u2rprop_new)**2)
#ax.plot(np.abs(u2rprop)**2)
```
%%%% Output: display_data
%%%% Output: execute_result
%% Cell type:code id: tags:
[<matplotlib.lines.Line2D at 0x7f636679ea00>]
``` python
```
%% Cell type:code id: tags:
``` python
```
......
......@@ -36,7 +36,7 @@ def hankelMatrix(N, order=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(order, N + 1))
jn = np.asarray(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)
......@@ -69,7 +69,7 @@ def hankelSamples(N, xmax=1, order=0):
real space sampling grid
"""
jn = np.array(scipy.special.jn_zeros(order, N + 1))
jn = np.asarray(scipy.special.jn_zeros(order, N + 1))
return jn[:-1] * xmax / jn[N]
......@@ -93,7 +93,7 @@ def hankelFreq(N, kmax=0.5, order=0):
frequency sampling grid
"""
jn = np.array(scipy.special.jn_zeros(order, N + 1))
jn = np.asarray(scipy.special.jn_zeros(order, N + 1))
return jn[:-1] * kmax / jn[N]
......@@ -142,7 +142,7 @@ class ResampleTransform:
def __init__(self, N, samplesout, xmax=1, order=0, eps=1e-10):
jn = np.array(scipy.special.jn_zeros(order, N + 1))
jn = np.asarray(scipy.special.jn_zeros(order, N + 1))
jN = jn[N]
K = jN / xmax
......@@ -168,3 +168,35 @@ class ResampleTransform:
def __call__(self, f):
return self.matrix @ f
class interp1d:
def __init__(self, values, N, xmax=1, order=0, eps=1e-10):
self.eps = eps
jn = np.asarray(scipy.special.jn_zeros(order, N + 1))
jN = jn[N]
self.x = jn[:-1] * xmax / jN
self.y = values
self.K = jN / xmax
self.order = order
self.coeff = 2 * self.x / self.K / scipy.special.jn(order + 1, self.K * self.x)
def __call__(self, xnew):
factor = scipy.special.jn(self.order, self.K * xnew)
ynew = np.zeros(shape=xnew.shape, dtype=self.y.dtype)
for i in range(ynew.size):
diff = self.x ** 2 - xnew[i] ** 2
min_idx = np.argmin(np.abs(diff))
if np.abs(diff[min_idx]) <= self.eps:
ynew[i] = self.y[min_idx]
else:
ynew[i] = factor[i] * np.dot(self.coeff, self.y / diff)
return ynew
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