Coverage for frank/hankel.py: 78%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Frankenstein: 1D disc brightness profile reconstruction from Fourier data
2# using non-parametric Gaussian Processes
3#
4# Copyright (C) 2019-2020 R. Booth, J. Jennings, M. Tazzari
5#
6# This program is free software: you can redistribute it and/or modify
7# it under the terms of the GNU General Public License as published by
8# the Free Software Foundation, either version 3 of the License, or
9# (at your option) any later version.
10#
11# This program is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14# GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License
17# along with this program. If not, see <https://www.gnu.org/licenses/>
18#
19"""This module contains functions for computing the
20discrete Hankel transform (DHT).
21"""
22import numpy as np
23from scipy.special import j0, j1, jn_zeros, jv
25class DiscreteHankelTransform(object):
26 r"""
27 Utilities for computing the discrete Hankel transform.
29 This class provides the necessary interface to compute
30 a discrete version of the Hankel transform (DHT):
32 H[f](q) = \int_0^R_{max} f(r) J_nu(2*pi*q*r) * 2*pi*r dr.
34 The DHT is based on [1].
36 Additionally this class provides coefficients of the DHT [1] transform
37 matrix
39 Parameters
40 ----------
41 Rmax : float
42 Maximum radius beyond which f(r) is zero
43 N : integer
44 Number of terms to use in the series
45 nu : integer, default = 0
46 Order of the Bessel function, J_nu(r)
48 References
49 ----------
50 [1] Baddour & Chouinard (2015)
51 DOI: https://doi.org/10.1364/JOSAA.32.000611
52 Note: the definition of the DHT used here differs by factors
53 of 2*pi.
54 """
55 def __init__(self, Rmax, N, nu=0):
57 # Select the fast Bessel functions, if available
58 if nu == 0:
59 self._jnu0 = j0
60 self._jnup = j1
61 elif nu == 1:
62 self._jnu0 = j1
63 self._jnup = lambda x: jv(2, x)
64 else:
65 self._jnu0 = lambda x: jv(nu, x)
66 self._jnup = lambda x: jv(nu + 1, x)
68 self._N = N
69 self._nu = nu
71 # Compute the collocation points
72 j_nk = jn_zeros(nu, N + 1)
73 j_nk, j_nN = j_nk[:-1], j_nk[-1]
75 Qmax = j_nN / (2 * np.pi * Rmax)
77 self._Rnk = Rmax * (j_nk / j_nN)
78 self._Qnk = Qmax * (j_nk / j_nN)
80 self._Rmax = Rmax
81 self._Qmax = Qmax
83 # Compute the weights matrix
84 Jnk = np.outer(np.ones_like(j_nk), self._jnup(j_nk))
86 self._Ykm = (2 / (j_nN * Jnk * Jnk)) * \
87 self._jnu0(np.prod(np.meshgrid(j_nk, j_nk/j_nN), axis=0))
89 self._scale_factor = 1 / self._jnup(j_nk) ** 2
91 # Store the extra data needed
92 self._j_nk = j_nk
93 self._j_nN = j_nN
95 @classmethod
96 def get_collocation_points(cls, Rmax, N, nu=0):
97 """
98 Compute the collocation points for a Hankel Transform.
100 Parameters
101 ----------
102 Rmax : float
103 Maximum radius beyond which f(r) is zero
104 N : integer
105 Number of terms to use in the series
106 nu : integer, default = 0
107 Order of the Bessel function, J_nu(r)
109 Returns
110 -------
111 Rnk : array, shape=(N,) unit=radians
112 Radial collocation points in
113 Qnk : array, shape=(N,) unit=1/radians
114 Frequency collocation points
115 """
117 j_nk = jn_zeros(nu, N + 1)
118 j_nk, j_nN = j_nk[:-1], j_nk[-1]
120 Qmax = j_nN / (2 * np.pi * Rmax)
122 Rnk = Rmax * (j_nk / j_nN)
123 Qnk = Qmax * (j_nk / j_nN)
125 return Rnk, Qnk
127 def transform(self, f, q=None, direction='forward'):
128 """
129 Compute the Hankel transform of an array
131 Parameters
132 ----------
133 f : array, size = N
134 Function to Hankel transform, evaluated at the collocation points:
135 f[k] = f(r_k) or f[k] = f(q_k)
136 q : array or None
137 The frequency points at which to evaluate the Hankel
138 transform. If not specified, the conjugate points of the
139 DHT will be used. For the backwards transform, q should be
140 the radius points
141 direction : { 'forward', 'backward' }, optional
142 Direction of the transform. If not supplied, the forward
143 transform is used
145 Returns
146 -------
147 H[f] : array, size = N or len(q) if supplied
148 The Hankel transform of the array f
150 """
151 if q is None:
152 Y = self._Ykm
154 if direction == 'forward':
155 norm = (2 * np.pi * self._Rmax ** 2) / self._j_nN
156 elif direction == 'backward':
157 norm = (2 * np.pi * self._Qmax ** 2) / self._j_nN
158 else:
159 raise AttributeError("direction must be one of {}"
160 "".format(['forward', 'backward']))
161 else:
162 Y = self.coefficients(q, direction=direction)
163 norm = 1.0
165 return norm * np.dot(Y, f)
167 def coefficients(self, q=None, direction='forward'):
168 """
169 Coefficients of the transform matrix, defined by
170 H[f](q) = np.dot(Y, f)
172 Parameters
173 ----------
174 q : array or None
175 Frequency points at which to evaluate the transform. If q = None,
176 the points of the DHT are used. If direction='backward', these
177 points should instead be the radius points
178 direction : { 'forward', 'backward' }, optional
179 Direction of the transform. If not supplied, the forward transform
180 is used
182 Returns
183 -------
184 Y : array, size = (len(q), N)
185 The transformation matrix
186 """
187 if direction == 'forward':
188 norm = 1 / (np.pi * self._Qmax ** 2)
189 k = 1. / self._Qmax
190 elif direction == 'backward':
191 norm = 1 / (np.pi * self._Rmax ** 2)
192 k = 1. / self._Rmax
193 else:
194 raise AttributeError("direction must be one of {}"
195 "".format(['forward', 'backward']))
197 # For the DHT points we can use the cached Ykm points
198 if q is None:
199 return 0.5 * self._j_nN * norm * self._Ykm
201 H = (norm * self._scale_factor) * \
202 self._jnu0(np.outer(k * q, self._j_nk))
204 return H
206 def interpolation_coefficients(self, q, space='Real'):
207 """
208 Coefficients of the interpolation matrix, defined by
209 f(q) = np.dot(Y, f)
211 Parameters
212 ----------
213 q : array or None
214 The points at which to evaluate the interpolation.
215 space : { 'Real', 'Fourier' }, optional
216 Space in which the interpolation is done. If not supplied,
217 'Real' is assumed.
219 Returns
220 -------
221 Y : array, size = (len(q), N)
222 The interpolation matrix
223 """
224 if space == 'Real':
225 x = np.atleast_1d(2*np.pi * q * self._Qmax)
226 elif space == 'Fourier':
227 x = np.atleast_1d(2*np.pi * q * self._Rmax)
228 else:
229 raise ValueError("Space must be one of 'Real' or 'Fourier', not "
230 f"{space}.")
232 coeff = np.outer(np.where(x < self._j_nN, self._jnu0(x), 0),
233 2*self._j_nk / self._jnup(self._j_nk))
234 Y = coeff / (self._j_nk.reshape(1,-1)**2 - x.reshape(-1, 1)**2)
236 return Y
238 def interpolate(self, f, q, space='Real'):
239 """
240 Interpolate f (evaluated at the collocation points) to the new points,
241 pts, using interpolation that is consistent with the Fourier-Bessel
242 Series / Discrete Hankel Transform.
244 Parameters
245 ----------
246 f : array, size = N
247 Function to interpolate, evaluated at the collocation points:
248 f[k] = f(r_k) or f[k] = f(q_k)
249 q : array or None
250 The points at which to evaluate the interpolation.
251 space : { 'Real', 'Fourier' }, optional
252 Space in which the interpolation is done. If not supplied,
253 'Real' is assumed.
256 Returns
257 -------
258 f_interp : array, size = len(q)
259 The interpolated results
260 """
261 Y = self.interpolation_coefficients(q, space)
263 return np.dot(Y, f)
266 @property
267 def r(self):
268 """Radius points"""
269 return self._Rnk
271 @property
272 def Rmax(self):
273 """Maximum radius"""
274 return self._Rmax
276 @property
277 def q(self):
278 """Frequency points"""
279 return self._Qnk
281 @property
282 def Qmax(self):
283 """Maximum frequency"""
284 return self._Qmax
286 @property
287 def size(self):
288 """Number of points used in the DHT"""
289 return self._N
291 @property
292 def order(self):
293 """Order of the Bessel function"""
294 return self._nu