From f43db340752cbfbb203e3ade71767d3b0b76f430 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 21 Feb 2019 13:04:59 +0200 Subject: [PATCH] T-matrix interpolator. Former-commit-id: 3ec46a4afa7b92dd928f0b2156b5403780edd827 --- qpms/scatsystem.c | 75 +++++++++++++++++++++++++++++++++++++++++++++-- qpms/scatsystem.h | 34 +++++++++++++++++---- qpms/vswf.c | 10 +++++++ qpms/vswf.h | 2 ++ 4 files changed, 113 insertions(+), 8 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 5001b05..91873c9 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2,6 +2,7 @@ #include #include "scatsystem.h" #include "indexing.h" +#include 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; +} diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 4ff4588..1664d23 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -5,6 +5,7 @@ #define QPMS_SCATSYSTEM_H #include "vswf.h" #include +#include /// 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 +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 */; + 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 -qpms_tmatrix_t *qpms_tmatrix_interpolator_get(const qpms_tmatrix_interpolator_t *interp); +/// Evaluate a T-matrix interpolated value. +qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq); -/// NOT IMPLEMENTED -qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t *n, const double *freqs, const qpms_tmatrix_t *tmatrices_array /*, ...? */); +/// Create a T-matrix interpolator from frequency and T-matrix arrays. +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 diff --git a/qpms/vswf.c b/qpms/vswf.c index da237b0..e853d19 100644 --- a/qpms/vswf.c +++ b/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); diff --git a/qpms/vswf.h b/qpms/vswf.h index 3bf4d4c..79e8bdb 100644 --- a/qpms/vswf.h +++ b/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;