qpms/materials.c

   1 #define _POSIX_C_SOURCE 200809L // for getline()
   2 #include <stdio.h>
   3 #include <stddef.h>
   4 #include <unistd.h>
   5 #include <string.h>
   6 #include "materials.h"
   7 #include "qpms_error.h"
   8
   9 #define SQ(x) ((x)*(x))
  10
  11 qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(
  12     const size_t incount, const double *wavelen_m, const double *n, const double *k,
  13     const gsl_interp_type *iptype)
  14 {
  15   if (incount <= 0) return NULL;
  16   qpms_permittivity_interpolator_t *ip;
  17   QPMS_CRASHING_MALLOC(ip, sizeof(qpms_permittivity_interpolator_t));
  18   ip->size = incount;
  19   QPMS_CRASHING_MALLOC(ip->wavelength_m, incount * sizeof(double));
  20   QPMS_CRASHING_MALLOC(ip->n, incount * sizeof(double));
  21   QPMS_CRASHING_MALLOC(ip->k, incount * sizeof(double));
  22   memcpy(ip->wavelength_m, wavelen_m, incount*sizeof(double));
  23   memcpy(ip->k, k, incount*sizeof(double));
  24   memcpy(ip->n, n, incount*sizeof(double));
  25   ip->interp_n = gsl_interp_alloc(iptype, incount);
  26   ip->interp_k = gsl_interp_alloc(iptype, incount);
  27   QPMS_ENSURE_SUCCESS(gsl_interp_init(ip->interp_n, ip->wavelength_m, ip->n, incount));
  28   QPMS_ENSURE_SUCCESS(gsl_interp_init(ip->interp_k, ip->wavelength_m, ip->k, incount));
  29   return ip;
  30 }
  31
  32 void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp)
  33 {
  34   if(interp) {
  35     gsl_interp_free(interp->interp_n);
  36     gsl_interp_free(interp->interp_k);
  37     free(interp->n);
  38     free(interp->k);
  39     free(interp->wavelength_m);
  40   }
  41   free(interp);
  42 }
  43
  44 qpms_errno_t qpms_read_refractiveindex_yml(
  45     FILE *f, ///< file handle
  46     size_t *const count, ///< Number of successfully loaded triples.
  47     double* *const lambdas_m, ///< Vacuum wavelengths in metres.
  48     double* *const n, ///< Read refraction indices.
  49     double* *const k ///< Read attenuation coeffs.
  50     )
  51 {
  52   QPMS_ENSURE(f && lambdas_m && n && k,"f, lambdas_m, n, k are mandatory arguments and must not be NULL.");
  53   int count_alloc = 128; // First chunk to allocate
  54   *count = 0;
  55   QPMS_CRASHING_MALLOC(*lambdas_m, count_alloc * sizeof(double));
  56   QPMS_CRASHING_MALLOC(*n, count_alloc * sizeof(double));
  57   QPMS_CRASHING_MALLOC(*k, count_alloc * sizeof(double));
  58   size_t linebufsz = 256;
  59   char *linebuf;
  60   QPMS_CRASHING_MALLOC(linebuf, linebufsz);
  61   ssize_t readchars;
  62   bool data_started = false;
  63   while((readchars = getline(&linebuf, &linebufsz, f)) != -1) {
  64     if (linebuf[0] == '#') continue;
  65     // We need to find the beginning of the tabulated data; everything before that is ignored.
  66     if (!data_started) {
  67       char *test = strstr(linebuf, "data: |");
  68       if(test) data_started = true;
  69       continue;
  70     }
  71
  72     if (3 == sscanf(linebuf, "%lf %lf %lf", *lambdas_m + *count, *n + *count , *k + *count)) {
  73       (*lambdas_m)[*count] *= 1e-6; // The original data is in micrometres.
  74       ++*count;
  75       if (*count > count_alloc) {
  76         count_alloc *= 2;
  77         QPMS_CRASHING_REALLOC(*lambdas_m, count_alloc * sizeof(double));
  78         QPMS_CRASHING_REALLOC(*n, count_alloc * sizeof(double));
  79         QPMS_CRASHING_REALLOC(*k, count_alloc * sizeof(double));
  80       }
  81     } else break;
  82   }
  83   QPMS_ENSURE(*count > 0, "Could not read any refractive index data; the format must be wrong!");
  84   free(linebuf);
  85   return QPMS_SUCCESS;
  86 }
  87
  88 qpms_errno_t qpms_load_refractiveindex_yml(
  89     const char *path,
  90     size_t *const count, ///< Number of successfully loaded triples.
  91     double* *const lambdas_m, ///< Vacuum wavelengths in metres.
  92     double* *const n, ///< Read refraction indices.
  93     double* *const k ///< Read attenuation coeffs.
  94     )
  95 {
  96   FILE *f = fopen(path, "r");
  97   QPMS_ENSURE(f, "Could not open refractive index file %s", path);
  98   qpms_errno_t retval =
  99     qpms_read_refractiveindex_yml(f, count, lambdas_m, n, k);
 100   QPMS_ENSURE_SUCCESS(fclose(f));
 101   return retval;
 102 }
 103
 104 qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
 105     const char *path, ///< Path to the yml file.
 106     const gsl_interp_type *iptype ///< GSL interpolator type
 107     )
 108 {
 109   size_t count;
 110   double *lambdas_m, *n, *k;
 111   QPMS_ENSURE_SUCCESS(qpms_load_refractiveindex_yml(path, &count, &lambdas_m, &n, &k));
 112   qpms_permittivity_interpolator_t *ip = qpms_permittivity_interpolator_create(
 113       count, lambdas_m, n, k, iptype);
 114   free(lambdas_m);
 115   free(n);
 116   free(k);
 117   return ip;
 118 }
 119
 120 complex double qpms_permittivity_interpolator_eps_at_omega(
 121     const qpms_permittivity_interpolator_t *ip,  double omega_SI)
 122 {
 123   double lambda, n, k;
 124   lambda = 2*M_PI*SPEED_OF_LIGHT/omega_SI;
 125   n = gsl_interp_eval(ip->interp_n, ip->wavelength_m, ip->n, lambda, NULL);
 126   k = gsl_interp_eval(ip->interp_k, ip->wavelength_m, ip->k, lambda, NULL);
 127   complex double epsilon = n*n - k*k + 2*n*k*I;
 128   return epsilon;
 129 }
 130
 131 qpms_epsmu_t qpms_permittivity_interpolator_epsmu_g(
 132     complex double omega, const void *p)
 133 {
 134   const qpms_permittivity_interpolator_t *interp = p;
 135   static bool imag_already_bitched = false;
 136   if(cimag(omega) && !imag_already_bitched)
 137     QPMS_WARN("Complex frequencies not supported by qpms_permittivity_interpolator_t. Imaginary parts will be discarded!");
 138   qpms_epsmu_t em;
 139   em.eps = qpms_permittivity_interpolator_eps_at_omega(interp, omega);
 140   em.mu = 1;
 141   return em;
 142 }
 143
 144 double qpms_permittivity_interpolator_omega_max(
 145     const qpms_permittivity_interpolator_t *ip)
 146 {
 147   return 2*M_PI*SPEED_OF_LIGHT / ip->wavelength_m[0];
 148 }
 149
 150 double qpms_permittivity_interpolator_omega_min(
 151     const qpms_permittivity_interpolator_t *ip)
 152 {
 153   return 2*M_PI*SPEED_OF_LIGHT / ip->wavelength_m[ip->size-1];
 154 }
 155
 156 complex double qpms_lorentzdrude_eps(complex double omega, const qpms_ldparams_t *p)
 157 {
 158   complex double eps = 0;
 159   for(size_t j = 0; j < p->n; ++j) {
 160     const qpms_ldparams_triple_t d = p->data[j];
 161     eps += d.f * SQ(p->omega_p) / (SQ(d.omega) - SQ(omega) - I*omega*d.gamma );
 162   }
 163   return eps;
 164 }
 165
 166 qpms_epsmu_t qpms_lorentzdrude_epsmu(complex double omega, const qpms_ldparams_t *p)
 167 {
 168   qpms_epsmu_t em;
 169   em.eps = qpms_lorentzdrude_eps(omega, p);
 170   em.mu = 1;
 171   return em;
 172 }
 173
 174 qpms_epsmu_t qpms_lorentzdrude_epsmu_g(complex double omega, const void *p)
 175 {
 176   return qpms_lorentzdrude_epsmu(omega, (const qpms_ldparams_t *)p);
 177 }
 178
 179 qpms_epsmu_t qpms_epsmu_const_g(complex double omega, const void *p)
 180 {
 181   return *(const qpms_epsmu_t *)p;
 182 }