Move material properties-related stuff to materials.[ch]

Former-commit-id: 0039b167488c5aca55ea1912e43ec3ab8b289c1e

qpms/CMakeLists.txt

@@ -12,7 +12,7 @@ include_directories(${DIRS})
 
 add_library (qpms translations.c tmatrices.c vecprint.c vswf.c wigner.c
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
-	bessel.c own_zgemm.c parsing.c scatsystem.c)
+	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c)
 use_c99()
 
 set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES})

qpms/materials.c

@@ -0,0 +1,142 @@
+#define _POSIX_C_SOURCE 200809L // for getline()
+#include <stdio.h>
+#include <stddef.h>
+#include <unistd.h>
+#include <string.h>
+#include "materials.h"
+#include "qpms_error.h"
+
+#define SPEED_OF_LIGHT (2.99792458e8)
+
+qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(
+    const size_t incount, const double *wavelen_m, const double *n, const double *k,
+    const gsl_interp_type *iptype)
+{
+  if (incount <= 0) return NULL;
+  qpms_permittivity_interpolator_t *ip;
+  QPMS_CRASHING_MALLOC(ip, sizeof(qpms_permittivity_interpolator_t));
+  ip->size = incount;
+  QPMS_CRASHING_MALLOC(ip->wavelength_m, incount * sizeof(double));
+  QPMS_CRASHING_MALLOC(ip->n, incount * sizeof(double));
+  QPMS_CRASHING_MALLOC(ip->k, incount * sizeof(double));
+  memcpy(ip->wavelength_m, wavelen_m, incount*sizeof(double));
+  memcpy(ip->k, k, incount*sizeof(double));
+  memcpy(ip->n, n, incount*sizeof(double));
+  ip->interp_n = gsl_interp_alloc(iptype, incount);
+  ip->interp_k = gsl_interp_alloc(iptype, incount);
+  QPMS_ENSURE_SUCCESS(gsl_interp_init(ip->interp_n, ip->wavelength_m, ip->n, incount));
+  QPMS_ENSURE_SUCCESS(gsl_interp_init(ip->interp_k, ip->wavelength_m, ip->k, incount));
+  return ip;
+}
+
+void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp)
+{
+  if(interp) {
+    gsl_interp_free(interp->interp_n);
+    gsl_interp_free(interp->interp_k);
+    free(interp->n);
+    free(interp->k);
+    free(interp->wavelength_m);
+  }
+  free(interp);
+}
+
+qpms_errno_t qpms_read_refractiveindex_yml(
+    FILE *f, ///< file handle
+    size_t *const count, ///< Number of successfully loaded triples.
+    double* *const lambdas_m, ///< Vacuum wavelengths in metres.
+    double* *const n, ///< Read refraction indices.
+    double* *const k ///< Read attenuation coeffs.
+    )
+{
+  QPMS_ENSURE(f && lambdas_m && n && k,"f, lambdas_m, n, k are mandatory arguments and must not be NULL.");
+  int count_alloc = 128; // First chunk to allocate
+  *count = 0;
+  QPMS_CRASHING_MALLOC(*lambdas_m, count_alloc * sizeof(double));
+  QPMS_CRASHING_MALLOC(*n, count_alloc * sizeof(double));
+  QPMS_CRASHING_MALLOC(*k, count_alloc * sizeof(double));
+  size_t linebufsz = 256;
+  char *linebuf;
+  QPMS_CRASHING_MALLOC(linebuf, linebufsz);
+  ssize_t readchars;
+  bool data_started = false;
+  while((readchars = getline(&linebuf, &linebufsz, f)) != -1) {
+    if (linebuf[0] == '#') continue;
+    // We need to find the beginning of the tabulated data; everything before that is ignored.
+    if (!data_started) {
+      char *test = strstr(linebuf, "data: |");
+      if(test) data_started = true;
+      continue;
+    }
+
+    if (3 == sscanf(linebuf, "%lf %lf %lf", *lambdas_m + *count, *n + *count , *k + *count)) {
+      (*lambdas_m)[*count] *= 1e-6; // The original data is in micrometres.
+      ++*count;
+      if (*count > count_alloc) {
+        count_alloc *= 2;
+        QPMS_CRASHING_REALLOC(*lambdas_m, count_alloc * sizeof(double));
+        QPMS_CRASHING_REALLOC(*n, count_alloc * sizeof(double));
+        QPMS_CRASHING_REALLOC(*k, count_alloc * sizeof(double));
+      }
+    } else break;
+  }
+  QPMS_ENSURE(*count > 0, "Could not read any refractive index data; the format must be wrong!");
+  free(linebuf);
+  return QPMS_SUCCESS;
+}
+
+qpms_errno_t qpms_load_refractiveindex_yml(
+    const char *path,
+    size_t *const count, ///< Number of successfully loaded triples.
+    double* *const lambdas_m, ///< Vacuum wavelengths in metres.
+    double* *const n, ///< Read refraction indices.
+    double* *const k ///< Read attenuation coeffs.
+    )
+{
+  FILE *f = fopen(path, "r");
+  QPMS_ENSURE(f, "Could not open refractive index file %s", path);
+  qpms_errno_t retval =
+    qpms_read_refractiveindex_yml(f, count, lambdas_m, n, k);
+  QPMS_ENSURE_SUCCESS(fclose(f));
+  return retval;
+}
+
+qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
+    const char *path, ///< Path to the yml file.
+    const gsl_interp_type *iptype ///< GSL interpolator type
+    )
+{
+  size_t count;
+  double *lambdas_m, *n, *k;
+  QPMS_ENSURE_SUCCESS(qpms_load_refractiveindex_yml(path, &count, &lambdas_m, &n, &k));
+  qpms_permittivity_interpolator_t *ip = qpms_permittivity_interpolator_create(
+      count, lambdas_m, n, k, iptype);
+  free(lambdas_m);
+  free(n);
+  free(k);
+  return ip;
+}
+
+complex double qpms_permittivity_interpolator_eps_at_omega(
+    const qpms_permittivity_interpolator_t *ip,  double omega_SI)
+{
+  double lambda, n, k;
+  lambda = 2*M_PI*SPEED_OF_LIGHT/omega_SI;
+  n = gsl_interp_eval(ip->interp_n, ip->wavelength_m, ip->n, lambda, NULL);
+  k = gsl_interp_eval(ip->interp_k, ip->wavelength_m, ip->k, lambda, NULL);
+  complex double epsilon = n*n - k*k + 2*n*k*I;
+  return epsilon;
+}
+
+double qpms_permittivity_interpolator_omega_max(
+    const qpms_permittivity_interpolator_t *ip)
+{
+  return 2*M_PI*SPEED_OF_LIGHT / ip->wavelength_m[0];
+}
+
+double qpms_permittivity_interpolator_omega_min(
+    const qpms_permittivity_interpolator_t *ip)
+{
+  return 2*M_PI*SPEED_OF_LIGHT / ip->wavelength_m[ip->size-1];
+}
+

qpms/materials.h

@@ -0,0 +1,61 @@
+/* \file materials.h
+ * \brief Optical properties of materials.
+ */
+#ifndef QPMS_MATERIALS_H
+#define QPMS_MATERIALS_H
+#include "qpms_types.h"
+#include <gsl/gsl_spline.h>
+
+/// Interpolator of tabulated optical properties.
+// TODO use gsl_interp instead of gsl_spline.
+typedef struct qpms_permittivity_interpolator_t {
+        double *wavelength_m; ///< Wavelength array (in meters).
+        double *n; ///< Refraction index array.
+        double *k; ///< Attenuation coefficient array.
+        gsl_interp *interp_n; ///< Refraction index interpolator object.
+        gsl_interp *interp_k; ///< Attenuation coeff interpolator object.
+        size_t size; ///< Size of n[], k[], and wavelength_m[].
+        // I could add gsl_interp_accel, but that is not necessary.
+} qpms_permittivity_interpolator_t;
+
+/// Creates a permittivity interpolator from tabulated wavelengths, refraction indices and extinction coeffs.
+qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(const size_t incount,
+                const double *wavelength_m, ///< Tabulated vacuum wavelength in metres, in strictly increasing order.
+                const double *n, ///< Tabulated refraction indices at omega.
+                const double *k, ///< Tabulated extinction coefficients.
+                const gsl_interp_type *iptype ///< GSL interpolator type
+                );
+
+/// Creates a permittivity interpolator from an yml file downloaded from refractiveindex.info website.
+qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
+                const char *path, ///< Path to the yml file.
+                const gsl_interp_type *iptype ///< GSL interpolator type
+                );
+
+/// Evaluates interpolated material permittivity at a given angular frequency.
+complex double qpms_permittivity_interpolator_eps_at_omega(
+                const qpms_permittivity_interpolator_t *interp, double omega_SI);
+
+/// Returns the minimum angular frequency supported by the interpolator.
+double qpms_permittivity_interpolator_omega_min(
+                const qpms_permittivity_interpolator_t *ip);
+
+/// Returns the minimum angular frequency supported by the interpolator.
+double qpms_permittivity_interpolator_omega_max(
+                const qpms_permittivity_interpolator_t *ip);
+
+/// Destroy a permittivity interpolator.
+void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp);
+
+/// Relative permittivity from the Drude model.
+static inline complex double qpms_drude_epsilon(
+                complex double eps_inf, ///< Relative permittivity "at infinity".
+                complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
+                complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
+                complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
+                ) {
+        return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p));
+}
+
+
+#endif //QPMS_MATERIALS_H

