solstice-solver

ssol_spectrum.c (6072B)

---

 /* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com)
  * Copyright (C) 2016, 2018 CNRS
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include "ssol.h"
 #include "ssol_spectrum_c.h"
 #include "ssol_device_c.h"
 
 #include <rsys/algorithm.h>
 #include <rsys/hash.h>
 #include <rsys/math.h>
 #include <rsys/mem_allocator.h>
 #include <rsys/ref_count.h>
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static void
 spectrum_release(ref_T* ref)
 {
   struct ssol_device* dev;
   struct ssol_spectrum* spectrum = CONTAINER_OF(ref, struct ssol_spectrum, ref);
   ASSERT(ref);
   dev = spectrum->dev;
   ASSERT(dev && dev->allocator);
   darray_double_release(&spectrum->wavelengths);
   darray_double_release(&spectrum->intensities);
   MEM_RM(dev->allocator, spectrum);
   SSOL(device_ref_put(dev));
 }
 
 static int
 eq_dbl(const void* key, const void* base)
 {
   const double k = *(const double*) key;
   const double b = *(const double*) base;
   if(k > b) return +1;
   if(k < b) return -1;
   return 0;
 }
 
 /*******************************************************************************
  * Local ssol_spectrum functions
  ******************************************************************************/
 double
 spectrum_interpolate
   (const struct ssol_spectrum* spectrum, const double wavelength)
 {
   const double* wls;
   const double* ints;
   const double* next;
   double slope;
   double intensity;
   size_t id_next, sz;
   ASSERT(spectrum);
 
   sz = darray_double_size_get(&spectrum->wavelengths);
   wls = darray_double_cdata_get(&spectrum->wavelengths);
   ints = darray_double_cdata_get(&spectrum->intensities);
   next = search_lower_bound(&wavelength, wls, sz, sizeof(double), &eq_dbl);
   if(!next) { /* Clamp to upper bound */
     return ints[sz-1];
   }
 
   id_next = (size_t)(next - wls);
   if(!id_next) { /* Clamp to lower bound */
     return ints[0];
   }
 
   ASSERT(id_next); /* because spectrum_includes_point */
   ASSERT(wls[id_next] >= wls[id_next - 1]);
 
   slope = (ints[id_next] - ints[id_next-1]) / (wls[id_next] - wls[id_next-1]);
   intensity = ints[id_next-1] + (wavelength - wls[id_next - 1]) * slope;
   ASSERT(intensity >= 0);
   return intensity;
 }
 
 int
 spectrum_check_data
   (const struct ssol_spectrum* spectrum, const double lower, const double upper)
 {
   size_t sz, i;
   double current_wl = 0;
   ASSERT(spectrum && lower <= upper);
   sz = darray_double_size_get(&spectrum->intensities);
   if(!sz) return 0;
   if(sz != darray_double_size_get(&spectrum->wavelengths)) return 0;
   FOR_EACH(i, 0, sz) {
     const double wl = darray_double_cdata_get(&spectrum->wavelengths)[i];
     const double data = darray_double_cdata_get(&spectrum->intensities)[i];
     if(data < lower || data > upper) return 0;
     if(wl <= 0) return 0;
     if(wl <= current_wl) return 0;
     current_wl = wl;
   }
   return 1;
 }
 
 /*******************************************************************************
  * Exported ssol_spectrum functions
  ******************************************************************************/
 res_T
 ssol_spectrum_create
   (struct ssol_device* dev, struct ssol_spectrum** out_spectrum)
 {
   struct ssol_spectrum* spectrum = NULL;
   res_T res = RES_OK;
 
   if(!dev || !out_spectrum) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   spectrum = MEM_CALLOC(dev->allocator, 1, sizeof(struct ssol_spectrum));
   if(!spectrum) {
     res = RES_MEM_ERR;
     goto error;
   }
 
   SSOL(device_ref_get(dev));
   spectrum->dev = dev;
   ref_init(&spectrum->ref);
   darray_double_init(dev->allocator, &spectrum->wavelengths);
   darray_double_init(dev->allocator, &spectrum->intensities);
 
 exit:
   if(out_spectrum) *out_spectrum = spectrum;
   return res;
 error:
   if(spectrum) {
     SSOL(spectrum_ref_put(spectrum));
     spectrum = NULL;
   }
   goto exit;
 }
 
 res_T
 ssol_spectrum_ref_get(struct ssol_spectrum* spectrum)
 {
   if(!spectrum) return RES_BAD_ARG;
   ref_get(&spectrum->ref);
   return RES_OK;
 }
 
 res_T
 ssol_spectrum_ref_put(struct ssol_spectrum* spectrum)
 {
   if(!spectrum) return RES_BAD_ARG;
   ref_put(&spectrum->ref, spectrum_release);
   return RES_OK;
 }
 
 SSOL_API res_T
 ssol_spectrum_setup
   (struct ssol_spectrum* spectrum,
    void (*get)(const size_t iwlen, double* wlen, double* data, void* ctx),
    const size_t nwlens,
    void* ctx)
 {
   double* wavelengths;
   double* intensities;
   double current_wl = 0;
   size_t i;
   res_T res = RES_OK;
 
   if(!spectrum || !nwlens || !get) {
     res = RES_BAD_ARG;
     goto error;
   }
 
   res = darray_double_resize(&spectrum->wavelengths, nwlens);
   if(res != RES_OK) goto error;
   res = darray_double_resize(&spectrum->intensities, nwlens);
   if(res != RES_OK) goto error;
 
   wavelengths = darray_double_data_get(&spectrum->wavelengths);
   intensities = darray_double_data_get(&spectrum->intensities);
   FOR_EACH(i, 0, nwlens) {
     get(i, wavelengths + i, intensities + i, ctx);
     if(wavelengths[i] <= current_wl || intensities[i] < 0) {
       res = RES_BAD_ARG;
       goto error;
     }
     current_wl = *(wavelengths + i);
   }
 
   spectrum->checksum[0] = hash_fnv64(wavelengths, nwlens*sizeof(double));
   spectrum->checksum[1] = hash_fnv64(intensities, nwlens*sizeof(double));
 
 exit:
   return res;
 error:
   if(spectrum) {
     darray_double_clear(&spectrum->wavelengths);
     darray_double_clear(&spectrum->intensities);
     spectrum->checksum[0] = 0;
     spectrum->checksum[1] = 0;
   }
   goto exit;
 }