From 12e0188c6dfd15b05b1fb8d63ac5c04cce2d95e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 5 Mar 2019 16:11:15 +0000 Subject: [PATCH] C function to read scuff t-matrices. Former-commit-id: 26a01ac3f9aaf5a7dce2b0cb2a722a713af2f63c --- qpms/scatsystem.h | 52 ++++++++++++++++++- qpms/tmatrix_io.c | 126 ++++++++++++++++++++++++++++++++++++++++++++++ qpms/vswf.h | 12 +++++ 3 files changed, 188 insertions(+), 2 deletions(-) create mode 100644 qpms/tmatrix_io.c diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index e90507f..4d6e540 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -6,6 +6,7 @@ #include "qpms_types.h" #include #include +#include // 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 ); diff --git a/qpms/tmatrix_io.c b/qpms/tmatrix_io.c new file mode 100644 index 0000000..3bfb4d9 --- /dev/null +++ b/qpms/tmatrix_io.c @@ -0,0 +1,126 @@ +#define _POSIX_C_SOURCE 200809L // for getline() +#include "scatsystem.h" +#include +#include +#include +#include "indexing.h" +#include "vswf.h" // qpms_vswf_set_spec_find_uvswfi() +//#include + + +#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", + ¤tfreq_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; +} + + diff --git a/qpms/vswf.h b/qpms/vswf.h index ba22b3a..065a735 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -7,6 +7,7 @@ */ #ifndef QPMS_VSWF_H #define QPMS_VSWF_H +#include // ssize_t #include "qpms_types.h" #include @@ -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.