Add guile API for get-eigenmode-coefficients (#477)

* Add get-eigenmode-coefficients to scheme API

* Add support for kpoint functions

* Return kpoints from get-eigenmode-coefficients

* Remove extra define

scheme/examples/mode-coeffs.ctl

@@ -0,0 +1,99 @@
+
+(set-param! resolution 15)
+
+(define-param w 1)  ; width of waveguide
+(define-param L 10)  ; length of waveguide
+
+(define-param dair 3.0)
+(define-param dpml 3.0)
+
+;; mode frequency
+(define-param fcen 0.20)  ; > 0.5/sqrt(11) to have at least 2 modes
+
+(define (run-test mode-num kpoint-func)
+  (reset-meep)
+
+  (let* ((sx (+ dpml L dpml))
+        (sy (+ dpml dair w dair dpml))
+        (prism_x (+ sx 1))
+        (prism_y (/ w 2))
+        (verts
+         (list
+          (vector3 (* prism_x -1) prism_y)
+          (vector3 prism_x prism_y)
+          (vector3 prism_x (* prism_y -1))
+          (vector3 (* prism_x -1) (* prism_y -1))))
+        (modes-to-check '(1 2))  ; indices of modes for which to compute expansion coefficients
+        (xm (- (* 0.5 sx) dpml)))  ; x-coordinate of monitor
+
+      (set! geometry-lattice (make lattice (size sx sy no-size)))
+      (set! geometry
+            (list
+             (make prism
+               (center (vector3 0 0))
+               (vertices verts)
+               (height infinity)
+               (material (make medium (epsilon 12.0))))))
+
+      (set! pml-layers (list (make pml (thickness dpml))))
+
+
+      (set! sources (list
+                     (make eigenmode-source
+                       (src (make gaussian-src (frequency fcen) (fwidth (* 0.5 fcen))))
+                       (center (vector3 (+ (* -0.5 sx) dpml) 0))
+                       (component ALL-COMPONENTS)
+                       (size (vector3 0 (- sy (* 2 dpml))))
+                       (eig-match-freq? true)
+                       (eig-band mode-num)
+                       (eig-resolution 32))))
+
+      (set! symmetries
+            (list
+             (make mirror-sym (direction Y) (phase (if (odd? mode-num) 1 -1)))))
+
+      (let ((mflux (add-mode-monitor fcen 0 1
+                      (make mode-region (center (vector3 xm 0)) (size (vector3 0 (- sy (* 2 dpml)))))))
+            (mode-flux (add-flux fcen 0 1
+              (make flux-region (center (vector3 xm 0)) (size (vector3 0 (- sy (* 2 dpml))))))))
+
+        (run-sources+ 100)
+
+        (let* ((alpha-vgrp-kpoints (get-eigenmode-coefficients mflux modes-to-check #:kpoint-func kpoint-func))
+               (alpha (first alpha-vgrp-kpoints))
+               (vgrp (second alpha-vgrp-kpoints))
+               (kpoints (third alpha-vgrp-kpoints))
+               (mode-power (car (get-fluxes mode-flux)))
+               (test-passed #t)
+               (tolerance 5.0e-3)
+               (c0 (array-ref alpha (- mode-num 1) 0 0)))  ; coefficient of forward-traveling wave for mode # mode_num
+
+          (if (or (not (vector3-close? (first kpoints) (vector3 0.604301 0 0) 1e-7))
+                  (not (vector3-close? (second kpoints) (vector3 0.494353 0 0) 1e-2)))
+              (begin
+                (print "Test failed")
+                (exit 1)))
+
+          (map (lambda (nm)
+                 (if (not (= mode-num nm))
+                     (let ((cfrel (/ (magnitude (array-ref alpha (- nm 1) 0 0)) (magnitude c0)))
+                           (cbrel (/ (magnitude (array-ref alpha (- nm 1) 0 1)) (magnitude c0))))
+                       (if (or (> cfrel tolerance) (> cbrel tolerance))
+                           (set! test-passed #f)))))
+               modes-to-check)
+
+          ;; test 1: coefficient of excited mode >> coeffs of all other modes
+          (if (not test-passed)
+              (begin
+                (print "Test failed")
+                (exit 1)))
+
+          ;; test 2: |mode coeff|^2 = power
+          (if (> (abs (- 1 (/ mode-power (* (magnitude c0) (magnitude c0))))) 0.1)
+              (begin
+                (print "Test failed")
+                (exit 1)))))))
+
+(run-test 1 '())
+(run-test 2 '())
+(run-test 1 (lambda (freq mode) (vector3 0 0)))

scheme/meep-ctl-swig.hpp

@@ -5,6 +5,11 @@
 #ifndef MEEP_CTL_SWIG_HPP
 #define MEEP_CTL_SWIG_HPP 1
 
+struct kpoint_list {
+    meep::vec *kpoints;
+    size_t n;
+};
+
 vector3 vec_to_vector3(const meep::vec &);
 meep::vec vector3_to_vec(const vector3 v3);
 void set_dimensions(int dims);
@@ -23,12 +28,17 @@ meep::structure *make_structure(int dims, vector3 size, vector3 center,
 				double global_D_conductivity_diag_,
 				double global_B_conductivity_diag_);
 
-ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt, 
+ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt,
 				double fmin, double fmax, int maxbands,
                                 double spectral_density, double Q_thresh,
                                 double rel_err_thresh, double err_thresh,
                                 double rel_amp_thresh, double amp_thresh);
 
+kpoint_list do_get_eigenmode_coefficients(meep::fields *f, meep::dft_flux flux, const meep::volume &eig_vol,
+                                   int *bands, int num_bands, int parity, std::complex<double> *coeffs,
+                                   double *vgrp, double eig_resolution, double eigensolver_tol,
+                                   meep::kpoint_func user_kpoint_func, void *user_kpoint_data);
+
 ctlio::number_list dft_flux_flux(meep::dft_flux *f);
 ctlio::number_list dft_force_force(meep::dft_force *f);
 ctlio::number_list dft_ldos_ldos(meep::dft_ldos *f);

scheme/meep.cpp

@@ -33,7 +33,7 @@ void ctl_export_hook(void)
 
 /**************************************************************************/
 
-ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt, 
+ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt,
 				double fmin, double fmax, int maxbands,
                                 double spectral_density, double Q_thresh,
                                 double rel_err_thresh, double err_thresh,
@@ -66,6 +66,22 @@ ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt,
   return res;
 }
 
+kpoint_list do_get_eigenmode_coefficients(fields *f, dft_flux flux, const volume &eig_vol, int *bands, int num_bands,
+                                   int parity, std::complex<double> *coeffs, double *vgrp, double eig_resolution,
+                                   double eigensolver_tol, meep::kpoint_func user_kpoint_func,
+                                   void *user_kpoint_data)
+{
+  size_t num_kpoints = num_bands * flux.Nfreq;
+  meep::vec *kpoints = new meep::vec[num_kpoints];
+
+  f->get_eigenmode_coefficients(flux, eig_vol, bands, num_bands, parity, eig_resolution, eigensolver_tol,
+                                coeffs, vgrp, user_kpoint_func, user_kpoint_data, kpoints);
+
+  kpoint_list res = {kpoints, num_kpoints};
+
+  return res;
+}
+
 /**************************************************************************/
 
 /* This is a wrapper function to fool SWIG...since our list constructor

scheme/meep.i

@@ -60,52 +60,60 @@ static inline std::complex<double> my_complex_func3(std::complex<double> x) {
   return std::complex<double>(cret.re, cret.im);
 }
 
+static meep::vec my_kpoint_func(double freq, int mode, void *user_data) {
+  SCM scm_res = gh_call2((SCM)user_data, ctl_convert_number_to_scm(freq),
+                         ctl_convert_number_to_scm(mode));
+  vector3 v3 = ctl_convert_vector3_to_c(scm_res);
+  meep::vec result(v3.x, v3.y, v3.z);
+  return result;
+}
+
 %}
 
 %typecheck(SWIG_TYPECHECK_COMPLEX) std::complex<double> {
   $1 = SwigComplex_Check($input);
 }
 
-%typemap(guile,out) complex, std::complex<double>, std::complex<double> {
+%typemap(out) complex, std::complex<double>, std::complex<double> {
   $result = scm_make_rectangular(ctl_convert_number_to_scm($1.real()),
 				 ctl_convert_number_to_scm($1.imag()));
 }
-%typemap(guile,in) complex, std::complex<double>, std::complex<double> {
+%typemap(in) complex, std::complex<double>, std::complex<double> {
   cnumber cnum = ctl_convert_cnumber_to_c($input);
   $1 = std::complex<double>(cnum.re, cnum.im);
 }
 
-%typemap(guile,in) std::complex<double>(*)(meep::vec const &) {
+%typemap(in) std::complex<double>(*)(meep::vec const &) {
   my_complex_func_scm = $input;
   $1 = my_complex_func;
 }
 %typecheck(SWIG_TYPECHECK_POINTER) std::complex<double>(*)(meep::vec const &) {
   $1 = SCM_NFALSEP(scm_procedure_p($input));
 }
 
-%typemap(guile,in) std::complex<double>(*)(std::complex<double>) {
+%typemap(in) std::complex<double>(*)(std::complex<double>) {
   my_complex_func3_scm = $input;
   $1 = my_complex_func3;
 }
 %typecheck(SWIG_TYPECHECK_POINTER) std::complex<double>(*)(std::complex<double>) {
   $1 = SCM_NFALSEP(scm_procedure_p($input));
 }
 
-%typemap(guile,in) (std::complex<double> (*func)(double t, void *), void *data) {
+%typemap(in) (std::complex<double> (*func)(double t, void *), void *data) {
   $1 = my_complex_func2;
   $2 = (void *) $input; // input is SCM pointer to Scheme function
 }
 %typecheck(SWIG_TYPECHECK_POINTER) (std::complex<double> (*func)(double t, void *), void *data) {
   $1 = SCM_NFALSEP(scm_procedure_p($input));
 }
 
-%typemap(guile,in) meep::vec {
+%typemap(in) meep::vec {
   $1 = vector3_to_vec(ctl_convert_vector3_to_c($input));
 }
-%typemap(guile,out) meep::vec {
+%typemap(out) meep::vec {
   $result = ctl_convert_vector3_to_scm(vec_to_vector3($1));
 }
-%typemap(guile,in) meep::vec const & %{
+%typemap(in) meep::vec const & %{
   meep::vec vec__$1 = vector3_to_vec(ctl_convert_vector3_to_c($input));
   $1 = &vec__$1;
 %}
@@ -115,7 +123,7 @@ static inline std::complex<double> my_complex_func3(std::complex<double> x) {
 
 /* field_function arguments are passed as a cons pair of (components . func)
    in order to set all four arguments at once. */
-%typemap(guile,in) (int num_fields, const meep::component *components, meep::field_function fun, void *fun_data_) (my_field_func_data data) {
+%typemap(in) (int num_fields, const meep::component *components, meep::field_function fun, void *fun_data_) (my_field_func_data data) {
   $1 = list_length(gh_car($input));
   $2 = new meep::component[$1];
   for (int i = 0; i < $1; ++i)
@@ -137,7 +145,7 @@ static inline std::complex<double> my_complex_func3(std::complex<double> x) {
 /* integrate2 arguments are passed as a cons pair of
    ((components1 . components2) . func)
    in order to set all six arguments at once. */
-%typemap(guile,in) (int num_fields1, const meep::component *components1, int num_fields2, const meep::component *components2, meep::field_function integrand, void *integrand_data_) (my_field_func_data data) {
+%typemap(in) (int num_fields1, const meep::component *components1, int num_fields2, const meep::component *components2, meep::field_function integrand, void *integrand_data_) (my_field_func_data data) {
   $1 = list_length(gh_car(gh_car($input)));
   $2 = new meep::component[$1];
   for (int i = 0; i < $1; ++i)
@@ -182,6 +190,64 @@ static inline std::complex<double> my_complex_func3(std::complex<double> x) {
   }
 }
 
+// do_get_eigenmode_coefficients
+
+%typemap(in) (int *bands, int num_bands) {
+  $2 = list_length($input);
+  $1 = new int[$2];
+
+  for (int i = 0; i < $2; ++i) {
+    $1[i] = integer_list_ref($input, i);
+  }
+}
+
+%typemap(freearg) (int *bands, int num_bands) {
+  if ($1) {
+    delete[] $1;
+  }
+}
+
+%typemap(in) (meep::kpoint_func user_kpoint_func, void *user_kpoint_data) {
+  if (SCM_NFALSEP(scm_procedure_p($input))) {
+    $1 = my_kpoint_func;
+    $2 = (void*)$input;
+  }
+  else {
+    $1 = NULL;
+    $2 = NULL;
+  }
+}
+
+%typemap(in, noblock=1) std::complex<double> *coeffs {
+  scm_t_array_handle coeffs_handle;
+  scm_array_get_handle($input, &coeffs_handle);
+  $1 = (std::complex<double>*)scm_array_handle_uniform_writable_elements(&coeffs_handle);
+}
+
+%typemap(in, noblock=1) double *vgrp {
+  scm_t_array_handle vgrp_handle;
+  scm_array_get_handle($input, &vgrp_handle);
+  $1 = (double*)scm_array_handle_uniform_writable_elements(&vgrp_handle);
+}
+
+%typemap(freearg) std::complex<double> *coeffs {
+  scm_array_handle_release(&coeffs_handle);
+}
+
+%typemap(freearg) double *vgrp {
+  scm_array_handle_release(&vgrp_handle);
+}
+
+%typemap(out) kpoint_list {
+  int i;
+  list cur_list = SCM_EOL;
+  for (i = $1.n - 1; i >= 0; --i) {
+    cur_list = gh_cons(vector32scm(vec_to_vector3($1.kpoints[i])), cur_list);
+  }
+  $result = cur_list;
+  delete[] $1.kpoints;
+}
+
 // Need to tell SWIG about any method that returns a new object
 // which needs to be garbage-collected.
 %newobject meep::fields::open_h5file;

scheme/meep.scm.in

@@ -579,6 +579,26 @@
   (load-flux fname flux)
   (meep-dft-flux-scale-dfts flux -1.0))
 
+; ****************************************************************
+; Mode monitor
+
+(define mode-region flux-region)
+
+(define (fields-add-mode-monitor fields fcen df nfreq . fluxes)
+  (if (not (= (length fluxes) 1)) (error "add-mode-monitor expected just one mode-region."))
+  (let* ((region (car fluxes))
+         (v (volume (center (object-property-value region 'center))
+                    (size (object-property-value region 'size))))
+         (d0 (object-property-value region 'direction))
+         (d (if (< d0 0)
+                (meep-fields-normal-direction fields v)
+                d0)))
+    (meep-fields-add-mode-monitor fields d v (- fcen (/ df 2)) (+ fcen (/ df 2)) nfreq)))
+
+(define (add-mode-monitor fcen df nfreq . fluxes)
+  (if (null? fields) (init-fields))
+  (apply fields-add-mode-monitor (append (list fields fcen df nfreq) fluxes)))
+
 ; ****************************************************************
 ; Force spectra (from stress tensor) - very similar interface to flux spectra
 
@@ -1084,6 +1104,32 @@
 (define (harminv c pt fcen df . mxbands)
   (harminv! harminv-data harminv-data-dt harminv-results c pt fcen df mxbands))
 
+
+; ****************************************************************
+; get-eigenmode-coefficients
+
+(define (get-keyword-value args keyword default)
+  (let ((kv (memq keyword args)))
+    (if (and kv (>= (length kv) 2))
+        (cadr kv)
+        default)))
+
+(define (get-eigenmode-coefficients flux bands . args)
+  (if (null? fields)
+      (error "init-fields is required before using get-eigenmode-coefficients"))
+  (let ((eig-parity (get-keyword-value args #:eig-parity NO-PARITY))
+        (eig-vol (get-keyword-value args #:eig-vol (meep-dft-flux-where-get flux)))
+        (eig-resolution (get-keyword-value args #:eig-resolution 0))
+        (eig-tolerance (get-keyword-value args #:eig-tolerance 1e-7))
+        (num-bands (length bands))
+        (kpoint-func (get-keyword-value args #:kpoint-func '())))
+    (begin
+      (define coeffs (make-typed-array 'c64 '0 num-bands (meep-dft-flux-Nfreq-get flux) 2))
+      (define vgrp (make-typed-array 'f64 '0 num-bands (meep-dft-flux-Nfreq-get flux)))
+      (define kpoints (do-get-eigenmode-coefficients fields flux eig-vol bands eig-parity coeffs vgrp
+                                                     eig-resolution eig-tolerance kpoint-func))
+      (list coeffs vgrp kpoints))))
+
 ; ****************************************************************
 ; dft-ldos step function