T-matrix interpolator.

Former-commit-id: 3ec46a4afa7b92dd928f0b2156b5403780edd827

qpms/scatsystem.c

@@ -2,6 +2,7 @@
 #include <cblas.h>
 #include "scatsystem.h"
 #include "indexing.h"
+#include <gsl/gsl_spline.h>
 
 qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
   qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t));
@@ -128,10 +129,10 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) {
   return T;
 }
 
-bool qpms_tmatrix_isnear(const qpms_tmatrix *A, const qpms_tmatrix *B,
+bool qpms_tmatrix_isnear(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
     const double rtol, const double atol)
 {
-  if (!qpms_vswf_set_specisidentical(A->spec, B->spec))
+  if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec))
     return false;
   if (A->m == B->m)
     return true;
@@ -145,5 +146,75 @@ bool qpms_tmatrix_isnear(const qpms_tmatrix *A, const qpms_tmatrix *B,
 }
     
 
+qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incount,
+    const double *freqs, const qpms_tmatrix_t *ta, const gsl_interp_type *iptype//, const bool copy_bspec
+    ) {
+  if (incount <= 0) return NULL;
+  qpms_tmatrix_interpolator_t *ip = malloc(sizeof(qpms_tmatrix_interpolator_t));
+  /*
+  if (copy_bspec) {
+    ip->bspec = qpms_vswf_set_spec_copy(ta[0].spec);
+    ip->owns_bspec = true;
+  }
+  else {
+  */
+    ip->bspec = ta[0].spec;
+  //  ip->owns_bspec = false;
+  //}
+  const size_t n = ip->bspec->n;
 
+  // check if all matrices have the same bspec
+  for (size_t i = 0; i < incount; ++i)
+    if (!qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec))
+      abort();
+
+  if (!(ip->splines_real = calloc(n*n,sizeof(gsl_spline *)))) abort();
+  if (!(ip->splines_imag = calloc(n*n,sizeof(gsl_spline *)))) abort();
+  for (size_t row = 0; row < n; ++row)
+    for (size_t col = 0; col < n; ++col) {
+      double y_real[incount], y_imag[incount];
+      bool n0_real = false, n0_imag = false;
+      for (size_t i = 0; i < incount; ++i) {
+        complex double telem = ta[i].m[n * row + col];
+        if (y_real[i] = creal(telem)) n0_real = true;
+        if (y_imag[i] = cimag(telem)) n0_imag = true;
+      }
+      if (n0_real) {
+        gsl_spline *s =
+        ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount);
+        if (gsl_spline_init(s, freqs, y_real, incount) != 0 /*GSL_SUCCESS*/) abort();
+      }
+      else ip->splines_real[n * row + col] = NULL;
+     if (n0_imag) {
+        gsl_spline *s =
+        ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount);
+        if (gsl_spline_init(s, freqs, y_imag, incount) != 0 /*GSL_SUCCESS*/) abort();
+      }
+      else ip->splines_imag[n * row + col] = NULL;
+    }
+  return ip;
+}
+
+void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *ip) {
+  if (ip) {
+    const size_t n = ip->bspec->n;
+    for (size_t i = 0; i < n*n; ++i) {
+      if (ip->splines_real[i]) gsl_spline_free(ip->splines_real[i]);
+      if (ip->splines_imag[i]) gsl_spline_free(ip->splines_imag[i]);
+    }
+    //if (ip->owns_bspec)
+    //  qpms_vswf_set_spec_free(ip->bspec);
+    free(ip);
+  }
+}
+
+qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *ip, double freq) {
+  qpms_tmatrix_t *t = qpms_tmatrix_init(ip->bspec);
+  const size_t n = ip->bspec->n;
+  for (size_t i = 0; i < n*n; ++i){
+    if (ip->splines_real[i]) t->m[i] = gsl_spline_eval(ip->splines_real[i], freq, NULL /*does this work?*/);
+    if (ip->splines_imag[i]) t->m[i] += I* gsl_spline_eval(ip->splines_imag[i], freq, NULL /*does this work?*/);
+  }
+  return t;
+}
 

qpms/scatsystem.h

@@ -5,6 +5,7 @@
 #define QPMS_SCATSYSTEM_H
 #include "vswf.h"
 #include <stdbool.h>
+#include <gsl/gsl_spline.h>
 
 /// A T-matrix.
 /** In the future, I might rather use a more abstract approach in which T-matrix
@@ -141,18 +142,39 @@ qpms_errno_t qpms_load_scuff_tmatrix(
 		complex double **tmdata ///< The t-matrices raw contents
 		);
 
+/* Fuck this, include the whole <gsl/gsl_spline.h>
+typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct
+typedef struct gsl_interp_type gsl_interp_type;
+extern const gsl_interp_type * gsl_interp_linear;
+extern const gsl_interp_type * gsl_interp_polynomial;
+extern const gsl_interp_type * gsl_interp_cspline;
+extern const gsl_interp_type * gsl_interp_cspline_periodic;
+extern const gsl_interp_type * gsl_interp_akima;
+extern const gsl_interp_type * gsl_interp_akima_periodic;
+extern const gsl_interp_type * gsl_interp_steffen;
+*/
 
+// struct gsl_interp_accel; // use if lookup shows to be too slow
 typedef struct qpms_tmatrix_interpolator_t {
-	/*TODO; probably use gsl_spline from <gsl/gsl_spline.h> */;
+	const qpms_vswf_set_spec_t *bspec;
+	//bool owns_bspec;
+	gsl_spline **splines_real; ///< There will be a spline object for each nonzero element
+	gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element
+	// gsl_interp_accel **accel_real;
+	// gsl_interp_accel **accel_imag;
 } qpms_tmatrix_interpolator_t;
 
-/// NOT IMPLEMENTED
+/// Free a T-matrix interpolator.
 void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp);
 
-/// NOT IMPLEMENTED
+/// Evaluate a T-matrix interpolated value.
-qpms_tmatrix_t *qpms_tmatrix_interpolator_get(const qpms_tmatrix_interpolator_t *interp);
+qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq);
 
-/// NOT IMPLEMENTED 
+/// Create a T-matrix interpolator from frequency and T-matrix arrays.
-qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t *n, const double *freqs, const qpms_tmatrix_t *tmatrices_array /*, ...? */);
+qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided.
+		const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers!
+		const gsl_interp_type *iptype
+	       //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
+       	       /*, ...? */);
 
 #endif //QPMS_SCATSYSTEM_H

qpms/vswf.c

@@ -59,6 +59,16 @@ bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
   return true;
 }
 
+qpms_vswf_set_spec_t *qpms_vswf_set_spec_copy(const qpms_vswf_set_spec_t *or){
+  qpms_vswf_set_spec_t *c = malloc(sizeof(qpms_vswf_set_spec_t));
+  if (!c) abort; // return NULL
+  *c = *or;
+  c->ilist = malloc(sizeof(qpms_uvswfi_t) * c->n);
+  memcpy(c->ilist, or->ilist, sizeof(qpms_vswfi_t)*c->n);
+  c->capacity = c->n;
+  return c;
+}
+
 void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *s) {
   if(s) free(s->ilist);
   free(s);

qpms/vswf.h

@@ -41,6 +41,8 @@ void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *);
  */
 bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
 		const qpms_vswf_set_spec_t *b);
+/// Copies an instance of qpms_vswf_set_spec_t
+qpms_vswf_set_spec_t *qpms_vswf_set_spec_copy(const qpms_vswf_set_spec_t *orig);
 
 /// NOT IMPLEMENTED Evaluates a set of VSWF basis functions at a given point.
 /** The list of basis wave indices is specified in \a setspec;