qpms/qpms_cdefs.pxd

@@ -366,6 +366,8 @@ cdef extern from "tmatrices.h":
             cdouble epsilon_fg, cdouble epsilon_bg)
     qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(const qpms_vswf_set_spec_t *bspec, double a,
             double omega, cdouble epsilon_fg, cdouble epsilon_bg)
+
+cdef extern from "materials.h":
     qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(const size_t incount,
             cdouble *wavelength_m, cdouble *n, cdouble *k, const gsl_interp_type *iptype)
     qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(const char *path,
@@ -375,6 +377,7 @@ cdef extern from "tmatrices.h":
     double qpms_permittivity_interpolator_omega_min(const qpms_permittivity_interpolator_t *interp)
     void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp)
 
+
 cdef extern from "pointgroups.h":
     bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls)
     double qpms_pg_quat_cmp_atol

qpms/qpms_types.h

@@ -399,6 +399,12 @@ typedef struct qpms_abstract_tmatrix_t {
 } qpms_abstract_tmatrix_t;
 
 
+/// A type holding electric permittivity and magnetic permeability of a material.
+typedef struct qpms_epsmu_t {
+	complex double eps; ///< Relative permittivity.
+	complex double mu; ///< Relative permeability.
+} qpms_epsmu_t;
+
 
 #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
 #endif // QPMS_TYPES

qpms/tmatrices.c

@@ -8,7 +8,6 @@
 #include "vswf.h"
 #include "groups.h"
 #include "symmetries.h"
-#include <gsl/gsl_spline.h>
 #include <assert.h>
 #include <unistd.h>
 #include "vectors.h"
@@ -557,138 +556,6 @@ qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose
   return QPMS_SUCCESS;
 }
 
-qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(
-    const size_t incount, const double *wavelen_m, const double *n, const double *k,
-    const gsl_interp_type *iptype) 
-{
-  if (incount <= 0) return NULL;
-  qpms_permittivity_interpolator_t *ip;
-  QPMS_CRASHING_MALLOC(ip, sizeof(qpms_permittivity_interpolator_t));
-  ip->size = incount;
-  QPMS_CRASHING_MALLOC(ip->wavelength_m, incount * sizeof(double));
-  QPMS_CRASHING_MALLOC(ip->n, incount * sizeof(double));
-  QPMS_CRASHING_MALLOC(ip->k, incount * sizeof(double));
-  memcpy(ip->wavelength_m, wavelen_m, incount*sizeof(double));
-  memcpy(ip->k, k, incount*sizeof(double));
-  memcpy(ip->n, n, incount*sizeof(double));
-  ip->interp_n = gsl_interp_alloc(iptype, incount);
-  ip->interp_k = gsl_interp_alloc(iptype, incount);
-  QPMS_ENSURE_SUCCESS(gsl_interp_init(ip->interp_n, ip->wavelength_m, ip->n, incount));
-  QPMS_ENSURE_SUCCESS(gsl_interp_init(ip->interp_k, ip->wavelength_m, ip->k, incount));
-  return ip;
-}
-
-void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp) 
-{
-  if(interp) {
-    gsl_interp_free(interp->interp_n);
-    gsl_interp_free(interp->interp_k);
-    free(interp->n);
-    free(interp->k);
-    free(interp->wavelength_m);
-  }
-  free(interp);
-}
-
-qpms_errno_t qpms_read_refractiveindex_yml(
-    FILE *f, ///< file handle
-    size_t *const count, ///< Number of successfully loaded triples.
-    double* *const lambdas_m, ///< Vacuum wavelengths in metres.
-    double* *const n, ///< Read refraction indices.
-    double* *const k ///< Read attenuation coeffs.
-    ) 
-{
-  QPMS_ENSURE(f && lambdas_m && n && k,"f, lambdas_m, n, k are mandatory arguments and must not be NULL.");
-  int count_alloc = 128; // First chunk to allocate
-  *count = 0;
-  QPMS_CRASHING_MALLOC(*lambdas_m, count_alloc * sizeof(double));
-  QPMS_CRASHING_MALLOC(*n, count_alloc * sizeof(double));
-  QPMS_CRASHING_MALLOC(*k, count_alloc * sizeof(double));
-  size_t linebufsz = 256;
-  char *linebuf;
-  QPMS_CRASHING_MALLOC(linebuf, linebufsz);
-  ssize_t readchars;
-  bool data_started = false;
-  while((readchars = getline(&linebuf, &linebufsz, f)) != -1) {
-    if (linebuf[0] == '#') continue;
-    // We need to find the beginning of the tabulated data; everything before that is ignored.
-    if (!data_started) {
-      char *test = strstr(linebuf, "data: |");
-      if(test) data_started = true;
-      continue;
-    }
-
-    if (3 == sscanf(linebuf, "%lf %lf %lf", *lambdas_m + *count, *n + *count , *k + *count)) {
-      (*lambdas_m)[*count] *= 1e-6; // The original data is in micrometres.
-      ++*count;
-      if (*count > count_alloc) {
-        count_alloc *= 2;
-        QPMS_CRASHING_REALLOC(*lambdas_m, count_alloc * sizeof(double));
-        QPMS_CRASHING_REALLOC(*n, count_alloc * sizeof(double));
-        QPMS_CRASHING_REALLOC(*k, count_alloc * sizeof(double));
-      }
-    } else break;
-  }
-  QPMS_ENSURE(*count > 0, "Could not read any refractive index data; the format must be wrong!");
-  free(linebuf);
-  return QPMS_SUCCESS; 
-}
-
-qpms_errno_t qpms_load_refractiveindex_yml(
-    const char *path,
-    size_t *const count, ///< Number of successfully loaded triples.
-    double* *const lambdas_m, ///< Vacuum wavelengths in metres.
-    double* *const n, ///< Read refraction indices.
-    double* *const k ///< Read attenuation coeffs.
-    ) 
-{
-  FILE *f = fopen(path, "r");
-  QPMS_ENSURE(f, "Could not open refractive index file %s", path); 
-  qpms_errno_t retval = 
-    qpms_read_refractiveindex_yml(f, count, lambdas_m, n, k);
-  QPMS_ENSURE_SUCCESS(fclose(f)); 
-  return retval;
-}
-
-qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
-		const char *path, ///< Path to the yml file.
-		const gsl_interp_type *iptype ///< GSL interpolator type
-		) 
-{
-  size_t count;
-  double *lambdas_m, *n, *k;
-  QPMS_ENSURE_SUCCESS(qpms_load_refractiveindex_yml(path, &count, &lambdas_m, &n, &k));
-  qpms_permittivity_interpolator_t *ip = qpms_permittivity_interpolator_create(
-      count, lambdas_m, n, k, iptype);
-  free(lambdas_m);
-  free(n);
-  free(k);
-  return ip;
-}
-
-complex double qpms_permittivity_interpolator_eps_at_omega(
-    const qpms_permittivity_interpolator_t *ip,  double omega_SI) 
-{ 
-  double lambda, n, k;
-  lambda = 2*M_PI*SPEED_OF_LIGHT/omega_SI;
-  n = gsl_interp_eval(ip->interp_n, ip->wavelength_m, ip->n, lambda, NULL);
-  k = gsl_interp_eval(ip->interp_k, ip->wavelength_m, ip->k, lambda, NULL);
-  complex double epsilon = n*n - k*k + 2*n*k*I;
-  return epsilon;
-}
-
-double qpms_permittivity_interpolator_omega_max(
-    const qpms_permittivity_interpolator_t *ip)
-{
-  return 2*M_PI*SPEED_OF_LIGHT / ip->wavelength_m[0];
-}
-
-double qpms_permittivity_interpolator_omega_min(
-    const qpms_permittivity_interpolator_t *ip)
-{
-  return 2*M_PI*SPEED_OF_LIGHT / ip->wavelength_m[ip->size-1];
-}
-
 /// Convenience function to calculate T-matrix of a non-magnetic spherical \
 particle using the permittivity values, replacing existing T-matrix data.
 qpms_errno_t qpms_tmatrix_spherical_mu0_fill(

qpms/tmatrices.h

@@ -3,8 +3,9 @@
  */
 #ifndef TMATRICES_H
 #define TMATRICES_H
-#include "qpms_types.h"
+// #include "qpms_types.h" // included via materials.h
-#include <gsl/gsl_spline.h>
+// #include <gsl/gsl_spline.h> // included via materials.h
+#include "materials.h"
 #include <stdio.h>
 
 struct qpms_finite_group_t;
@@ -276,48 +277,6 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
 	       //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
        	       /*, ...? */);
 
