add hankel resample matrix by Lars Melchior

fresnel/hankel.py

@@ -47,7 +47,7 @@ def hankelMatrix(N, order=0):
     return Y
 
 
-def hankelSamples(N, xmax=0.5, order=0):
+def hankelSamples(N, xmax=1, order=0):
     """
     Returns the real space sampling grid for the DHT
 
@@ -65,10 +65,6 @@ def hankelSamples(N, xmax=0.5, order=0):
     -------
     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))
@@ -108,3 +104,66 @@ class DiscreteHankelTransform:
     def __call__(self, x):
 
         return self._matrix @ x
+
+
+def hankel_resample_matrix(Nin, samplesout, xmax=1, order=0, eps=1e-10):
+    """
+    Create a matrix to resample a function sampled on the sampling grid of the DHT as in [2]_.
+
+    Parameters
+    ----------
+    Nin : int
+        number of samples
+    samplesout : array
+        new sampling grid
+    xmax : float
+        radius of the DHT sampling grid
+    order : int
+        order of the DHT
+    eps : float
+        regularization of the resampling
+
+    Returns
+    -------
+    ndarray
+        transformation matrix
+
+    See also
+    --------
+    fresnel.hankel.hankelSamples
+
+    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
+    """
+
+    Nout = len(samplesout)
+
+    jn = np.array(sci.special.jn_zeros(order, Nin + 1))
+    jN = jn[Nin]
+
+    K = jN / xmax
+
+    samplesin = jn[:-1] / K
+
+    samplesin = np.expand_dims(samplesin, axis=0)
+    samplesout = np.expand_dims(samplesout, axis=1)
+
+    k = np.arange(Nin)
+    m = np.arange(Nout)
+
+    diff = samplesin[k] ** 2 - samplesout[m] ** 2
+    singular_points = diff < eps
+
+    samen = (singular_points).sum(axis=1) > 0
+    same = samen[m]
+
+    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