diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 4e18059..259eafd 100644 --- a/qpms/CMakeLists.txt +++ b/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}) diff --git a/qpms/materials.c b/qpms/materials.c new file mode 100644 index 0000000..0be5467 --- /dev/null +++ b/qpms/materials.c @@ -0,0 +1,142 @@ +#define _POSIX_C_SOURCE 200809L // for getline() +#include +#include +#include +#include +#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]; +} + diff --git a/qpms/materials.h b/qpms/materials.h new file mode 100644 index 0000000..55cc012 --- /dev/null +++ b/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 + +/// 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 diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 682f92f..9e9acbf 100644 --- a/qpms/qpms_cdefs.pxd +++ b/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 diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 3a99f12..4d7955f 100644 --- a/qpms/qpms_types.h +++ b/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 diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index b58767b..8925193 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -8,7 +8,6 @@ #include "vswf.h" #include "groups.h" #include "symmetries.h" -#include #include #include #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( diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index fb2b42b..e0a9bf9 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -3,8 +3,9 @@ */ #ifndef TMATRICES_H #define TMATRICES_H -#include "qpms_types.h" -#include +// #include "qpms_types.h" // included via materials.h +// #include // included via materials.h +#include "materials.h" #include 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, @@ -332,16 +291,16 @@ void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *inter * TODO better doc. */ complex double *qpms_mie_coefficients_reflection( - complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated. - const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated. - double a, ///< Radius of the sphere. - complex double k_i, ///< Wave number of the internal material of the sphere. - complex double k_e, ///< Wave number of the surrounding medium. - complex double mu_i, ///< Relative permeability of the sphere material. - complex double mu_e, ///< Relative permeability of the surrounding medium. - qpms_bessel_t J_ext, ///< Kind of the "incoming" waves. Most likely QPMS_BESSEL_REGULAR. - qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS. - ); + complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated. + const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated. + double a, ///< Radius of the sphere. + complex double k_i, ///< Wave number of the internal material of the sphere. + complex double k_e, ///< Wave number of the surrounding medium. + complex double mu_i, ///< Relative permeability of the sphere material. + complex double mu_e, ///< Relative permeability of the surrounding medium. + qpms_bessel_t J_ext, ///< Kind of the "incoming" waves. Most likely QPMS_BESSEL_REGULAR. + qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS. + ); /// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory. qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL. @@ -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. diff --git a/setup.py b/setup.py index 40b0a87..1e3ad14 100755 --- a/setup.py +++ b/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',