commit 350e81af171467d07413abd4d400c8a3304720ec
parent 181d4d63c4b31082ac5113f7e55554bb3d0e24ed
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 16 May 2022 16:45:40 +0200
Manage pressure and molecul cutoffs in lines view
Diffstat:
3 files changed, 165 insertions(+), 169 deletions(-)
diff --git a/src/shtr.h b/src/shtr.h
@@ -45,7 +45,7 @@
struct shtr_isotope {
double abundance; /* in ]0, 1] */
double Q296K; /* Partition function at Tref = 296K */
- double molar_mass; /* In g */
+ double molar_mass; /* In g.mol^-1 */
/* Local idx of the molecule to which the isotope belongs */
size_t molecule_id_local;
@@ -133,8 +133,9 @@ struct shtr_isotope_selection {
size_t nisotopes; /* 0 <=> select all the isotopes */
int32_t id; /* Molecule to which the isotopes belongs */
+ double cutoff; /* In cm^-1 */
};
-#define SHTR_ISOTOPE_SELECTION_NULL__ {0}
+#define SHTR_ISOTOPE_SELECTION_NULL__ {{0}, 0, 0, 0}
static const struct shtr_isotope_selection
SHTR_ISOTOPE_SELECTION_NULL = SHTR_ISOTOPE_SELECTION_NULL__;
@@ -144,8 +145,11 @@ struct shtr_lines_view_create_args {
/* List of molecule be selected */
struct shtr_isotope_selection molecules[SHTR_MAX_MOLECULES_COUNT];
size_t nmolecules;
+
+ double pressure; /* In atm. Used to compute the line center */
};
-#define SHTR_LINES_VIEW_CREATE_ARGS_NULL__ {0}
+#define SHTR_LINES_VIEW_CREATE_ARGS_NULL__ \
+ {{0,0}, {SHTR_ISOTOPE_SELECTION_NULL__}, 0, 0}
static const struct shtr_lines_view_create_args
SHTR_LINES_VIEW_CREATE_ARGS_NULL =
SHTR_LINES_VIEW_CREATE_ARGS_NULL__;
diff --git a/src/shtr_lines_view.c b/src/shtr_lines_view.c
@@ -23,7 +23,6 @@
#include "shtr_lines_list_c.h"
#include "shtr_log.h"
-#include <rsys/algorithm.h>
#include <rsys/cstr.h>
#include <rsys/dynamic_array_size_t.h>
@@ -35,200 +34,155 @@ struct shtr_lines_view {
ref_T ref;
};
+struct molecule_selection {
+ /* Map the isotope local identifier to a boolean defining if the isotope is
+ * selected or not */
+ char isotopes[SHTR_MAX_ISOTOPES_COUNT];
+ double cutoff; /* Molecule cutoff */
+};
+
/*******************************************************************************
* Helper functions
******************************************************************************/
static res_T
-check_shtr_lines_view_create_args
+check_shtr_isotope_selection
(struct shtr* shtr,
- const struct shtr_lines_view_create_args* args)
+ const char* caller,
+ const struct shtr_isotope_selection* molecule)
{
- if(!args) return RES_BAD_ARG;
- if(args->wavenumber_range[0] > args->wavenumber_range[1]) {
- log_err(shtr, "Invalid lines view spectral range [%g, %g].\n",
- args->wavenumber_range[0], args->wavenumber_range[1]);
- return RES_BAD_ARG;
- }
- return RES_OK;
-}
+ size_t i;
+ ASSERT(caller && molecule);
-static int
-cmp_line_wavenumber(const void* key, const void* entry)
-{
- const double wavenumber = *((double*)key);
- const struct shtr_line* line = entry;
- if(wavenumber < line->wavenumber) {
- return -1;
- } else if(wavenumber > line->wavenumber) {
- return 1;
- } else {
- return 0;
+ if((size_t)molecule->id >= SHTR_MAX_MOLECULES_COUNT) {
+ log_err(shtr,
+ "%s: molecule %d: invalid molecule identifier. "
+ "It must be less than %d.\n",
+ caller, molecule->id, SHTR_MAX_MOLECULES_COUNT);
+ return RES_BAD_ARG;
}
-}
-
-static res_T
-find_line_id_range
- (const struct shtr_lines_view* view,
- const double wavenumber_range[2],
- size_t line_id_range[2])
-{
- const struct shtr_line* lines = NULL;
- const struct shtr_line* l0 = NULL;
- const struct shtr_line* l1 = NULL;
- size_t nlines = 0;
- double exclusive_bound = 0;
- res_T res = RES_OK;
- ASSERT(view && wavenumber_range && line_id_range);
-
- res = shtr_lines_list_get(view->list, &lines);
- if(res != RES_OK) goto error;
- res = shtr_lines_list_get_size(view->list, &nlines);
- if(res != RES_OK) goto error;
- l0 = search_lower_bound(&wavenumber_range[0], lines, nlines,
- sizeof(struct shtr_line), cmp_line_wavenumber);
- if(l0) { /* No lower bound is found */
- line_id_range[0] = (size_t)(l0 - lines);
- } else {
- line_id_range[0] = SIZE_MAX;
- line_id_range[1] = 0;
- goto exit;
+ if(molecule->cutoff <= 0) {
+ log_err(shtr, "%s: molecule %d: invalid cutoff %g.\n",
+ caller, molecule->id, molecule->cutoff);
+ return RES_BAD_ARG;
}
- /* The search_lower_bound function returns the first entry that is greater
- * than _or_ equal to the search value. This is fine for finding an inclusive
- * lower bound, but not for finding an inclusive upper bound like below. In
- * the latter case, it is instead necessary to find the first entry greater
- * than the value sought and return the previous one. To do this, we
- * therefore seek the next representable floating point number directly above
- * the upper bound wavenumber.*/
- exclusive_bound = nextafter(wavenumber_range[1], INF);
- l1 = search_lower_bound(&exclusive_bound, lines, nlines,
- sizeof(struct shtr_line), cmp_line_wavenumber);
- if(!l1) { /* No upper bound is found. Take all the lines to the end */
- line_id_range[1] = nlines - 1;
- } else {
- line_id_range[1] = (size_t)(l1 - lines) - 1;
+ FOR_EACH(i, 0, molecule->nisotopes) {
+ if(molecule->isotope_ids_local[i] >= SHTR_MAX_ISOTOPES_COUNT) {
+ log_err(shtr,
+ "%s: molecule %d: isotope %d: invalid isotope local identifier. "
+ "It must be less than %d.\n",
+ caller,
+ molecule->id,
+ molecule->isotope_ids_local[i],
+ SHTR_MAX_ISOTOPES_COUNT);
+ return RES_BAD_ARG;
+ }
}
- ASSERT(line_id_range[0] <= line_id_range[1]);
- ASSERT(line_id_range[1] < nlines);
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static int
-cmp_isotope_selection(const void* a, const void* b)
-{
- const struct shtr_isotope_selection* molecule0 = a;
- const struct shtr_isotope_selection* molecule1 = b;
- ASSERT(a && b);
- return molecule0->id - molecule1->id;
-}
-
-static int
-cmp_isotope_selection_molecule_id(const void* key, const void* val)
-{
- const int32_t* molecule_id = key;
- const struct shtr_isotope_selection* molecule = val;
- ASSERT(key && val);
- return *molecule_id - molecule->id;
+ return RES_OK;
}
-static int
-cmp_int32(const void* a, const void* b)
-{
- ASSERT(a && b);
- return *((int32_t*)a) - *((int32_t*)b);
-}
-/* Return 1 if the line belongs to an isotope to be selected and 0 otherwise */
-static INLINE int
-is_the_line_to_select
- (const struct shtr_line* line,
- const struct shtr_isotope_selection* molecules,
- const size_t nmolecules)
+static res_T
+check_shtr_lines_view_create_args
+ (struct shtr* shtr,
+ const char* caller,
+ const struct shtr_lines_view_create_args* args)
{
- const struct shtr_isotope_selection* molecule = NULL;
- int32_t* isotope = NULL;
- ASSERT(molecules);
-
- molecule = search_lower_bound(&line->molecule_id, molecules, nmolecules,
- sizeof(*molecules), cmp_isotope_selection_molecule_id);
+ size_t i;
+ ASSERT(caller);
- if(molecule == NULL || molecule->id != line->molecule_id)
- return 0; /* The lines does not belongs to a molecule to be selected */
+ if(!args) return RES_BAD_ARG;
- if(molecule->nisotopes == 0) /* <=> All the isotopes are accepted */
- return 1;
+ if(args->wavenumber_range[0] > args->wavenumber_range[1]) {
+ log_err(shtr, "%s: invalid lines view spectral range [%g, %g].\n",
+ caller, args->wavenumber_range[0], args->wavenumber_range[1]);
+ return RES_BAD_ARG;
+ }
- isotope = search_lower_bound(&line->isotope_id_local,
- molecule->isotope_ids_local, molecule->nisotopes,
- sizeof(*molecule->isotope_ids_local), cmp_int32);
+ if(args->pressure < 0) {
+ log_err(shtr, "%s: invalid pressure %g.\n", caller, args->pressure);
+ return RES_BAD_ARG;
+ }
- if(isotope == NULL || *isotope != line->isotope_id_local)
- return 0; /* The lines does not belongs to an isotope to be selected */
+ FOR_EACH(i, 0, args->nmolecules) {
+ const res_T res = check_shtr_isotope_selection
+ (shtr, caller, &args->molecules[i]);
+ if(res != RES_OK) return res;
- return 1;
+ }
+ return RES_OK;
}
static res_T
select_lines
(struct shtr_lines_view* view,
+ const char* caller,
const struct shtr_lines_view_create_args* args)
{
const struct shtr_line* lines;
- struct shtr_isotope_selection molecules[SHTR_MAX_MOLECULES_COUNT];
- size_t line_id_range[2];
- size_t i;
+ struct molecule_selection selection[SHTR_MAX_MOLECULES_COUNT];
+ size_t imol;
+ size_t iiso;
+ size_t iline;
+ size_t nlines;
res_T res = RES_OK;
- ASSERT(view && args);
+ ASSERT(view && caller && args);
/* Nothing to do */
if(args->nmolecules == 0) goto exit;
- /* Sort in ascending order the ids of the molecules to include */
- memcpy(molecules, args->molecules, args->nmolecules*sizeof(molecules[0]));
- qsort
- (molecules,
- args->nmolecules,
- sizeof(molecules[0]),
- cmp_isotope_selection);
+ /* Setup the selection lookup table that map the isotope of a molecule to a
+ * boolean defining if it is selected */
+ memset(selection, 0, sizeof(selection));
+ FOR_EACH(imol, 0, args->nmolecules) {
+ const int32_t mol_id = args->molecules[imol].id;
+ ASSERT(mol_id < SHTR_MAX_MOLECULES_COUNT);
+ selection[mol_id].cutoff = args->molecules[imol].cutoff;
+
+ if(args->molecules[imol].nisotopes == 0) { /* Select all isotopes */
+ FOR_EACH(iiso, 0, SHTR_MAX_ISOTOPES_COUNT) {
+ selection[mol_id].isotopes[iiso] = 1;
+ }
- /* Sort in ascending order the ids of the par molecule isotopes to include */
- FOR_EACH(i, 0, args->nmolecules) {
- qsort
- (molecules[i].isotope_ids_local,
- molecules[i].nisotopes,
- sizeof(molecules[i].isotope_ids_local[0]),
- cmp_int32);
+ } else {
+ FOR_EACH(iiso, 0, args->molecules[imol].nisotopes) {
+ const int32_t iso_id = args->molecules[imol].isotope_ids_local[iiso];
+ ASSERT(iso_id < SHTR_MAX_ISOTOPES_COUNT);
+ selection[mol_id].isotopes[iso_id] = 1;
+ }
+ }
}
-#ifndef NDEBUG
- FOR_EACH(i, 1, args->nmolecules) {
- size_t j;
- ASSERT(molecules[i].id >= molecules[i-1].id);
- FOR_EACH(j, 1, molecules[i].nisotopes) {
- ASSERT(molecules[i].isotope_ids_local[j] >= molecules[i].isotope_ids_local[j-1]);
+ lines = darray_line_cdata_get(&view->list->lines);
+ nlines = darray_line_size_get(&view->list->lines);
+
+ /* Iterate through list of lines to find the ones to select based on spectral
+ * range and isotope selection */
+ FOR_EACH(iline, 0, nlines) {
+ const struct shtr_line* line = lines + iline;
+ double nu = 0;
+
+ /* The line is not selected */
+ if(selection[line->molecule_id].isotopes[line->isotope_id_local] == 0) {
+ continue;
}
- }
-#endif
- res = find_line_id_range(view, args->wavenumber_range, line_id_range);
- if(res != RES_OK) goto error;
+ /* Compute the line center for the submitted pressure */
+ nu = line->wavenumber + line->delta_air * args->pressure;
- lines = darray_line_cdata_get(&view->list->lines);
- FOR_EACH(i, line_id_range[0], line_id_range[1]+1) {
- if(is_the_line_to_select(&lines[i], molecules, args->nmolecules)) {
- res = darray_size_t_push_back(&view->line_ids, &i);
- if(res != RES_OK) {
- log_err(view->list->shtr,
- "Could not register the line into the view -- %s.\n",
- res_to_cstr(res));
- goto error;
- }
+ /* The line is out of the spectral range */
+ if(nu + selection[line->molecule_id].cutoff < args->wavenumber_range[0]
+ || nu - selection[line->molecule_id].cutoff > args->wavenumber_range[1]) {
+ continue;
+ }
+
+ res = darray_size_t_push_back(&view->line_ids, &iline);
+ if(res != RES_OK) {
+ log_err(view->list->shtr,
+ "%s: could not register the line into the view -- %s.\n",
+ caller, res_to_cstr(res));
+ goto error;
}
}
@@ -263,7 +217,7 @@ shtr_lines_view_create
res_T res = RES_OK;
if(!list || !out_view) { res = RES_BAD_ARG; goto error; }
- res = check_shtr_lines_view_create_args(list->shtr, args);
+ res = check_shtr_lines_view_create_args(list->shtr, FUNC_NAME, args);
if(res != RES_OK) goto error;
view = MEM_CALLOC(list->shtr->allocator, 1, sizeof(*view));
if(!view) {
@@ -276,7 +230,7 @@ shtr_lines_view_create
view->list = list;
darray_size_t_init(list->shtr->allocator, &view->line_ids);
- res = select_lines(view, args);
+ res = select_lines(view, FUNC_NAME, args);
if(res != RES_OK) goto error;
exit:
diff --git a/src/test_shtr_lines.c b/src/test_shtr_lines.c
@@ -251,12 +251,9 @@ check_view
CHK(shtr_lines_list_get(list, &lines) == RES_OK);
FOR_EACH(i, 0, n) {
const struct shtr_line* line = NULL;
+ double nu = 0;
size_t imolecule;
- /* Discard lines that are out of ranges */
- if(lines[i].wavenumber < args->wavenumber_range[0]) continue;
- if(lines[i].wavenumber > args->wavenumber_range[1]) break;
-
/* Discard lines that does not belongs to a listed molecule */
FOR_EACH(imolecule, 0, args->nmolecules) {
if(lines[i].molecule_id == args->molecules[imolecule].id) break;
@@ -264,6 +261,15 @@ check_view
if(imolecule >= args->nmolecules)
continue; /* Skip the molecule */
+ /* Compute the line center for the given pressure */
+ nu = lines[i].wavenumber + lines[i].delta_air * args->pressure;
+
+ /* Discard lines that are out of ranges */
+ if(nu + args->molecules[imolecule].cutoff < args->wavenumber_range[0])
+ continue;
+ if(nu - args->molecules[imolecule].cutoff > args->wavenumber_range[1])
+ continue;
+
/* nisotope == 0 <=> All isotopes must be selected */
if(args->molecules[imolecule].nisotopes) {
size_t iisotope;
@@ -355,6 +361,8 @@ test_view(struct shtr* shtr)
};
const size_t nlines = sizeof(l) / sizeof(struct shtr_line);
+ const double cutoff = 0.01;
+
struct shtr_lines_view_create_args args = SHTR_LINES_VIEW_CREATE_ARGS_NULL;
struct shtr_lines_list* list = NULL;
struct shtr_lines_view* view = NULL;
@@ -374,7 +382,11 @@ test_view(struct shtr* shtr)
args.molecules[0].id = 1;
args.molecules[1].id = 2;
args.molecules[2].id = 3;
+ args.molecules[0].cutoff = cutoff;
+ args.molecules[1].cutoff = cutoff;
+ args.molecules[2].cutoff = cutoff;
args.nmolecules = 3;
+ args.pressure = 0;
CHK(shtr_lines_view_create(NULL, &args, &view) == RES_BAD_ARG);
CHK(shtr_lines_view_create(list, NULL, &view) == RES_BAD_ARG);
CHK(shtr_lines_view_create(list, &args, NULL) == RES_BAD_ARG);
@@ -412,9 +424,16 @@ test_view(struct shtr* shtr)
CHK(n == 0);
CHK(shtr_lines_view_ref_put(view) == RES_OK);
- args.wavenumber_range[0] = 2;
- args.wavenumber_range[1] = INF;
+ args.pressure = -1;
+ CHK(shtr_lines_view_create(list, &args, &view) == RES_BAD_ARG);
+ args.pressure = 0;
+ args.molecules[0].cutoff = -1;
args.nmolecules = 3;
+ CHK(shtr_lines_view_create(list, &args, &view) == RES_BAD_ARG);
+ args.molecules[0].cutoff = cutoff;
+
+ args.wavenumber_range[0] = nextafter(l[nlines-1].wavenumber, INF) + cutoff;
+ args.wavenumber_range[1] = INF;
CHK(shtr_lines_view_create(list, &args, &view) == RES_OK);
CHK(shtr_lines_view_get_size(view, &n) == RES_OK);
CHK(n == 0);
@@ -436,7 +455,7 @@ test_view(struct shtr* shtr)
args.wavenumber_range[1] = 0.29477;
CHK(shtr_lines_view_create(list, &args, &view) == RES_OK);
CHK(shtr_lines_view_get_size(view, &n) == RES_OK);
- CHK(n == 1);
+ CHK(n == 6);
check_view(list, view, &args);
CHK(shtr_lines_view_ref_put(view) == RES_OK);
@@ -444,7 +463,7 @@ test_view(struct shtr* shtr)
args.wavenumber_range[1] = 0.21283;
CHK(shtr_lines_view_create(list, &args, &view) == RES_OK);
CHK(shtr_lines_view_get_size(view, &n) == RES_OK);
- CHK(n == 10);
+ CHK(n == 14);
check_view(list, view, &args);
CHK(shtr_lines_view_ref_put(view) == RES_OK);
@@ -467,7 +486,7 @@ test_view(struct shtr* shtr)
args.nmolecules = 2;
CHK(shtr_lines_view_create(list, &args, &view) == RES_OK);
CHK(shtr_lines_view_get_size(view, &n) == RES_OK);
- CHK(n == 6);
+ CHK(n == 10);
check_view(list, view, &args);
CHK(shtr_lines_view_ref_put(view) == RES_OK);
@@ -479,15 +498,34 @@ test_view(struct shtr* shtr)
check_view(list, view, &args);
CHK(shtr_lines_view_ref_put(view) == RES_OK);
+ args.molecules[0].id = 1;
+ args.molecules[0].cutoff = 1e-6;
+ args.molecules[0].nisotopes = 0;
+ args.nmolecules = 1;
+ args.wavenumber_range[0] = 0.1899;
+ args.wavenumber_range[1] = 0.195;
+ CHK(shtr_lines_view_create(list, &args, &view) == RES_OK);
+ CHK(shtr_lines_view_get_size(view, &n) == RES_OK);
+ CHK(n == 1);
+ CHK(shtr_lines_view_ref_put(view) == RES_OK);
+ args.pressure = 1;
+ CHK(shtr_lines_view_create(list, &args, &view) == RES_OK);
+ CHK(shtr_lines_view_get_size(view, &n) == RES_OK);
+ CHK(n == 2);
+ check_view(list, view, &args);
+ CHK(shtr_lines_view_ref_put(view) == RES_OK);
+
args.wavenumber_range[0] = 0;
args.wavenumber_range[1] = INF;
args.molecules[0].id = 3;
+ args.molecules[0].cutoff = cutoff;
args.molecules[0].isotope_ids_local[0] = 4;
args.molecules[0].isotope_ids_local[1] = 1;
args.molecules[0].isotope_ids_local[2] = 2;
args.molecules[0].isotope_ids_local[3] = 0;
args.molecules[0].nisotopes = 4;
args.molecules[1].id = 1;
+ args.molecules[0].cutoff = cutoff;
args.molecules[1].isotope_ids_local[0] = 4;
args.molecules[1].isotope_ids_local[1] = 3;
args.molecules[1].nisotopes = 2;