-
-/// Interpolator of tabulated optical properties.
-// TODO use gsl_interp instead of gsl_spline.
-typedef struct qpms_permittivity_interpolator_t {
-	double *wavelength_m; ///< Wavelength array (in meters).
-	double *n; ///< Refraction index array.
-	double *k; ///< Attenuation coefficient array.
-	gsl_interp *interp_n; ///< Refraction index interpolator object.
-	gsl_interp *interp_k; ///< Attenuation coeff interpolator object.
-	size_t size; ///< Size of n[], k[], and wavelength_m[].
-	// I could add gsl_interp_accel, but that is not necessary.
-} qpms_permittivity_interpolator_t;
-
-/// Creates a permittivity interpolator from tabulated wavelengths, refraction indices and extinction coeffs.
-qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(const size_t incount,
-		const double *wavelength_m, ///< Tabulated vacuum wavelength in metres, in strictly increasing order. 
-		const double *n, ///< Tabulated refraction indices at omega.
-	       	const double *k, ///< Tabulated extinction coefficients.
-		const gsl_interp_type *iptype ///< GSL interpolator type
-		);
-
-/// Creates a permittivity interpolator from an yml file downloaded from refractiveindex.info website.
-qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
-		const char *path, ///< Path to the yml file.
-		const gsl_interp_type *iptype ///< GSL interpolator type
-		);
-
-/// Evaluates interpolated material permittivity at a given angular frequency.
-complex double qpms_permittivity_interpolator_eps_at_omega(
-		const qpms_permittivity_interpolator_t *interp, double omega_SI);
-
-/// Returns the minimum angular frequency supported by the interpolator.
-double qpms_permittivity_interpolator_omega_min(
-		const qpms_permittivity_interpolator_t *ip);
-
-/// Returns the minimum angular frequency supported by the interpolator.
-double qpms_permittivity_interpolator_omega_max(
-		const qpms_permittivity_interpolator_t *ip);
-
-/// Destroy a permittivity interpolator.
-void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp);
-
 /// Calculates the reflection Mie-Lorentz coefficients for a spherical particle.
 /**
  * This function is based on the previous python implementation mie_coefficients() from qpms_p.py,
@@ -366,16 +325,6 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical(
 	return t;
 }
 
-/// Relative permittivity from the Drude model.
-static inline complex double qpms_drude_epsilon(
-		complex double eps_inf, ///< Relative permittivity "at infinity".
-		complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
-		complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
-		complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
-		) {
-	return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p));
-}
-
 /// Convenience function to calculate T-matrix of a non-magnetic spherical \
 particle using the permittivity values, replacing existing T-matrix data.
 qpms_errno_t qpms_tmatrix_spherical_mu0_fill(
@@ -399,8 +348,6 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(
 };
 
 
-
-
 #if 0
 // Abstract types that describe T-matrix/particle/scatsystem symmetries
 // To be implemented later. See also the thoughts in the beginning of groups.h.

setup.py

@@ -72,6 +72,7 @@ qpms_c = Extension('qpms_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/tmatrices.c',
+            'qpms/materials.c',
             'qpms/error.c',
             'qpms/bessel.c',
             'qpms/own_zgemm.c',