T-matrix reading function fix wrong nesting & other fixes.

Former-commit-id: e0a81510ec618fcc554f5ca9eee0b0a006970414

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, <cdouble[:len(self.spec),:len(self.spec)]>(t[0].m))

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
+    const gsl_interp_type *gsl_interp_linear
-    gsl_interp_type gsl_interp_cspline
+    const gsl_interp_type *gsl_interp_cspline
     # ^^^ These are probably the only relevant ones.
 
 cdef extern from "scatsystem.h":

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, 

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));
+        *freqs = realloc(*freqs, n_alloc * sizeof(double));
-          if (freqs_su) *freqs_su = realloc(*freqs_su, 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);
+        *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc);
-          if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
+        if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
-            qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
+          qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
-                  "realloc() failed.");
+              "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;
       }
+      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) {

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',