Coverage for yasfpy/functions/cpu_numba.py: 13%
167 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-02-15 20:36 +0100
1from numba import jit, prange, complex128, float64, int64
3import numpy as np
4from scipy.special import spherical_jn, hankel1, lpmv
7@jit(nopython=True, parallel=True, nogil=True, fastmath=True, cache=True)
8def particle_interaction(
9 lmax: int,
10 particle_number: int,
11 idx: np.ndarray,
12 x: np.ndarray,
13 translation_table: np.ndarray,
14 plm: np.ndarray,
15 sph_h: np.ndarray,
16 e_j_dm_phi,
17):
18 """Calculates the interaction between particles based on their properties and returns the result.
20 Args:
21 lmax (int): The maximum value of the angular momentum quantum number `l`. It determines the size of the arrays `plm` and `sph_h`.
22 particle_number (int): The number of particles in the system.
23 idx (np.ndarray): A numpy array of shape `(jmax, 5)`, where `jmax` is the total number of interactions between particles. Each row of `idx` represents an interaction and contains the following information:
24 - s1 (int): The index of the first particle.
25 - n1 (int): The index of the first particle's property.
26 - tau1 (int): The tau value of the first particle.
27 - l1 (int): The l value of the first particle.
28 - m1 (int): The m value of the first particle.
29 x (np.ndarray): A numpy array representing the positions of the particles. It has shape `(particle_number,)` and contains the x-coordinates of the particles.
30 translation_table (np.ndarray): A 3-dimensional numpy array that stores the translation coefficients used in the calculation. It has shape `(n2, n1, p)` where `n2` and `n1` are the indices of the translation coefficients, and `p` is the maximum.
31 plm (np.ndarray): A numpy array representing the associated Legendre polynomials. It has shape `(pmax * (pmax + 1) // 2, s1max, s2max)`, where `pmax` is the maximum degree of the Legendre polynomials.
32 sph_h (np.ndarray): A numpy array representing the spherical harmonics. It has shape `(jmax, jmax, channels)`, where `jmax` is the maximum number of particles, `channels` is the number of channels.
33 e_j_dm_phi (np.ndarray): The parameter `e_j_dm_phi` is not defined in the code snippet you provided. Could you please provide the definition or explanation of what `e_j_dm_phi` represents?
35 Returns:
36 wx (np.ndarray): The array `wx`, which represents the result of the particle interaction calculations.
37 """
38 jmax = particle_number * 2 * lmax * (lmax + 2)
39 channels = sph_h.shape[-1]
41 wx = np.zeros(x.size * channels, dtype=complex128).reshape(x.shape + (channels,))
43 for w_idx in prange(jmax * jmax * channels):
44 w = w_idx % channels
45 j_idx = w_idx // channels
46 j1 = j_idx // jmax
47 j2 = j_idx % jmax
48 s1, n1, tau1, l1, m1 = idx[j1, :]
49 s2, n2, tau2, l2, m2 = idx[j2, :]
51 if s1 == s2:
52 continue
54 delta_tau = np.absolute(tau1 - tau2)
55 delta_l = np.absolute(l1 - l2)
56 delta_m = np.absolute(m1 - m2)
58 val = 0j
59 for p in range(np.maximum(delta_m, delta_l + delta_tau), l1 + l2 + 1):
60 val += (
61 translation_table[n2, n1, p]
62 * plm[p * (p + 1) // 2 + delta_m, s1, s2]
63 * sph_h[p, s1, s2, w]
64 )
65 val *= e_j_dm_phi[m2 - m1 + 2 * lmax, s1, s2] * x[j2]
67 wx[j1, w] += val
69 return wx
72@jit(nopython=True, parallel=True, fastmath=True)
73def compute_idx_lookups(lmax: int, particle_number: int):
74 """
75 The function `compute_idx_lookups` generates an index lookup table for a given `lmax` and
76 `particle_number` using parallel processing.
78 Args:
79 lmax (int): The maximum value of the angular momentum quantum number `l`. It determines the range of values for `l` in the nested loop.
80 particle_number (int): The number of particles in the system.
82 Returns:
83 idx (np.ndarray): A NumPy array `idx` which contains the computed index lookups.
84 """
85 nmax = 2 * lmax * (lmax + 2)
86 idx = np.zeros(nmax * particle_number * 5, dtype=int64).reshape(
87 (nmax * particle_number, 5)
88 )
90 for s in prange(particle_number):
91 for tau in range(1, 3):
92 for l in range(1, lmax + 1):
93 for m in range(-l, l + 1):
94 n = (tau - 1) * lmax * (lmax + 2) + (l - 1) * (l + 1) + l + m
95 i = n + s * nmax
96 idx[i, 0] = s
97 idx[i, 1] = n
98 idx[i, 2] = tau
99 idx[i, 3] = l
100 idx[i, 4] = m
102 return idx
105@jit(nopython=True, parallel=True, nogil=True, fastmath=True)
106def compute_scattering_cross_section(
107 lmax: int,
108 particle_number: int,
109 idx: np.ndarray,
110 sfc: np.ndarray,
111 translation_table: np.ndarray,
112 plm: np.ndarray,
113 sph_h: np.ndarray,
114 e_j_dm_phi: np.ndarray,
115):
116 """Calculates the scattering cross section for a given set of input parameters.
118 Args:
119 lmax (int): The maximum angular momentum quantum number. It determines the maximum value of `l` in the calculations.
120 particle_number (int): The number of particles in the system.
121 idx (np.ndarray): A numpy array of shape `(jmax, 5)`, where `jmax` is the total number of particle pairs. Each row of `idx` represents a particle pair and contains the following information:
122 - s (int): The index of the first particle.
123 - n (int): The index of the second particle.
124 - tau (int): The tau value.
125 - l (int): The l value.
126 - m (int): The m value.
127 sfc (np.ndarray): A numpy array of shape `(s, n, channels)`, where:
128 translation_table (np.ndarray): A 3-dimensional numpy array that stores the translation coefficients used in the computation of the scattering cross section. It has shape `(n2, n1, p)` where `n2` and `n1` are the number of radial functions for the second and first particles, respectively, and `p` is the order of the Legendre polynomial.
129 plm (np.ndarray): A numpy array representing the associated Legendre polynomials. It has shape `(pmax * (pmax + 1) // 2, 2, 2)`, where `pmax` is the maximum value of `p` in the loop.
130 sph_h (np.ndarray): A numpy array of shape `(pmax, s1max, s2max, channels)`. It represents the scattering matrix elements for each combination of `s1`, `s2`, and `p`, where `p` is the order of the Legendre polynomial.
131 e_j_dm_phi (np.ndarray): A numpy array representing the scattering phase function. It has shape `(2*lmax+1, channels, channels)` and contains complex values. The indices `(j, s1, s2)` represent the angular momentum index `j`, and the spin indices `s1` and `s2`.
133 Returns:
134 c_sca_complex (np.ndarray): The complex scattering cross section `c_sca_complex`.
135 """
136 jmax = particle_number * 2 * lmax * (lmax + 2)
137 channels = sph_h.shape[-1]
139 c_sca_complex = np.zeros(channels, dtype=complex128)
141 for j_idx in prange(jmax * jmax):
142 j1 = j_idx // jmax
143 j2 = j_idx % jmax
144 s1, n1, _, _, m1 = idx[j1, :]
145 s2, n2, _, _, m2 = idx[j2, :]
147 delta_m = np.absolute(m1 - m2)
149 p_dependent = np.zeros(channels, dtype=complex128)
150 for p in range(delta_m, 2 * lmax + 1):
151 p_dependent += (
152 translation_table[n2, n1, p]
153 * plm[p * (p + 1) // 2 + delta_m, s1, s2]
154 * sph_h[p, s1, s2, :]
155 )
156 p_dependent *= (
157 np.conj(sfc[s1, n1, :])
158 * e_j_dm_phi[m2 - m1 + 2 * lmax, s1, s2]
159 * sfc[s2, n2, :]
160 )
162 c_sca_complex += p_dependent
164 return c_sca_complex
167@jit(nopython=True, parallel=True, nogil=True, fastmath=True)
168def compute_radial_independent_scattered_field_legacy(
169 lmax: int,
170 particles_position: np.ndarray,
171 idx: np.ndarray,
172 sfc: np.ndarray,
173 k_medium: np.ndarray,
174 azimuthal_angles: np.ndarray,
175 e_r: np.ndarray,
176 e_phi: np.ndarray,
177 e_theta: np.ndarray,
178 pilm: np.ndarray,
179 taulm: np.ndarray,
180):
181 """Calculates the scattered field for a given set of parameters and returns the result.
183 Args:
184 lmax (int): The maximum value of the angular momentum quantum number `l`. It determines the maximum order of the spherical harmonics used in the computation.
185 particles_position (np.ndarray): An array representing the positions of particles. It has shape `(num_particles, 3)`, where `num_particles` is the number of particles and each row represents the x, y, and z coordinates of a particle.
186 idx (np.ndarray): An array containing the indices of the particles. It has shape `(jmax, 5)` where `jmax` is the total number of particles. Each row of `idx` represents a particle and contains the following information:
187 sfc (np.ndarray): A 3-dimensional array representing the scattering form factors. It has dimensions `(s, n, w)`, where:
188 k_medium (np.ndarray): An array representing the wave number in the medium. It is used in the calculation of the scattered field.
189 azimuthal_angles (np.ndarray): An array of azimuthal angles, representing the angles at which the scattered field is computed.
190 e_r (np.ndarray): An array representing the radial component of the electric field. It has shape `(azimuthal_angles.size, 3)`, where `azimuthal_angles.size` is the number of azimuthal angles and 3 represents the three Cartesian components of the electric field.
191 e_phi (np.ndarray): An array representing the electric field component in the azimuthal direction. It has a shape of `(azimuthal_angles.size, 3)`, where `azimuthal_angles.size` is the number of azimuthal angles and `3` represents the three components of the electric field.
192 e_theta (np.ndarray): An array representing the electric field component in the theta direction. It has a shape of `(azimuthal_angles.size, 3)`, where `azimuthal_angles.size` is the number of azimuthal angles and `3` represents the three components of the electric field.
193 pilm (np.ndarray): An array representing the matrix of spherical harmonics coefficients. It has a shape of `(lmax+1, lmax+1, azimuthal_angles.size)`. Each element `pilm[l, m, a]` represents the coefficient of the spherical harmonics for a given `l`, `m`, and azimuthal angle `a`.
194 taulm (np.ndarray): An array representing the scattering coefficients for each combination of `l`, `m`, and azimuthal angle `a`. It has a shape of `(lmax+1, lmax+1, azimuthal_angles.size)`. The values in `taulm` represent the scattering coefficients.
196 Returns:
197 e_1_sca (np.ndarray): An array of complex numbers representing the scattered field.
199 """
200 e_1_sca = np.zeros(
201 azimuthal_angles.size * 3 * k_medium.size, dtype=complex128
202 ).reshape((azimuthal_angles.size, 3, k_medium.size))
203 jmax = particles_position.shape[0] * 2 * lmax * (lmax + 2)
205 for global_idx in prange(jmax * azimuthal_angles.size * k_medium.size):
206 w_idx = global_idx % (jmax * k_medium.size)
207 g_idx = global_idx // (jmax * k_medium.size)
209 a = g_idx
211 w = w_idx % k_medium.size
212 j_idx = w_idx // k_medium.size
213 s, n, tau, l, m = idx[j_idx, :]
215 t = (
216 np.power(1j, tau - l - 2)
217 * sfc[s, n, w]
218 / np.sqrt(2 * l * (l + 1))
219 * np.exp(
220 1j
221 * (
222 m * azimuthal_angles[a]
223 - k_medium[w] * np.sum(particles_position[s, :] * e_r[a, :])
224 )
225 )
226 )
228 if tau == 1:
229 e_1_sca[a, :, w] += t * (
230 e_theta[a, :] * pilm[l, np.abs(m), a] * 1j * m
231 - e_phi[a, :] * taulm[l, np.abs(m), a]
232 )
233 else:
234 e_1_sca[a, :, w] += t * (
235 e_phi[a, :] * pilm[l, np.abs(m), a] * 1j * m
236 + e_theta[a, :] * taulm[l, np.abs(m), a]
237 )
239 return e_1_sca
242@jit(nopython=True, parallel=True, nogil=True, fastmath=True)
243def compute_electric_field_angle_components(
244 lmax: int,
245 particles_position: np.ndarray,
246 idx: np.ndarray,
247 sfc: np.ndarray,
248 k_medium: np.ndarray,
249 azimuthal_angles: np.ndarray,
250 e_r: np.ndarray,
251 pilm: np.ndarray,
252 taulm: np.ndarray,
253):
254 """Calculates the electric field components in the theta and phi directions for given input parameters.
256 Args:
257 lmax (int): The maximum value of the angular momentum quantum number `l`. It determines the maximum value of `l` for which the calculations will be performed.
258 particles_position (np.ndarray): The positions of particles. It has shape `(num_particles, 3)`, where `num_particles` is the number of particles and each particle has 3 coordinates (x, y, z).
259 idx (np.ndarray): A numpy array of shape `(jmax, 5)`, where `jmax` is the total number of particles multiplied by `2 * lmax * (lmax + 2)`. Each row of `idx` represents the indices `(s, n, tau, l, m)`.
260 sfc (np.ndarray): A 3-dimensional numpy array representing the scattering form factors. It has dimensions `(s, n, w)`.
261 k_medium (np.ndarray): The wave vector in the medium. It is a numpy array that contains the wave vector values for different frequencies or wavelengths.
262 azimuthal_angles (np.ndarray): An array representing the azimuthal angles at which the electric field components are computed. It specifies the angles at which the electric field is measured in the azimuthal direction.
263 e_r (np.ndarray): The unit vector pointing in the direction of the electric field. It is a numpy array of shape `(azimuthal_angles.size, 3)`, where each row corresponds to a different azimuthal angle and the three columns represent the x, y, and z components.
264 pilm (np.ndarray): A 3-dimensional numpy array of shape `(lmax+1, lmax+1, azimuthal_angles.size)`. It represents the matrix elements of the electric field expansion coefficients for the theta component. The indices `(l, m, a)` correspond to the spherical harmonics.
265 taulm (np.ndarray): A numpy array that represents the angular momentum coupling coefficients. It has a shape of `(lmax+1, lmax+1, azimuthal_angles.size)`. The first dimension represents the value of `l`, the second dimension represents the value of `m`, and the third dimension represents the azimuthal angle.
267 Returns:
268 e_field_theta (np.ndarray): The electric field component in the theta direction.
269 e_field_phi (np.ndarray): The electric field component in the phi direction.
270 """
271 e_field_theta = np.zeros(
272 azimuthal_angles.size * k_medium.size, dtype=complex128
273 ).reshape((azimuthal_angles.size, k_medium.size))
274 e_field_phi = np.zeros_like(e_field_theta)
276 jmax = particles_position.shape[0] * 2 * lmax * (lmax + 2)
278 for global_idx in prange(jmax * azimuthal_angles.size * k_medium.size):
279 w_idx = global_idx % (jmax * k_medium.size)
280 g_idx = global_idx // (jmax * k_medium.size)
282 a = g_idx
284 w = w_idx % k_medium.size
285 j_idx = w_idx // k_medium.size
286 s, n, tau, l, m = idx[j_idx, :]
288 t = (
289 np.power(1j, tau - l - 2)
290 * sfc[s, n, w]
291 / np.sqrt(2 * l * (l + 1))
292 * np.exp(
293 1j
294 * (
295 m * azimuthal_angles[a]
296 - k_medium[w] * np.sum(particles_position[s, :] * e_r[a, :])
297 )
298 )
299 )
301 if tau == 1:
302 e_field_theta[a, w] += t * pilm[l, np.abs(m), a] * 1j * m
303 e_field_phi[a, w] -= t * taulm[l, np.abs(m), a]
304 else:
305 e_field_theta[a, w] += t * taulm[l, np.abs(m), a]
306 e_field_phi[a, w] += t * pilm[l, np.abs(m), a] * 1j * m
308 return e_field_theta, e_field_phi
311@jit(nopython=True, parallel=True, nogil=True, fastmath=True)
312def compute_polarization_components(
313 number_of_wavelengths: int,
314 number_of_angles: int,
315 e_field_theta: np.ndarray,
316 e_field_phi: np.ndarray,
317):
318 """
319 Compute the polarization components of electromagnetic fields.
321 Args:
322 number_of_wavelengths (int): The number of wavelengths.
323 number_of_angles (int): The number of angles.
324 e_field_theta (np.ndarray): The electric field component in the theta direction.
325 e_field_phi (np.ndarray): The electric field component in the phi direction.
327 Returns:
328 degree_of_polarization_tuple (tuple): A tuple containing the following polarization components:
329 - I (np.ndarray): The total intensity.
330 - degree_of_polarization (np.ndarray): The degree of polarization.
331 - degree_of_linear_polarization (np.ndarray): The degree of linear polarization.
332 - degree_of_linear_polarization_q (np.ndarray): The degree of linear polarization in the Q direction.
333 - degree_of_linear_polarization_u (np.ndarray): The degree of linear polarization in the U direction.
334 - degree_of_circular_polarization (np.ndarray): The degree of circular polarization.
335 """
336 # Stokes components
337 # S = np.zeros(4 * number_of_angles * number_of_wavelengths, dtype=complex128).reshape((4, number_of_angles, number_of_wavelengths))
338 I = np.zeros(number_of_angles * number_of_wavelengths, dtype=float64).reshape(
339 (number_of_angles, number_of_wavelengths)
340 )
341 Q = np.zeros_like(I)
342 U = np.zeros_like(I)
343 V = np.zeros_like(I)
345 for global_idx in prange(number_of_angles * number_of_wavelengths):
346 w_idx = global_idx % number_of_wavelengths
347 a_idx = global_idx // number_of_wavelengths
349 e_field_theta_abs = (
350 e_field_theta[a_idx, w_idx].real ** 2
351 + e_field_theta[a_idx, w_idx].imag ** 2
352 )
353 e_field_phi_abs = (
354 e_field_phi[a_idx, w_idx].real ** 2 + e_field_phi[a_idx, w_idx].imag ** 2
355 )
356 e_field_angle_interaction = (
357 e_field_theta[a_idx, w_idx] * e_field_phi[a_idx, w_idx].conjugate()
358 )
360 I[a_idx, w_idx] = e_field_theta_abs + e_field_phi_abs
361 Q[a_idx, w_idx] = e_field_theta_abs - e_field_phi_abs
362 U[a_idx, w_idx] = -2 * e_field_angle_interaction.real
363 V[a_idx, w_idx] = 2 * e_field_angle_interaction.imag
365 degree_of_polarization = np.sqrt(Q**2 + U**2 + V**2) / I
366 degree_of_linear_polarization = np.sqrt(Q**2 + U**2) / I
367 degree_of_linear_polarization_q = -Q / I
368 degree_of_linear_polarization_u = U / I
369 degree_of_circular_polarization = V / I
371 return (
372 I,
373 degree_of_polarization,
374 degree_of_linear_polarization,
375 degree_of_linear_polarization_q,
376 degree_of_linear_polarization_u,
377 degree_of_circular_polarization,
378 )
381@jit(nopython=True, parallel=True, nogil=True, fastmath=True)
382def compute_radial_independent_scattered_field(
383 number_of_wavelengths: int,
384 number_of_angles: int,
385 e_phi: np.ndarray,
386 e_theta: np.ndarray,
387 e_field_theta: np.ndarray,
388 e_field_phi: np.ndarray,
389):
390 """
391 Compute the radial independent scattered field.
393 Args:
394 number_of_wavelengths (int): The number of wavelengths.
395 number_of_angles (int): The number of angles.
396 e_phi (np.ndarray): The electric field in the phi direction.
397 e_theta (np.ndarray): The electric field in the theta direction.
398 e_field_theta (np.ndarray): The electric field theta component.
399 e_field_phi (np.ndarray): The electric field phi component.
401 Returns:
402 e_1_sca (np.ndarray): The computed radial independent scattered field.
403 """
404 e_1_sca = np.zeros(
405 number_of_angles * 3 * number_of_wavelengths, dtype=complex128
406 ).reshape((number_of_angles, 3, number_of_wavelengths))
408 for global_idx in prange(number_of_angles * number_of_wavelengths):
409 w = global_idx % number_of_wavelengths
410 a = global_idx // number_of_wavelengths
412 e_1_sca[a, :, w] = (
413 e_field_theta[a, w] * e_theta[a, :] + e_field_phi[a, w] * e_phi[a, :]
414 )
416 return e_1_sca
419@jit(parallel=True, forceobj=True)
420def compute_lookup_tables(
421 lmax: int, size_parameter: np.ndarray, phi: np.ndarray, cosine_theta: np.ndarray
422):
423 """
424 Compute lookup tables for spherical computations.
426 Args:
427 lmax (int): The maximum degree of the spherical harmonics.
428 size_parameter (np.ndarray): Array of size parameters.
429 phi (np.ndarray): Array of azimuthal angles.
430 cosine_theta (np.ndarray): Array of cosine of polar angles.
432 Returns:
433 spherical_bessel (np.ndarray): Array of spherical Bessel functions.
434 spherical_hankel (np.ndarray): Array of spherical Hankel functions.
435 e_j_dm_phi (np.ndarray): Array of exponential terms.
436 p_lm (np.ndarray): Array of associated Legendre polynomials.
437 """
438 spherical_hankel = np.zeros(
439 (2 * lmax + 1) * np.prod(size_parameter.shape), dtype=complex
440 ).reshape((2 * lmax + 1,) + size_parameter.shape)
441 spherical_bessel = np.zeros_like(spherical_hankel)
442 e_j_dm_phi = np.zeros(
443 (4 * lmax + 1) * np.prod(size_parameter.shape[:2]), dtype=complex
444 ).reshape((4 * lmax + 1,) + size_parameter.shape[:2])
445 p_lm = np.zeros(
446 (lmax + 1) * (2 * lmax + 1) * np.prod(size_parameter.shape[:2])
447 ).reshape(((lmax + 1) * (2 * lmax + 1),) + size_parameter.shape[:2])
449 for p in prange(2 * lmax + 1):
450 spherical_hankel[p, :, :, :] = np.sqrt(
451 np.divide(
452 np.pi / 2,
453 size_parameter,
454 out=np.zeros_like(size_parameter),
455 where=size_parameter != 0,
456 )
457 ) * hankel1(p + 1 / 2, size_parameter)
458 spherical_bessel[p, :, :, :] = spherical_jn(p, size_parameter)
459 e_j_dm_phi[p, :, :] = np.exp(1j * (p - 2 * lmax) * phi)
460 e_j_dm_phi[p + 2 * lmax, :, :] = np.exp(1j * p * phi)
461 for absdm in range(p + 1):
462 cml = np.sqrt(
463 (2 * p + 1) / 2 / np.prod(np.arange(p - absdm + 1, p + absdm + 1))
464 )
465 p_lm[p * (p + 1) // 2 + absdm, :, :] = (
466 cml * np.power(-1.0, absdm) * lpmv(absdm, p, cosine_theta)
467 )
469 return spherical_bessel, spherical_hankel, e_j_dm_phi, p_lm
472@jit(nopython=True, parallel=True, nogil=True, fastmath=True, cache=True)
473def compute_field(
474 lmax: int,
475 idx: np.ndarray,
476 size_parameter: np.ndarray,
477 sph_h: np.ndarray,
478 derivative: np.ndarray,
479 e_j_dm_phi: np.ndarray,
480 p_lm: np.ndarray,
481 pi_lm: np.ndarray,
482 tau_lm: np.ndarray,
483 e_r: np.ndarray,
484 e_theta: np.ndarray,
485 e_phi: np.ndarray,
486 scattered_field_coefficients: np.ndarray = None,
487 initial_field_coefficients: np.ndarray = None,
488 scatter_to_internal: np.ndarray = None,
489):
490 """
491 Compute the field using the given parameters and coefficients.
493 Parameters:
494 lmax (int): The maximum degree of the spherical harmonics.
495 idx (np.ndarray): The index array containing the values of s, n, tau, l, and m.
496 size_parameter (np.ndarray): The size parameter array.
497 sph_h (np.ndarray): The spherical harmonics array.
498 derivative (np.ndarray): The derivative array.
499 e_j_dm_phi (np.ndarray): The e_j_dm_phi array.
500 p_lm (np.ndarray): The p_lm array.
501 pi_lm (np.ndarray): The pi_lm array.
502 tau_lm (np.ndarray): The tau_lm array.
503 e_r (np.ndarray): The e_r array.
504 e_theta (np.ndarray): The e_theta array.
505 e_phi (np.ndarray): The e_phi array.
506 scattered_field_coefficients (np.ndarray, optional): The scattered field coefficients array. Defaults to None.
507 initial_field_coefficients (np.ndarray, optional): The initial field coefficients array. Defaults to None.
508 scatter_to_internal (np.ndarray, optional): The scatter to internal array. Defaults to None.
510 Returns:
511 field (np.ndarray): The computed field array.
512 """
513 jmax = sph_h.shape[1] * 2 * lmax * (lmax + 2)
514 channels = sph_h.shape[-1]
516 field = np.zeros(channels * sph_h.shape[2] * 3, dtype=complex128).reshape(
517 (channels, sph_h.shape[2], 3)
518 )
520 if (scattered_field_coefficients is None) and (initial_field_coefficients is None):
521 print(
522 "At least one, scattered field or initial field coefficients, need to be given."
523 )
524 print("Returning a zero array")
525 return field
527 for w_idx in prange(2 * lmax * (lmax + 2) * np.prod(np.array(sph_h.shape[1:]))):
528 w = w_idx % channels
529 j_idx = w_idx // channels
530 sampling_idx = j_idx // jmax
531 j_idx = j_idx % jmax
532 s, n, tau, l, m = idx[j_idx, :]
534 invariant = (
535 1 / np.sqrt(2 * (l + 1) * l) * e_j_dm_phi[m + 2 * lmax, s, sampling_idx]
536 )
537 # Calculate M
538 if tau == 1:
539 c_term_1 = (
540 1j
541 * m
542 * pi_lm[l, np.abs(m), s, sampling_idx]
543 * e_theta[s, sampling_idx, :]
544 )
545 c_term_2 = tau_lm[l, np.abs(m), s, sampling_idx] * e_phi[s, sampling_idx, :]
546 c_term = sph_h[l, s, sampling_idx, w] * (c_term_1 - c_term_2)
548 field[w, sampling_idx, :] += (
549 scattered_field_coefficients[s, n, w] * invariant * c_term
550 )
551 # Calculate N
552 else:
553 p_term = (
554 l
555 * (l + 1)
556 / size_parameter[s, sampling_idx, w]
557 * sph_h[l, s, sampling_idx, w]
558 * p_lm[l, np.abs(m), s, sampling_idx]
559 * e_r[s, sampling_idx, :]
560 )
561 # p_term = l * (l + 1) / size_parameter[s, sampling_idx, w]
562 # p_term *= sph_h[l, s, sampling_idx, w]
563 # p_term *= p_lm[l, np.abs(m), s, sampling_idx]
564 # p_term *= e_r[s, sampling_idx, :]
566 b_term_1 = (
567 derivative[l, s, sampling_idx, w] / size_parameter[s, sampling_idx, w]
568 )
569 b_term_2 = (
570 tau_lm[l, np.abs(m), s, sampling_idx] * e_theta[s, sampling_idx, :]
571 )
572 b_term_3 = (
573 1j
574 * m
575 * pi_lm[l, np.abs(m), s, sampling_idx]
576 * e_phi[s, sampling_idx, :]
577 )
578 b_term = b_term_1 * (b_term_2 + b_term_3)
580 field[w, sampling_idx, :] += (
581 scattered_field_coefficients[s, n, w] * invariant * (p_term + b_term)
582 )
584 return field