add C++ function for discrete hankel transform and some cleanup

.gitlab-ci.yml

@@ -20,7 +20,7 @@ cpplint:
   stage: lint
   script:
     - pip install cpplint
-    - cpplint --linelength=120 --filter=-legal/copyright,-build/include_order,-build/include_subdir source/*cpp source/*h
+    - cpplint --linelength=120 --filter=-legal/copyright source/*cpp source/*h
 
 pages:
   stage: deploy

CMakeLists.txt

 cmake_minimum_required (VERSION 3.15..3.17)
 project (fresnel VERSION 0.1)
 
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(default_build_type "Release")
 
@@ -15,21 +15,30 @@ add_subdirectory (extern/pybind11)
 find_package (Eigen3 3.3 REQUIRED NO_MODULE)
 find_package (OpenMP)
 
-set(python_module_name _solver)
+set(solver_module_name _solver)
 
-pybind11_add_module (${python_module_name} MODULE 
+pybind11_add_module (${solver_module_name} MODULE 
   source/solver_py.cpp
   )
-target_link_libraries (${python_module_name} PRIVATE Eigen3::Eigen)
-target_compile_options(${python_module_name} PRIVATE -Wpedantic -Wall -O3 -ffast-math)
+target_link_libraries (${solver_module_name} PRIVATE Eigen3::Eigen)
+target_compile_options(${solver_module_name} PRIVATE -Wpedantic -Wall -O3 -ffast-math)
 
 if(OpenMP_CXX_FOUND)
-    target_link_libraries(${python_module_name} PRIVATE OpenMP::OpenMP_CXX)
+    target_link_libraries(${solver_module_name} PRIVATE OpenMP::OpenMP_CXX)
 endif()
 
-install(TARGETS ${python_module_name} DESTINATION .)
-#install(DIRECTORY fresnel/ DESTINATION fresnel
-#        FILES_MATCHING PATTERN "*.py")
+set(hankel_module_name _hankel)
+pybind11_add_module (${hankel_module_name} MODULE 
+  source/hankel_py.cpp
+  )
+target_link_libraries (${hankel_module_name} PRIVATE Eigen3::Eigen)
+target_compile_options(${hankel_module_name} PRIVATE -Wpedantic -Wall -O3 -ffast-math)
+
+if(OpenMP_CXX_FOUND)
+    target_link_libraries(${hankel_module_name} PRIVATE OpenMP::OpenMP_CXX)
+endif()
+
+install(TARGETS ${solver_module_name} ${hankel_module_name} DESTINATION .)
 
 # Quiet a warning
 set(ignoreMe "${SKBUILD}")

fresnel/hankel.py

@@ -6,6 +6,8 @@ Discrete Hankel transform.
 import numpy as np
 import scipy.special
 
+from . import _hankel
+
 
 def hankelMatrix(N, order=0):
     """
@@ -97,13 +99,14 @@ def hankelFreq(N, kmax=0.5, order=0):
 
 
 class DiscreteHankelTransform:
-    def __init__(self, N, kmax=0.5, order=0):
+    def __init__(self, N, order=0):
 
-        self._matrix = hankelMatrix(N, order)
+        self._jn = np.asarray(scipy.special.jn_zeros(order, N + 1))
+        self._order = order
 
     def __call__(self, x):
 
-        return self._matrix @ x
+        return _hankel.dht(x, self._order, self._jn)
 
 
 class ResampleTransform:

fresnel/misc.py

 import numpy as np
-from pyfftw.interfaces.numpy_fft import fftfreq
+from scipy.fft import fftfreq
 
 
 def fftfreqn(N, d=1.0):

fresnel/vacuum.py

@@ -19,15 +19,14 @@ class FresnelPropagatorCS:
 
         self._N = N
 
-        Y = hankel.hankelMatrix(N, 0)
-        kern = fresnelKernelCS(N, fresnelNumber)
-
-        # TODO: express this nicer (kern is a diagonal matrix)
-        self._matrix = Y @ (kern[:, None] * Y)
+        self._kern = fresnelKernelCS(N, fresnelNumber)
+        self._dht = hankel.DiscreteHankelTransform(N, 0)
 
     def __call__(self, u):
 
-        uprop = self._matrix @ u
+        hf_u = self._dht(u)
+        hf_u *= self._kern
+        uprop = self._dht(hf_u)
 
         return uprop
 

source/finite_difference.h

 #pragma once
 
-#include <complex>
 #include <Eigen/Dense>
+#include <complex>
 #include <utility>  // std::move
 
 namespace algebra {

source/hankel.h

+#pragma once
+
+#include <Eigen/Dense>
+#include <cmath>
+#include <algorithm>
+
+
+namespace hankel {
+
+  template <typename T>
+  using array1d = Eigen::Array<T, Eigen::Dynamic, 1>;
+
+  template <typename T>
+    array1d<T> dht(const Eigen::Ref<const array1d<T>> in,
+                   const unsigned order,
+                   const Eigen::Ref<const array1d<double>> jn) {
+      const int n = std::min(in.size(), jn.size()-1);
+      const double jn_n = jn[n];
+      auto jn_v = jn.segment(0, n);
+
+      array1d<double> bessel_o1 = jn_v.unaryExpr(
+        [order](double x){
+          return std::cyl_bessel_j(order + 1, x);
+        });
+      array1d<T> in_mult = 2.0 / jn_n * bessel_o1.square().inverse() * in;
+
+      array1d<T> out(n, 1);
+      for (int i=0; i < n; ++i) {
+          array1d<double> row(n);
+          #pragma omp parallel for
+          for (int j=0; j < n; ++j) {
+              row[j] = std::cyl_bessel_j(order, jn_v[i] * jn_v[j] / jn_n);
+          }
+          out[i] = row.matrix().dot(in_mult.matrix());
+      }
+
+      return out;
+    }
+}  // namespace hankel

source/hankel_py.cpp

+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+#include <complex>
+
+#include "./hankel.h"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(_hankel, m) {
+  m.def("dht",
+        &hankel::dht<std::complex<double>>,
+        py::arg("in").noconvert(),
+        py::arg("order").noconvert(),
+        py::arg("jn").noconvert(),
+        py::return_value_policy::move);
+}

source/solver_py.cpp

 #include <pybind11/pybind11.h>
 #include <pybind11/eigen.h>
-#include <Eigen/Dense>
 
-#include "finite_difference.h"
+#include "./finite_difference.h"
 
 namespace py = pybind11;
-
 namespace fd = finite_differences;
 
 PYBIND11_MODULE(_solver, m) {