star-hitran

Load line-by-line data from the HITRAN database
git clone git://git.meso-star.fr/star-hitran.git
Log | Files | Refs | README | LICENSE

commit c67215c15a7603efc011203a477c8dff0d9b5f46
parent 9f23dda14b1405092614e70a9778e723a14ed28d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 11 Mar 2022 09:31:52 +0100

Add isotope selection per molecule to the lines view

Diffstat:
Msrc/shtr.h | 58+++++++++++++++++++++++++++++++++++++---------------------
Msrc/shtr_lines_view.c | 79+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------
Msrc/test_shtr_lines.c | 51++++++++++++++++++++++++++++++++++++++++++---------
3 files changed, 144 insertions(+), 44 deletions(-)

diff --git a/src/shtr.h b/src/shtr.h @@ -40,27 +40,7 @@ #endif #define SHTR_MAX_MOLECULES_COUNT 100 - -struct shtr_create_args { - struct logger* logger; /* May be NULL <=> default logger */ - struct mem_allocator* allocator; /* NULL <=> use default allocator */ - int verbose; /* Verbosity level */ -}; -#define SHTR_CREATE_ARGS_DEFAULT__ {NULL, NULL, 0} -static const struct shtr_create_args SHTR_CREATE_ARGS_DEFAULT = - SHTR_CREATE_ARGS_DEFAULT__; - -struct shtr_lines_view_create_args { - double wavenumber_range[2]; /* Spectral range */ - - /* List of molecule identifiers to consider sorted in ascending order */ - int32_t molecules[SHTR_MAX_MOLECULES_COUNT]; - size_t nmolecules; -}; -#define SHTR_LINES_VIEW_CREATE_ARGS_NULL__ {{0,0},{0}, 0} -static const struct shtr_lines_view_create_args -SHTR_LINES_VIEW_CREATE_ARGS_NULL = - SHTR_LINES_VIEW_CREATE_ARGS_NULL__; +#define SHTR_MAX_ISOTOPES_COUNT 10 struct shtr_isotope { double abundance; /* in ]0, 1] */ @@ -119,6 +99,42 @@ struct shtr_isotope_metadata; struct shtr_lines_list; struct shtr_lines_view; +/******************************************************************************* + * Input arguments for API functions + ******************************************************************************/ +struct shtr_create_args { + struct logger* logger; /* May be NULL <=> default logger */ + struct mem_allocator* allocator; /* NULL <=> use default allocator */ + int verbose; /* Verbosity level */ +}; +#define SHTR_CREATE_ARGS_DEFAULT__ {0} +static const struct shtr_create_args SHTR_CREATE_ARGS_DEFAULT = + SHTR_CREATE_ARGS_DEFAULT__; + +struct shtr_isotope_selection { + /* List of isotopes identifier to consider for this molecule. The listed + * idenfiers are _local_ to the molecule */ + int32_t isotope_ids_local[SHTR_MAX_ISOTOPES_COUNT]; + size_t nisotopes; /* 0 <=> select all the isotopes */ + + int32_t id; /* Molecule to which the isotopes belongs */ +}; +#define SHTR_ISOTOPE_SELECTION_NULL__ {0} +static const struct shtr_isotope_selection +SHTR_ISOTOPE_SELECTION_NULL = SHTR_ISOTOPE_SELECTION_NULL__; + +struct shtr_lines_view_create_args { + double wavenumber_range[2]; /* Spectral range */ + + /* List of molecule be selected */ + struct shtr_isotope_selection molecules[SHTR_MAX_MOLECULES_COUNT]; + size_t nmolecules; +}; +#define SHTR_LINES_VIEW_CREATE_ARGS_NULL__ {0} +static const struct shtr_lines_view_create_args +SHTR_LINES_VIEW_CREATE_ARGS_NULL = + SHTR_LINES_VIEW_CREATE_ARGS_NULL__; + BEGIN_DECLS /******************************************************************************* diff --git a/src/shtr_lines_view.c b/src/shtr_lines_view.c @@ -120,25 +120,58 @@ error: } 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; +} + +static int cmp_int32(const void* a, const void* b) { ASSERT(a && b); return *((int32_t*)a) - *((int32_t*)b); } -/* Return 1 if `molecule' is in `molecules' and 0 otherwise */ +/* Return 1 if the line belongs to an isotope to be selected and 0 otherwise */ static INLINE int -include_molecule - (const int32_t molecules[], - const size_t nmolecules, - const int32_t molecule) +is_the_line_to_select + (const struct shtr_line* line, + const struct shtr_isotope_selection* molecules, + const size_t nmolecules) { - int32_t* i = NULL; + const struct shtr_isotope_selection* molecule = NULL; + int32_t* isotope = NULL; ASSERT(molecules); - i = search_lower_bound - (&molecule, molecules, nmolecules, sizeof(int32_t), cmp_int32); - return i != NULL && *i == molecule; + molecule = search_lower_bound(&line->molecule_id, molecules, nmolecules, + sizeof(*molecules), cmp_isotope_selection_molecule_id); + + if(molecule == NULL || molecule->id != line->molecule_id) + return 0; /* The lines does not belongs to a molecule to be selected */ + + if(molecule->nisotopes == 0) /* <=> All the isotopes are accepted */ + return 1; + + isotope = search_lower_bound(&line->isotope_id_local, + molecule->isotope_ids_local, molecule->nisotopes, + sizeof(*molecule->isotope_ids_local), cmp_int32); + + if(isotope == NULL || *isotope != line->isotope_id_local) + return 0; /* The lines does not belongs to an isotope to be selected */ + + return 1; } static res_T @@ -147,7 +180,7 @@ select_lines const struct shtr_lines_view_create_args* args) { const struct shtr_line* lines; - int32_t molecules[SHTR_MAX_MOLECULES_COUNT]; + struct shtr_isotope_selection molecules[SHTR_MAX_MOLECULES_COUNT]; size_t line_id_range[2]; size_t i; res_T res = RES_OK; @@ -156,12 +189,30 @@ select_lines /* Nothing to do */ if(args->nmolecules == 0) goto exit; - /* Sort in ascending order the indices of the molecules to include */ + /* 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_int32); + qsort + (molecules, + args->nmolecules, + sizeof(molecules[0]), + cmp_isotope_selection); + + /* 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); + } + #ifndef NDEBUG FOR_EACH(i, 1, args->nmolecules) { - ASSERT(molecules[i] >= molecules[i-1]); + 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]); + } } #endif @@ -170,7 +221,7 @@ select_lines lines = darray_line_cdata_get(&view->list->lines); FOR_EACH(i, line_id_range[0], line_id_range[1]+1) { - if(include_molecule(molecules, args->nmolecules, lines[i].molecule_id)) { + 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, diff --git a/src/test_shtr_lines.c b/src/test_shtr_lines.c @@ -232,9 +232,22 @@ check_view /* Discard lines that does not belongs to a listed molecule */ FOR_EACH(imolecule, 0, args->nmolecules) { - if(lines[i].molecule_id == args->molecules[imolecule]) break; + if(lines[i].molecule_id == args->molecules[imolecule].id) break; + } + if(imolecule >= args->nmolecules) + continue; /* Skip the molecule */ + + /* nisotope == 0 <=> All isotopes must be selected */ + if(args->molecules[imolecule].nisotopes) { + size_t iisotope; + FOR_EACH(iisotope, 0, args->molecules[imolecule].nisotopes) { + if(lines[i].isotope_id_local + == args->molecules[imolecule].isotope_ids_local[iisotope]) + break; + } + if(iisotope >= args->molecules[imolecule].nisotopes) + continue; /* Skip the isotope */ } - if(imolecule >= args->nmolecules) continue; CHK(shtr_lines_view_get_line(view, nlines_selected, &line) == RES_OK); CHK(line_eq(line, lines+i)); @@ -331,9 +344,9 @@ test_view(struct shtr* shtr) args.wavenumber_range[0] = 0; args.wavenumber_range[1] = INF; - args.molecules[0] = 1; - args.molecules[1] = 2; - args.molecules[2] = 3; + args.molecules[0].id = 1; + args.molecules[1].id = 2; + args.molecules[2].id = 3; args.nmolecules = 3; CHK(shtr_lines_view_create(NULL, &args, &view) == RES_BAD_ARG); CHK(shtr_lines_view_create(list, NULL, &view) == RES_BAD_ARG); @@ -408,14 +421,14 @@ test_view(struct shtr* shtr) check_view(list, view, &args); CHK(shtr_lines_view_ref_put(view) == RES_OK); - args.molecules[0] = 2; + args.molecules[0].id = 2; args.nmolecules = 1; CHK(shtr_lines_view_create(list, &args, &view) == RES_OK); CHK(shtr_lines_view_get_size(view, &n) == RES_OK); CHK(n == 0); CHK(shtr_lines_view_ref_put(view) == RES_OK); - args.molecules[1] = 1; + args.molecules[1].id = 1; args.nmolecules = 2; CHK(shtr_lines_view_create(list, &args, &view) == RES_OK); CHK(shtr_lines_view_get_size(view, &n) == RES_OK); @@ -423,7 +436,7 @@ test_view(struct shtr* shtr) check_view(list, view, &args); CHK(shtr_lines_view_ref_put(view) == RES_OK); - args.molecules[1] = 3; + args.molecules[1].id = 3; args.nmolecules = 2; CHK(shtr_lines_view_create(list, &args, &view) == RES_OK); CHK(shtr_lines_view_get_size(view, &n) == RES_OK); @@ -433,12 +446,32 @@ test_view(struct shtr* shtr) args.wavenumber_range[0] = 0; args.wavenumber_range[1] = INF; - args.molecules[0] = 2; + args.molecules[0].id = 2; args.nmolecules = 1; CHK(shtr_lines_view_create(list, &args, &view) == RES_OK); 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; /* 24 */ + 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; /* 11 */ + args.molecules[1].isotope_ids_local[0] = 4; + args.molecules[1].isotope_ids_local[1] = 3; + args.molecules[1].nisotopes = 2; + 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 == 35); + check_view(list, view, &args); + CHK(shtr_lines_view_ref_put(view) == RES_OK); + + CHK(shtr_lines_list_ref_put(list) == RES_OK); CHK(fclose(fp) == 0); }