star-line

sln_line.c (18022B)

---

 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2026 Université de Lorraine
  * Copyright (C) 2022 Centre National de la Recherche Scientifique
  * Copyright (C) 2022 Université Paul Sabatier
  *
  * 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/>. */
 
 #define _POSIX_C_SOURCE 200112L /* nextafterf support */
 
 #include "sln_device_c.h"
 #include "sln_line.h"
 #include "sln_tree_c.h"
 
 #include <lblu.h>
 #include <star/shtr.h>
 
 #include <rsys/algorithm.h>
 #include <rsys/cstr.h>
 #include <rsys/dynamic_array_double.h>
 #include <rsys/math.h>
 
 #include <math.h> /* nextafterf */
 
 #define T_REF 296.0 /* K */
 #define AVOGADRO_NUMBER 6.02214076e23 /* molec.mol^-1 */
 #define PERFECT_GAZ_CONSTANT 8.2057e-5 /* m^3.atm.mol^-1.K^-1 */
 
 #define MIN_NVERTICES_HINT 8
 #define MAX_NVERTICES_HINT 128
 STATIC_ASSERT(IS_POW2(MIN_NVERTICES_HINT), MIN_NVERTICES_HINT_is_not_a_pow2);
 STATIC_ASSERT(IS_POW2(MIN_NVERTICES_HINT), MAX_NVERTICES_HINT_is_not_a_pow2);
 
 /*******************************************************************************
  * Helper function
  ******************************************************************************/
 static INLINE double
 line_intensity
   (const double intensity_ref, /* Reference intensity [cm^-1/(molec.cm^2)] */
    const double lower_state_energy, /* [cm^-1] */
    const double partition_function,
    const double temperature, /* [K] */
    const double temperature_ref, /* [K] */
    const double wavenumber) /* [cm^-1] */
 {
   const double C2 = 1.4388; /* 2nd Planck constant [K.cm] */
 
   const double fol = /* TODO ask to Yaniss why this variable is named fol */
     (1-exp(-C2*wavenumber/temperature))
   / (1-exp(-C2*wavenumber/temperature_ref));
 
   const double tmp =
     exp(-C2*lower_state_energy/temperature)
   / exp(-C2*lower_state_energy/temperature_ref);
 
   return intensity_ref * partition_function * tmp * fol ;
 }
 
 static res_T
 line_profile_factor
   (const struct sln_tree* tree,
    const struct shtr_line* shtr_line,
    double* out_profile_factor)
 {
   /* Star-HITRAN data */
   struct shtr_molecule molecule = SHTR_MOLECULE_NULL;
   const struct shtr_isotope* isotope = NULL;
 
   /* Mixture parameters */
   const struct sln_molecule* mol_params = NULL;
 
   /* Miscellaneous */
   double iso_abundance;
   double density; /* In molec.cm^-3 */
   double intensity, intensity_ref; /* In cm^-1/(molec.cm^2) */
   double Q, Q_T, Q_Tref; /* Partition function */
   double nu_c; /* In cm^-1 */
   double profile_factor; /* In m^-1.cm^-1 */
   double gj; /* State independant degeneracy factor */
   double Ps; /* In atm */
   double T; /* Temperature */
   int molid; /* Molecule id */
   int isoid; /* Isotope id local to its molecule */
 
   res_T res = RES_OK;
   ASSERT(tree && shtr_line && out_profile_factor);
 
   /* Fetch the molecule data */
   mol_params = tree->args.molecules + shtr_line->molecule_id;
   SHTR(isotope_metadata_find_molecule
     (tree->args.metadata, shtr_line->molecule_id, &molecule));
   ASSERT(!SHTR_MOLECULE_IS_NULL(&molecule));
   ASSERT(molecule.nisotopes > (size_t)shtr_line->isotope_id_local);
   isotope = molecule.isotopes + shtr_line->isotope_id_local;
 
   nu_c = line_center(shtr_line, tree->args.pressure);
 
   /* Compute the intensity */
   Ps = tree->args.pressure * mol_params->concentration;
   density = (AVOGADRO_NUMBER * Ps);
   density = density / (PERFECT_GAZ_CONSTANT * tree->args.temperature);
   density = density * 1e-6; /* Convert in molec.cm^-3 */
 
   /* Compute the partition function. TODO precompute it for molid/isoid */
   Q_Tref = isotope->Q296K;
   molid = shtr_line->molecule_id;
   isoid = shtr_line->isotope_id_local+1/*Local indices start at 1 in BD_TIPS*/;
   T = tree->args.temperature;
   BD_TIPS_2017(&molid, &T, &isoid, &gj, &Q_T);
   if(Q_T <= 0) {
     ERROR(tree->sln,
       "molecule %d: isotope %d: invalid partition function at T=%g\n",
       molid, isoid, T);
     res = RES_BAD_ARG;
     goto error;
   }
 
   Q = Q_Tref/Q_T;
 
   /* Compute the intensity */
   if(!mol_params->non_default_isotope_abundances) { /* Use default abundance */
     intensity_ref = shtr_line->intensity;
   } else {
     iso_abundance = mol_params->isotope_abundances[shtr_line->isotope_id_local];
     intensity_ref = shtr_line->intensity/isotope->abundance*iso_abundance;
   }
   intensity = line_intensity(intensity_ref, shtr_line->lower_state_energy, Q,
     tree->args.temperature, T_REF, nu_c);
 
   profile_factor = 1.e2 * density * intensity; /* In m^-1.cm^-1 */
 
 exit:
   *out_profile_factor = profile_factor;
   return res;
 error:
   profile_factor = NaN;
   goto exit;
 }
 
 /* Regularly mesh the interval [wavenumber, wavenumber+spectral_length[. Note
  * that the upper bound is excluded, this means that the last vertex of the
  * interval is not emitted */
 static INLINE res_T
 regular_mesh
   (const double wavenumber, /* Wavenumber where the mesh begins [cm^-1] */
    const double spectral_length, /* Size of the spectral interval to mesh [cm^-1] */
    const size_t nvertices, /* #vertices to issue */
    struct darray_double* wavenumbers) /* List of issued vertices */
 {
   /* Do not issue the vertex on the upper bound of the spectral range. That's
    * why we assume that the number of steps is equal to the number of vertices
    * and not to the number of vertices minus 1 */
   const double step = spectral_length / (double)nvertices;
   size_t ivtx;
   res_T res = RES_OK;
   ASSERT(wavenumber >= 0  && spectral_length > 0 && wavenumbers);
 
   FOR_EACH(ivtx, 0, nvertices) {
     const double nu = wavenumber + (double)ivtx*step;
     res = darray_double_push_back(wavenumbers, &nu);
     if(res != RES_OK) goto error;
   }
 exit:
   return res;
 error:
   goto exit;
 }
 
 /* The line is regularly discretized into a set of fragments of fixed size.
  * Their discretization is finer for the fragments around the center of the
  * line and becomes coarser as the fragments move away from it. Note that a
  * line is symmetrical in its center. As a consequence, the returned list is
  * only the set of wavenumbers from the line center to its upper bound. */
 static res_T
 regular_mesh_fragmented
   (const struct sln_tree* tree,
    const struct line* line,
    const size_t nvertices,
    struct darray_double* wavenumbers) /* List of issued vertices */
 {
   /* Fragment parameters */
   double fragment_length = 0;
   double fragment_nu_min = 0; /* Lower bound of the fragment */
   size_t fragment_nvtx = 0; /* #vertices into the fragment */
 
   /* Miscellaneous */
   const struct sln_molecule* mol_params = NULL;
   double line_nu_min = 0; /* In cm^-1 */
   double line_nu_max = 0; /* In cm^-1 */
   res_T res = RES_OK;
 
   ASSERT(tree && line && wavenumbers);
   ASSERT(IS_POW2(nvertices));
 
   /* TODO check mol params */
   mol_params = tree->args.molecules + line->molecule_id;
 
   /* Compute the spectral range of the line from its center to its cutoff */
   line_nu_min = line->wavenumber;
   line_nu_max = line->wavenumber + mol_params->cutoff;
 
   /* Define the size of a fragment as the width of the line at mid-height for a
    * Lorentz profile */
   fragment_length = line->gamma_l;
 
   /* Define the number of vertices for the first interval in [nu, gamma_l] */
   fragment_nu_min = line_nu_min;
   fragment_nvtx = MMAX(nvertices/2, 2);
 
   while(fragment_nu_min < line_nu_max) {
     const double spectral_length =
       MMIN(fragment_length, line_nu_max - fragment_nu_min);
 
     res = regular_mesh
       (fragment_nu_min, spectral_length, fragment_nvtx, wavenumbers);
     if(res != RES_OK) goto error;
 
     fragment_nu_min += fragment_length;
     fragment_nvtx = MMAX(fragment_nvtx/2, 2);
   }
 
   /* Register the last vertex, i.e. the upper bound of the spectral range */
   res = darray_double_push_back(wavenumbers, &line_nu_max);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   ERROR(tree->sln, "Error meshing the line -- %s.\n", res_to_cstr(res));
   goto exit;
 }
 
 /* Calculate line values for a set of wave numbers */
 static res_T
 eval_mesh
   (const struct sln_tree* tree,
    const struct line* line,
    const struct darray_double* wavenumbers,
    struct darray_double* values)
 {
   const double* nu = NULL;
   double* ha = NULL;
   size_t ivertex, nvertices;
   res_T res = RES_OK;
   ASSERT(tree && line && wavenumbers && values);
 
   nvertices = darray_double_size_get(wavenumbers);
   ASSERT(nvertices);
 
   res = darray_double_resize(values, nvertices);
   if(res != RES_OK) goto error;
 
   nu = darray_double_cdata_get(wavenumbers);
   ha = darray_double_data_get(values);
   FOR_EACH(ivertex, 0, nvertices) {
     ha[ivertex] = line_value(tree, line, nu[ivertex]);
   }
 
 exit:
   return res;
 error:
   ERROR(tree->sln, "Error evaluating the line mesh -- %s.\n", res_to_cstr(res));
   goto exit;
 }
 
 static void
 snap_mesh_to_upper_bound
   (const struct darray_double* wavenumbers,
    struct darray_double* values)
 {
   double* ha = NULL;
   size_t ivertex, nvertices;
   size_t ivertex_1st;
   ASSERT(wavenumbers && values);
   ASSERT(darray_double_size_get(wavenumbers) == darray_double_size_get(values));
   (void)wavenumbers;
 
   ha = darray_double_data_get(values);
   nvertices = darray_double_size_get(wavenumbers);
 
   /* Search for the first (non clipped) vertex */
   FOR_EACH(ivertex_1st, 0, nvertices) {
     if(ha[ivertex_1st] > 0) break;
   }
   ASSERT(ivertex_1st < nvertices);
 
   /* Ensure that the stored vertex value is an exclusive upper bound of the
    * original value. We do this by storing a value in single precision that is
    * strictly greater than its encoding in double precision */
   if(ha[ivertex_1st] != (float)ha[ivertex_1st]) {
     ha[ivertex_1st] = nextafterf((float)ha[ivertex_1st], FLT_MAX);
   }
 
   /* We have meshed the upper half of the line which is a strictly decreasing
    * function. To ensure that the mesh is an upper limit of this function,
    * simply align the value of each vertex with the value of the preceding
    * vertex */
   FOR_EACH_REVERSE(ivertex, nvertices-1, ivertex_1st) {
     if(ha[ivertex] < 0) continue; /* The vertex is clipped */
     ha[ivertex] = ha[ivertex-1];
   }
 }
 
 static INLINE int
 cmp_dbl(const void* a, const void* b)
 {
   const double key = *((const double*)a);
   const double item = *((const double*)b);
   if(key < item) return -1;
   if(key > item) return +1;
   return 0;
 }
 
 /* Return the value of the vertex whose wavenumber is greater than 'nu' */
 static INLINE double
 next_vertex_value
   (const double nu,
    const struct darray_double* wavenumbers,
    const struct darray_double* values)
 {
   const double* wnum = NULL;
   size_t ivertex = 0;
   ASSERT(wavenumbers && values);
 
   wnum = search_lower_bound
     (&nu,
      darray_double_cdata_get(wavenumbers),
      darray_double_size_get(wavenumbers),
      sizeof(double),
      cmp_dbl);
   ASSERT(wnum); /* It necessary exists */
 
   ivertex = (size_t)(wnum - darray_double_cdata_get(wavenumbers));
   ASSERT(ivertex < darray_double_size_get(values));
 
   return darray_double_cdata_get(values)[ivertex];
 }
 
 /* Append the line mesh into the vertices array */
 static res_T
 save_line_mesh
   (struct sln_tree* tree,
    const struct line* line,
    const struct darray_double* wavenumbers,
    const struct darray_double* values,
    size_t vertices_range[2]) /* Range into which the line vertices are saved */
 {
   const double* wnums = NULL;
   const double* vals = NULL;
   size_t nvertices = 0;
   size_t nwavenumbers = 0;
   size_t line_nvertices = 0;
   size_t ivertex = 0;
   size_t i = 0;
   res_T res = RES_OK;
 
   ASSERT(tree && line && wavenumbers && values && vertices_range);
   ASSERT(darray_double_size_get(wavenumbers) == darray_double_size_get(values));
 
   nvertices = darray_vertex_size_get(&tree->vertices);
   nwavenumbers = darray_double_size_get(wavenumbers);
 
   /* Compute the overall number of vertices of the line */
   line_nvertices = nwavenumbers
     * 2 /* The line is symmetrical in its center */
     - 1;/* Do not duplicate the line center */
 
   /* Allocate the list of line vertices */
   res = darray_vertex_resize(&tree->vertices, nvertices + line_nvertices);
   if(res != RES_OK) goto error;
 
   wnums = darray_double_cdata_get(wavenumbers);
   vals = darray_double_cdata_get(values);
 
   i = nvertices;
 
   #define MIRROR(Nu) (2*line->wavenumber - (Nu))
 
   /* Copy the vertices of the line for its lower half */
   FOR_EACH_REVERSE(ivertex, nwavenumbers-1, 0) {
     struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i++;
     const double nu = MIRROR(wnums[ivertex]);
     const double ha = vals[ivertex];
 
     vtx->wavenumber = (float)nu;
     vtx->ka = (float)ha;
   }
 
   /* Copy the vertices of the line for its upper half */
   FOR_EACH(ivertex, 0, nwavenumbers) {
     struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i++;
     const double nu = wnums[ivertex];
     const double ha = vals[ivertex];
 
     vtx->wavenumber = (float)nu;
     vtx->ka = (float)ha;
   }
 
   #undef MIRROR
 
   ASSERT(i == nvertices + line_nvertices);
 
   /* Setup the range of the line vertices */
   vertices_range[0] = nvertices;
   vertices_range[1] = i-1; /* Make the bound inclusive */
 
 exit:
   return res;
 error:
   darray_vertex_resize(&tree->vertices, nvertices);
   ERROR(tree->sln, "Error while recording line vertices -- %s.\n",
     res_to_cstr(res));
   goto exit;
 }
 
 /*******************************************************************************
  * Local function
  ******************************************************************************/
 res_T
 line_setup
   (const struct sln_tree* tree,
    const size_t iline,
    struct line* line)
 {
   struct shtr_molecule molecule = SHTR_MOLECULE_NULL;
   struct shtr_line shtr_line = SHTR_LINE_NULL;
   double molar_mass = 0; /* In kg.mol^-1 */
   const struct sln_molecule* mol_params = NULL;
   res_T res = RES_OK;
 
   ASSERT(tree && line);
 
   SHTR(line_list_at(tree->args.lines, iline, &shtr_line));
   SHTR(isotope_metadata_find_molecule
     (tree->args.metadata, shtr_line.molecule_id, &molecule));
   ASSERT(!SHTR_MOLECULE_IS_NULL(&molecule));
   ASSERT(molecule.nisotopes > (size_t)shtr_line.isotope_id_local);
 
   mol_params = tree->args.molecules + shtr_line.molecule_id;
 
   /* Convert the molar mass of the line from g.mol^-1 to kg.mol^-1 */
   molar_mass = molecule.isotopes[shtr_line.isotope_id_local].molar_mass*1e-3;
 
   /* Setup the line */
   res = line_profile_factor(tree, &shtr_line, &line->profile_factor);
   if(res != RES_OK) goto error;
 
   line->wavenumber = line_center(&shtr_line, tree->args.pressure);
   line->gamma_d = sln_compute_line_half_width_doppler
     (line->wavenumber, molar_mass, tree->args.temperature);
   line->gamma_l = sln_compute_line_half_width_lorentz(shtr_line.gamma_air,
     shtr_line.gamma_self, tree->args.pressure, mol_params->concentration);
   line->molecule_id = shtr_line.molecule_id;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 double
 line_value
   (const struct sln_tree* tree,
    const struct line* line,
    const double wavenumber)
 {
   const struct sln_molecule* mol_params = NULL;
   double profile = 0;
   ASSERT(tree && line);
 
   /* Retrieve the molecular parameters of the line to be mesh */
   mol_params = tree->args.molecules + line->molecule_id;
 
   if(wavenumber < line->wavenumber - mol_params->cutoff
   || wavenumber > line->wavenumber + mol_params->cutoff) {
     return 0;
   }
 
   switch(tree->args.line_profile) {
     case SLN_LINE_PROFILE_VOIGT:
       profile = sln_compute_voigt_profile
         (wavenumber, line->wavenumber, line->gamma_d, line->gamma_l);
       break;
     default: FATAL("Unreachable code.\n"); break;
   }
   return line->profile_factor * profile;
 }
 
 res_T
 line_mesh
   (struct sln_tree* tree,
    const size_t iline,
    const size_t nvertices_hint,
    size_t vertices_range[2]) /* out. Bounds are inclusive */
 {
   /* The line */
   struct line line = LINE_NULL;
 
   /* Temporary mesh */
   struct darray_double values; /* List of evaluated values */
   struct darray_double wavenumbers; /* List of considered wavenumbers */
   size_t nvertices_adjusted = 0; /* computed from nvertices_hint */
 
   /* Miscellaneous */
   res_T res = RES_OK;
 
   /* Pre-conditions */
   ASSERT(tree && nvertices_hint);
 
   darray_double_init(tree->sln->allocator, &values);
   darray_double_init(tree->sln->allocator, &wavenumbers);
 
   /* Setup the line wrt molecule concentration, isotope abundance, temperature
    * and pressure */
   res = line_setup(tree, iline, &line);
   if(res != RES_OK) goto error;
 
   /* Adjust the hint on the number of vertices. This is not actually the real
    * number of vertices but an adjusted hint on it. This new value ensures that
    * it is a power of 2 included in [MIN_NVERTICES_HINT, MAX_NVERTICES_HINT] */
   nvertices_adjusted = CLAMP
     (nvertices_hint, MIN_NVERTICES_HINT, MAX_NVERTICES_HINT);
   nvertices_adjusted = round_up_pow2(nvertices_adjusted);
 
   /* Emit the vertex coordinates, i.e. the wavenumbers */
   res = regular_mesh_fragmented(tree, &line, nvertices_adjusted, &wavenumbers);
   if(res != RES_OK) goto error;
 
   /* Evaluate the mesh vertices, i.e. define the line value for the list of
    * wavenumbers */
   eval_mesh(tree, &line, &wavenumbers, &values);
 
   switch(tree->args.mesh_type) {
     case SLN_MESH_UPPER:
       snap_mesh_to_upper_bound(&wavenumbers, &values);
       break;
     case SLN_MESH_FIT: /* Do nothing */ break;
     default: FATAL("Unreachable code.\n"); break;
   }
 
   res = save_line_mesh(tree, &line, &wavenumbers, &values, vertices_range);
   if(res != RES_OK) goto error;
 
 exit:
   darray_double_release(&values);
   darray_double_release(&wavenumbers);
   return res;
 error:
   goto exit;
 }