QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Reimplementation of qpms_permittivity_interpolator_t. 

Former-commit-id: 740c89d441f3970bdd2aa4d64d8707c74e990b0c
This commit is contained in:
Marek Nečada 2019-03-24 14:28:09 +00:00
parent 115d1d07fe
commit d7a64f5928
2 changed files with 48 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
qpms/tmatrices.c
View File

@ -547,18 +547,28 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(
if (incount <= 0) return NULL; if (incount <= 0) return NULL;
qpms_permittivity_interpolator_t *ip; qpms_permittivity_interpolator_t *ip;
QPMS_CRASHING_MALLOC(ip, sizeof(qpms_permittivity_interpolator_t)); QPMS_CRASHING_MALLOC(ip, sizeof(qpms_permittivity_interpolator_t));
ip->spline_n = gsl_spline_alloc(iptype, incount); ip->size = incount;
ip->spline_k = gsl_spline_alloc(iptype, incount); QPMS_CRASHING_MALLOC(ip->wavelength_m, incount * sizeof(double));
QPMS_ENSURE_SUCCESS(gsl_spline_init(ip->spline_n, wavelen_m, n, incount)); QPMS_CRASHING_MALLOC(ip->n, incount * sizeof(double));
QPMS_ENSURE_SUCCESS(gsl_spline_init(ip->spline_k, wavelen_m, k, incount)); 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; return ip;
} }
void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp) void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp)
{ {
if(interp) { if(interp) {
gsl_spline_free(interp->spline_n); gsl_interp_free(interp->interp_n);
gsl_spline_free(interp->spline_k); gsl_interp_free(interp->interp_k);
free(interp->n);
free(interp->k);
free(interp->wavelength_m);
} }
free(interp); free(interp);
} }
@ -641,15 +651,26 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
complex double qpms_permittivity_interpolator_eps_at_omega( complex double qpms_permittivity_interpolator_eps_at_omega(
const qpms_permittivity_interpolator_t *ip, double omega_SI) const qpms_permittivity_interpolator_t *ip, double omega_SI)
{ {
double lambda, n, k; double lambda, n, k;
lambda = 2*M_PI*SPEED_OF_LIGHT/omega_SI; lambda = 2*M_PI*SPEED_OF_LIGHT/omega_SI;
n = gsl_spline_eval(ip->spline_n, lambda, NULL); n = gsl_interp_eval(ip->interp_n, ip->wavelength_m, ip->n, lambda, NULL);
k = gsl_spline_eval(ip->spline_k, 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; complex double epsilon = n*n - k*k + 2*n*k*I;
return epsilon; 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 \ /// Convenience function to calculate T-matrix of a non-magnetic spherical \
particle using the permittivity values, replacing existing T-matrix data. particle using the permittivity values, replacing existing T-matrix data.

21
qpms/tmatrices.h
View File

@ -252,9 +252,15 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
/// Interpolator of tabulated optical properties. /// Interpolator of tabulated optical properties.
// TODO use gsl_interp instead of gsl_spline.
typedef struct qpms_permittivity_interpolator_t { typedef struct qpms_permittivity_interpolator_t {
gsl_spline *spline_n; double *wavelength_m; ///< Wavelength array (in meters).
gsl_spline *spline_k; 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; } qpms_permittivity_interpolator_t;
/// Creates a permittivity interpolator from tabulated wavelengths, refraction indices and extinction coeffs. /// 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. /// 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. /// Destroy a permittivity interpolator.
void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp); void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp);