From d7a64f59287c592bb66392e02d4466cb06f34ba7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 24 Mar 2019 14:28:09 +0000 Subject: [PATCH] Reimplementation of qpms_permittivity_interpolator_t. Former-commit-id: 740c89d441f3970bdd2aa4d64d8707c74e990b0c --- qpms/tmatrices.c | 39 ++++++++++++++++++++++++++++++--------- qpms/tmatrices.h | 21 ++++++++++++++++++--- 2 files changed, 48 insertions(+), 12 deletions(-) diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 4946146..9ee21f1 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -547,18 +547,28 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create( if (incount <= 0) return NULL; qpms_permittivity_interpolator_t *ip; QPMS_CRASHING_MALLOC(ip, sizeof(qpms_permittivity_interpolator_t)); - ip->spline_n = gsl_spline_alloc(iptype, incount); - ip->spline_k = gsl_spline_alloc(iptype, incount); - QPMS_ENSURE_SUCCESS(gsl_spline_init(ip->spline_n, wavelen_m, n, incount)); - QPMS_ENSURE_SUCCESS(gsl_spline_init(ip->spline_k, wavelen_m, k, incount)); + 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_spline_free(interp->spline_n); - gsl_spline_free(interp->spline_k); + gsl_interp_free(interp->interp_n); + gsl_interp_free(interp->interp_k); + free(interp->n); + free(interp->k); + free(interp->wavelength_m); } free(interp); } @@ -641,15 +651,26 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml( 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_spline_eval(ip->spline_n, lambda, NULL); - k = gsl_spline_eval(ip->spline_k, lambda, NULL); + 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. diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index c0e114a..ae1208e 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -252,9 +252,15 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num /// Interpolator of tabulated optical properties. +// TODO use gsl_interp instead of gsl_spline. typedef struct qpms_permittivity_interpolator_t { - gsl_spline *spline_n; - gsl_spline *spline_k; + 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. @@ -272,7 +278,16 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml( ); /// 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); +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);