star-phor-input

commit 47828fa38176281f606a1419e457b8cfc306350b
parent 283f82b38a7b3c867accf2f0a3a1b2e2132ebd11
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Mon,  8 Sep 2025 14:52:44 +0200

Normalize source emission spectrum during parsing

Ensure that source spectra are normalized so the integral between lambda
min and lambda max equals 1.

Normalization is performed using the trapezoidal method.

Diffstat:
M
src/sphin_source_surface.c
 | 
40
+++++++++++++++++++++++++++++++++++++---
M
src/sphin_spectral_property.c
 | 
8
++++----
M
src/test_sphin_load_source.c
 | 
6
+++---

3 files changed, 44 insertions(+), 10 deletions(-)
diff --git a/src/sphin_source_surface.c b/src/sphin_source_surface.c
@@ -294,6 +294,37 @@ error:
 }
 
 static res_T
+normalize_spectrum
+  (struct sphin_spectral_property* spectrum)
+{
+  double sum = 0;
+  double* wavelengths = NULL;
+  double* values = NULL;
+  size_t data_count = 0;
+  size_t i = 0;
+  res_T res = RES_OK;
+
+  wavelengths = darray_double_data_get(&spectrum->wavelengths);
+  values = darray_double_data_get(&spectrum->values);
+  data_count = darray_double_size_get(&spectrum->wavelengths);
+
+  /* Compute the integral of the spectrum over the wavelengths domain using
+   * the trapeze method */
+
+  FOR_EACH(i, 0, data_count-1) {
+    double dx = wavelengths[i+1] - wavelengths[i]; /* Interval size */
+    sum += 0.5 * dx * (values[i] + values[i+1]);
+  }
+
+  /* Normalize spectrum */
+  FOR_EACH(i, 0, data_count-1) {
+    values[i] = values[i] / sum;
+  }
+
+  return res;
+}
+
+static res_T
 parse_flux_density
   (struct sphin_source_surface* source,
    struct txtrdr* txtrdr,
@@ -354,11 +385,14 @@ parse_flux_density
       (source->sphin, filename, &source->flux_density.emission_spectrum);
     if (RES_OK != res) { goto error; }
 
-    /* wavelengths will be in nm after conversion */
-    res = convert_wavelengths(&source->flux_density.emission_spectrum->wavelengths, spectral_unit); 
+    /* Wavelengths will be in nm after conversion */
+    res = convert_wavelengths
+      (&source->flux_density.emission_spectrum->wavelengths, spectral_unit);
     if (RES_OK != res) { goto error; }
 
-    /* TODO Convert and property according to units */
+    /* Normalize spectrum */
+    res = normalize_spectrum(source->flux_density.emission_spectrum);
+    if (RES_OK != res) { goto error; }
   }
 
   res = txtrdr_read_line(txtrdr);
diff --git a/src/sphin_spectral_property.c b/src/sphin_spectral_property.c
@@ -223,7 +223,7 @@ convert_wavelengths
   double* wls = NULL;
   size_t i, wl_count = 0;
   res_T res = RES_OK;
-  
+
   ASSERT(NULL != wavelengths);
   ASSERT(NULL != spectral_unit);
 
@@ -231,10 +231,10 @@ convert_wavelengths
   else if (0 == strcmp(spectral_unit, "cm")) {scaling_factor = 1e07;}
   else if (0 == strcmp(spectral_unit, "m")) {scaling_factor = 1e09;}
   else if (0 == strcmp(spectral_unit, "cm^-1")) {
-    scaling_factor = 1e07; 
+    scaling_factor = 1e07;
     invert = 1;}
   else if (0 == strcmp(spectral_unit, "1/cm")) {
-    scaling_factor = 1e07; 
+    scaling_factor = 1e07;
     invert = 1;}
   else { res = RES_BAD_ARG; goto error; }
 
@@ -246,7 +246,7 @@ convert_wavelengths
       if (0.0 == wls[i]) { res = RES_BAD_ARG; goto error; }
       wls[i] = scaling_factor / wls[i];
     }
-    else {wls[i] = wls[i] * scaling_factor;} 
+    else {wls[i] = wls[i] * scaling_factor;}
   }
 
 exit:
diff --git a/src/test_sphin_load_source.c b/src/test_sphin_load_source.c
@@ -151,12 +151,12 @@ test_source_api
   CHK(sphin_spectral_property_get_desc(flux_density.emission_spectrum, &emission_desc) == RES_OK);
   CHK(sphin_spectral_property_interpolate_at_wavelength
       (flux_density.emission_spectrum, 100, SPHIN_INTERPOLATION_LINEAR, &value) == RES_OK);
-  CHK(eq_eps(value, 1., 1e-15));
+  CHK(eq_eps(value, 1./750, 1e-15));
   CHK(sphin_spectral_property_interpolate_at_wavelength
       (flux_density.emission_spectrum, 175, SPHIN_INTERPOLATION_LINEAR, &value) == RES_OK);
-  CHK(eq_eps(value, 1.75, 1e-15));
+  CHK(eq_eps(value, 1.75/750, 1e-15));
   CHK(eq_eps(emission_desc.wavelengths[2], 300, 1e-15));
-  CHK(eq_eps(emission_desc.values[2], 3, 1e-15));
+  CHK(eq_eps(emission_desc.values[2], 3./750, 1e-15));
   CHK(emission_desc.data_count == 4);
 
   CHK(sphin_config_get_surface(config, 2, &surface) == RES_OK);