C function to read scuff t-matrices.

Former-commit-id: 26a01ac3f9aaf5a7dce2b0cb2a722a713af2f63c

qpms/scatsystem.h

@@ -6,6 +6,7 @@
 #include "qpms_types.h"
 #include <stdbool.h>
 #include <gsl/gsl_spline.h>
+#include <stdio.h> // only because of qpms_read_scuff_tmatrix()
 
 /// A T-matrix.
 /** In the future, I might rather use a more abstract approach in which T-matrix
@@ -137,7 +138,54 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(
 		int N ///< number of z-axis rotations in the group
 		);
 
-/// NOT IMPLEMENTED Loads scuff-tmatrix generated files.
+/// Reads an open scuff-tmatrix generated file.
+/** 
+ * \a *freqs, \a *freqs_su, \a *tmatrices_array and \a *tmdata 
+ * arrays are allocated by this function
+ * and have to be freed by the caller after use.
+ * \a freqs_su and \a tmatrices_array can be NULL, in that case
+ * the respective arrays are not filled nor allocated.
+ *
+ * The contents of tmatrices_array is NOT 
+ * supposed to be freed element per element.
+ *
+ * TODO more checks and options regarding NANs etc.
+ *
+ */
+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_su, ///< Frequencies in SCUFF units (optional).
+		/// The resulting T-matrices (optional).
+		qpms_tmatrix_t **tmatrices_array, 
+		complex double **tmdata ///< The T-matrices raw contents
+		);
+
+/// Loads a scuff-tmatrix generated file.
+/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a 
+ * path instead of open FILE.
+ */
+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_su, ///< Frequencies in SCUFF units (optional).
+		/// The resulting T-matrices (optional).
+		qpms_tmatrix_t **tmatrices_array, 
+		/// The T-matrices raw contents. 
+		/** The coefficient of outgoing wave defined by 
+		 * \a bspec->ilist[desti] as a result of incoming wave
+		 * \a bspec->ilist[srci] at frequency \a (*freqs)[fi]
+		 * is accessed via
+		 * (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci].
+		 */
+		complex double ** tmdata 
+		);
+
+/// Loads scuff-tmatrix generated files.
 /** 
  * freqs, freqs_su, tmatrices_array and tmdata arrays are allocated by this function
  * and have to be freed by the caller after use.
@@ -149,7 +197,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
 		size_t *n, ///< Number of successfully loaded t-matrices
 		double **freqs, ///< Frequencies in SI units
 		double **freqs_su, ///< Frequencies in SCUFF units (optional)
-		qpms_tmatrix_t *tmatrices_array, ///< The resulting T-matrices.
+		qpms_tmatrix_t **tmatrices_array, ///< The resulting T-matrices.
 		complex double **tmdata ///< The t-matrices raw contents
 		);
 

qpms/tmatrix_io.c

@@ -0,0 +1,126 @@
+#define _POSIX_C_SOURCE 200809L // for getline()
+#include "scatsystem.h"
+#include <unistd.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "indexing.h" 
+#include "vswf.h" // qpms_vswf_set_spec_find_uvswfi()
+//#include <regex.h>
+
+
+#define HBAR (1.05457162825e-34)
+#define ELECTRONVOLT (1.602176487e-19)
+#define SCUFF_OMEGAUNIT (3e14)
+
+double qpms_SU2eV(double e_SU) {
+  return e_SU * SCUFF_OMEGAUNIT / (ELECTRONVOLT / HBAR);
+}
+
+/// TODO doc and more checks
+qpms_errno_t qpms_read_scuff_tmatrix(
+    FILE *f, ///< file handle
+    const qpms_vswf_set_spec_t * bs, ///< VSWF set spec
+    size_t *const n, ///< Number of successfully loaded t-matrices
+    double* *const freqs, ///< Frequencies in SI units
+    double* *const freqs_su, ///< Frequencies in SCUFF units (optional)
+    qpms_tmatrix_t* *const tmatrices_array, ///< The resulting T-matrices (optional).
+    complex double* *const tmdata
+    ) {
+  if (!(freqs && n && tmdata)) abort(); // These are mandatory
+  int n_alloc = 128; // First chunk to allocate
+  *n = 0;
+  *freqs = malloc(n_alloc * sizeof(double));
+  if (freqs_su) *freqs_su = malloc(n_alloc * sizeof(double));
+  *tmdata = malloc(sizeof(complex double) * bs->n * bs->n * n_alloc);
+  if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
+            abort(); // allocation failed
+  size_t linebufsz = 256;
+  char *linebuf = malloc(linebufsz);
+  ssize_t readchars;
+  double lastfreq_su = NAN;
+  while((readchars = getline(&linebuf, &linebufsz, f)) != -1) {
+    if (linebuf[0] == '#') continue;
+    int Alpha, LAlpha, MAlpha, PAlpha, Beta, LBeta, MBeta, PBeta;
+    double currentfreq_su, tr, ti;
+    if (11 != sscanf(linebuf, "%lf %d %d %d %d %d %d %d %d %lf %lf",
+          &currentfreq_su, &Alpha, &LAlpha, &MAlpha, &PAlpha,
+          &Beta, &LBeta, &MBeta, &PBeta, &tr, &ti))
+      abort(); // Malformed T-matrix file
+    if (currentfreq_su != lastfreq_su) { // New frequency -> new T-matrix
+      ++*n;
+      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)
+            abort(); // allocation 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;
+      }
+    }
+    qpms_vswf_type_t TAlpha, TBeta;
+    switch(PAlpha) {
+      case 0: TAlpha = QPMS_VSWF_MAGNETIC; break;
+      case 1: TAlpha = QPMS_VSWF_ELECTRIC; break;
+      default: assert(0);
+    }
+    switch(PBeta) {
+      case 0: TBeta = QPMS_VSWF_MAGNETIC; break;
+      case 1: TBeta = QPMS_VSWF_ELECTRIC; break;
+      default: assert(0);
+    }
+    qpms_uvswfi_t srcui = qpms_tmn2uvswfi(TAlpha, MAlpha, LAlpha),
+                  destui = qpms_tmn2uvswfi(TBeta, MBeta, LBeta);
+    ssize_t srci = qpms_vswf_set_spec_find_uvswfi(bs, srcui),
+            desti = qpms_vswf_set_spec_find_uvswfi(bs, destui);
+    if (srci == -1 || desti == -1)
+      /* This element has not been requested in bs->ilist. */
+      continue;
+    else
+      (*tmdata)[(*n-1)*bs->n*bs->n + desti*bs->n + srci] = tr + I*ti;
+  }
+  free(linebuf);
+  // free some more memory
+  n_alloc = *n;
+  *freqs = realloc(*freqs, n_alloc * sizeof(double));
+  if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double));
+  if (tmatrices_array) *tmatrices_array = realloc(*tmatrices_array, n_alloc * sizeof(qpms_tmatrix_t));
+  *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc);
+  if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
+    abort(); // allocation failed
+  if (tmatrices_array) {
+    *tmatrices_array = malloc(n_alloc * sizeof(qpms_tmatrix_t));
+    if (!*tmatrices_array) abort(); // allocation failed
+    for (size_t i = 0; i < *n; ++i) {
+      (*tmatrices_array)[i].spec = bs;
+      (*tmatrices_array)[i].m = (*tmdata) + i * bs->n * bs->n;
+      (*tmatrices_array)[i].owns_m = false;
+    }
+  }
+  return QPMS_SUCCESS;
+}
+
+qpms_errno_t qpms_load_scuff_tmatrix(
+    const char *path, ///< file path
+    const qpms_vswf_set_spec_t * bs, ///< VSWF set spec
+    size_t *const n, ///< Number of successfully loaded t-matrices
+    double **const freqs, ///< Frequencies in SI units
+    double ** const freqs_su, ///< Frequencies in SCUFF units (optional)
+    qpms_tmatrix_t ** const tmatrices_array, ///< The resulting T-matrices (optional).
+    complex double ** const tmdata
+    ) {
+  FILE *f = fopen(path, "r");
+  if (!f) abort(); // Could not open file.
+  qpms_errno_t retval = 
+    qpms_read_scuff_tmatrix(f, bs, n, freqs, freqs_su, tmatrices_array, tmdata);
+  if(fclose(f)) abort(); // Well, that would be weird. A read-only file failing to close.
+  return retval;
+}
+
+

qpms/vswf.h

@@ -7,6 +7,7 @@
  */
 #ifndef QPMS_VSWF_H
 #define QPMS_VSWF_H
+#include <unistd.h> // ssize_t 
 #include "qpms_types.h"
 #include <gsl/gsl_sf_legendre.h>
 
@@ -27,6 +28,17 @@ bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
 /// 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);
 
+/// Finds the position of a given index in the bspec's ilist.
+/** If not found, returns -1. */
+// TODO more consistency in types (here size_t vs. ptrdiff_t).
+static inline ssize_t qpms_vswf_set_spec_find_uvswfi(const qpms_vswf_set_spec_t *bspec,
+		const qpms_uvswfi_t index) {
+	for(size_t i = 0; i < bspec->n; ++i) 
+		if (bspec->ilist[i] == index)
+			return i;
+	return -1;
+}
+
 /// NOT IMPLEMENTED Evaluates a set of VSWF basis functions at a given point.
 /** The list of basis wave indices is specified in \a setspec; 
  *  \a setspec->norm must be set as well.