From 04923d0b452e2ab87f91c28059b13a8594625f8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 6 Mar 2019 21:05:52 +0000 Subject: [PATCH] T-matrix reading function fix wrong nesting & other fixes. Former-commit-id: e0a81510ec618fcc554f5ca9eee0b0a006970414 --- qpms/qpms_c.pyx | 7 +++++-- qpms/qpms_cdefs.pxd | 4 ++-- qpms/scatsystem.h | 4 ++-- qpms/tmatrix_io.c | 31 +++++++++++++++++-------------- setup.py | 4 +++- 5 files changed, 29 insertions(+), 21 deletions(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 1456af3..6dcb3e5 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1085,10 +1085,13 @@ cdef class TMatrixInterpolator: &(self.tmatrices_array), &(self.tmdata)): raise IOError("Could not read T-matrix from %s" % filename) self.interp = qpms_tmatrix_interpolator_create(self.nfreqs, - self.freqs, self.tmatrices_array, &gsl_interp_cspline) + self.freqs, self.tmatrices_array, gsl_interp_cspline) if not self.interp: raise Exception("Unexpected NULL at interpolator creation.") - def __call__(self, freq): + def __call__(self, double freq): '''Returns a TMatrix instance, corresponding to a given frequency.''' + if freq < self.freqs[0] or freq > self.freqs[self.nfreqs-1]:# FIXME here I assume that the input is already sorted + raise ValueError("input frequency %g is outside the interpolator domain (%g, %g)" + % (freq, self.freqs[0], self.freqs[self.nfreqs-1])) # This is a bit stupid, I should rethink the CTMatrix constuctors cdef qpms_tmatrix_t *t = qpms_tmatrix_interpolator_eval(self.interp, freq) cdef CTMatrix res = CTMatrix(self.spec, (t[0].m)) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 6dbf3e6..e7f34fe 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -217,8 +217,8 @@ cdef extern from "translations.h": cdef extern from "gsl/gsl_interp.h": struct gsl_interp_type: pass - gsl_interp_type gsl_interp_linear - gsl_interp_type gsl_interp_cspline + const gsl_interp_type *gsl_interp_linear + const gsl_interp_type *gsl_interp_cspline # ^^^ These are probably the only relevant ones. cdef extern from "scatsystem.h": diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 7ae946b..3ce6b37 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -166,7 +166,7 @@ qpms_errno_t qpms_load_scuff_tmatrix( const char *path, ///< Path to the TMatrix file const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec size_t *n, ///< Number of successfully loaded t-matrices - double **freqs, ///< Frequencies / (eV / ħ). + double **freqs, ///< Frequencies in SI units.. double **freqs_su, ///< Frequencies in SCUFF units (optional). /// The resulting T-matrices (optional). qpms_tmatrix_t **tmatrices_array, @@ -180,7 +180,7 @@ qpms_errno_t qpms_read_scuff_tmatrix( FILE *f, ///< An open stream with the T-matrix data. const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec size_t *n, ///< Number of successfully loaded t-matrices - double **freqs, ///< Frequencies / (eV / ħ). + double **freqs, ///< Frequencies in SI units. double **freqs_su, ///< Frequencies in SCUFF units (optional). /// The resulting T-matrices (optional). qpms_tmatrix_t **tmatrices_array, diff --git a/qpms/tmatrix_io.c b/qpms/tmatrix_io.c index e9972ae..2108de9 100644 --- a/qpms/tmatrix_io.c +++ b/qpms/tmatrix_io.c @@ -16,6 +16,10 @@ double qpms_SU2eV(double e_SU) { return e_SU * SCUFF_OMEGAUNIT / (ELECTRONVOLT / HBAR); } +double qpms_SU2SI(double e_SU) { + return e_SU * SCUFF_OMEGAUNIT; +} + /// TODO doc and more checks qpms_errno_t qpms_read_scuff_tmatrix( FILE *f, ///< file handle @@ -29,9 +33,9 @@ qpms_errno_t qpms_read_scuff_tmatrix( if (!(freqs && n && tmdata)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, "freqs, n, and tmdata are mandatory arguments and must not be NULL."); - if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT? + if (bs->norm != QPMS_NORMALISATION_POWER) // CHECKME CORRECT? qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)" + "Not implemented; only QPMS_NORMALISATION_POWER (CHECKME)" " norm supported right now."); int n_alloc = 128; // First chunk to allocate *n = 0; @@ -58,19 +62,18 @@ qpms_errno_t qpms_read_scuff_tmatrix( lastfreq_su = currentfreq_su; if(*n > n_alloc) { n_alloc *= 2; - *freqs = realloc(*freqs, n_alloc * sizeof(double)); - if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double)); - *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc); - if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) - qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "realloc() failed."); - - if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su; - (*freqs)[*n-1] = qpms_SU2eV(currentfreq_su); - - for(size_t i = 0; i < bs->n * bs->n; ++i) - (*tmdata)[(*n-1)*bs->n*bs->n + i] = NAN + I*NAN; + *freqs = realloc(*freqs, n_alloc * sizeof(double)); + if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double)); + *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc); + if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "realloc() failed."); } + if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su; + (*freqs)[*n-1] = qpms_SU2SI(currentfreq_su); + + for(size_t i = 0; i < bs->n * bs->n; ++i) + (*tmdata)[(*n-1)*bs->n*bs->n + i] = NAN + I*NAN; } qpms_vswf_type_t TAlpha, TBeta; switch(PAlpha) { diff --git a/setup.py b/setup.py index a8a26dc..068187d 100755 --- a/setup.py +++ b/setup.py @@ -26,7 +26,9 @@ qpms_c = Extension('qpms_c', 'qpms/wigner.c', 'qpms/scatsystem.c', 'qpms/vswf.c', # FIXME many things from vswf.c are not required by this module, but they have many dependencies (following in this list); maybe I want to move all the "basespec stuff" - 'qpms/legendre.c'], + 'qpms/legendre.c', + 'qpms/tmatrix_io.c', + 'qpms/error.c'], extra_compile_args=['-std=c99','-ggdb', '-O0', '-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required #'-DQPMS_USE_OMP',