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

test_shtr_lines.c (13449B)


      1 /* Copyright (C) 2022, 2025, 2026 |Méso|Star> (contact@meso-star.com)
      2  * Copyright (C) 2025, 2026 Université de Lorraine
      3  * Copyright (C) 2022 Centre National de la Recherche Scientifique
      4  * Copyright (C) 2022 Université Paul Sabatier
      5  *
      6  * This program is free software: you can redistribute it and/or modify
      7  * it under the terms of the GNU General Public License as published by
      8  * the Free Software Foundation, either version 3 of the License, or
      9  * (at your option) any later version.
     10  *
     11  * This program is distributed in the hope that it will be useful,
     12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     14  * GNU General Public License for more details.
     15  *
     16  * You should have received a copy of the GNU General Public License
     17  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     18 
     19 #define _POSIX_C_SOURCE 200112L
     20 
     21 #include "shtr.h"
     22 
     23 #include <rsys/clock_time.h>
     24 #include <rsys/mem_allocator.h>
     25 #include <rsys/math.h>
     26 
     27 #include <math.h>
     28 #include <string.h>
     29 
     30 static void
     31 print_lines
     32   (FILE* fp,
     33    const struct shtr_line* lines,
     34    const size_t nlines)
     35 {
     36   size_t i;
     37 
     38   CHK(fp && (!nlines || lines));
     39   FOR_EACH(i, 0, nlines) {
     40     fprintf(fp,
     41       "%2d%1d%12.6f%10.3e 0.000E-00.%04d%5.3f%10.4f%4.2f%8.6f"
     42       /* Dummy remaining data */
     43       "          0 0 0" /* Global upper quanta */
     44       "          0 0 0" /* Global upper quanta */
     45       "  5  5  0      " /* Local upper quanta */
     46       "  5  5  1      " /* Local lower quanta */
     47       "562220" /* Error indices */
     48       "5041 7833348" /* References */
     49       " " /* Line mixing flag */
     50       "   66.0" /* g' */
     51       "   66.0" /* g'' */
     52       "\n",
     53       lines[i].molecule_id,
     54       lines[i].isotope_id_local == 9 ? 0 : lines[i].isotope_id_local+1,
     55       lines[i].wavenumber,
     56       lines[i].intensity,
     57       (int)(lines[i].gamma_air*10000+0.5/*round*/),
     58       lines[i].gamma_self,
     59       lines[i].lower_state_energy,
     60       lines[i].n_air,
     61       lines[i].delta_air);
     62   }
     63 }
     64 
     65 static void
     66 test_line_eq(void)
     67 {
     68   const struct shtr_line l0 = {
     69     0.000134, 2.672E-38, 0.0533, 0.410, 608.4727, 0.79, 0.000060, 1, 4
     70   };
     71   struct shtr_line l1 = l0;
     72 
     73   CHK(shtr_line_eq(&l0, &l1));
     74 
     75   l1.wavenumber = nextafter(l1.wavenumber, INF);
     76   CHK(!shtr_line_eq(&l0, &l1));
     77   l1.wavenumber = l0.wavenumber;
     78   l1.intensity = nextafter(l1.intensity, INF);
     79   CHK(!shtr_line_eq(&l0, &l1));
     80   l1.intensity = l0.intensity;
     81   l1.gamma_air = nextafter(l1.gamma_air, INF);
     82   CHK(!shtr_line_eq(&l0, &l1));
     83   l1.gamma_air = l0.gamma_air;
     84   l1.gamma_self = nextafter(l1.gamma_self, INF);
     85   CHK(!shtr_line_eq(&l0, &l1));
     86   l1.gamma_self = l0.gamma_self;
     87   l1.lower_state_energy = nextafter(l1.lower_state_energy, INF);
     88   CHK(!shtr_line_eq(&l0, &l1));
     89   l1.lower_state_energy = l0.lower_state_energy;
     90   l1.n_air = nextafter(l1.n_air, INF);
     91   CHK(!shtr_line_eq(&l0, &l1));
     92   l1.n_air = l0.n_air;
     93   l1.delta_air = nextafter(l1.delta_air, INF);
     94   CHK(!shtr_line_eq(&l0, &l1));
     95   l1.delta_air = l0.delta_air;
     96   l1.molecule_id += 1;
     97   CHK(!shtr_line_eq(&l0, &l1));
     98   l1.molecule_id = l0.molecule_id;
     99   l1.isotope_id_local += 1;
    100   CHK(!shtr_line_eq(&l0, &l1));
    101   l1.isotope_id_local = l0.isotope_id_local;
    102 
    103   CHK(shtr_line_eq(&l0, &l1));
    104 }
    105 
    106 static void
    107 line_eq_eps
    108   (const struct shtr_line* ln0,
    109    const struct shtr_line* ln1,
    110    struct shtr_line_list_info* info)
    111 {
    112   #define EQ_EPS(Name) CHK(eq_eps(ln0->Name, ln1->Name, info->Name.err))
    113   EQ_EPS(wavenumber);
    114   EQ_EPS(intensity);
    115   EQ_EPS(gamma_air);
    116   EQ_EPS(gamma_self);
    117   EQ_EPS(lower_state_energy);
    118   EQ_EPS(n_air);
    119   EQ_EPS(delta_air);
    120   #undef EQ_EPS
    121   CHK(ln0->molecule_id == ln1->molecule_id);
    122   CHK(ln0->isotope_id_local == ln1->isotope_id_local);
    123 }
    124 
    125 static void
    126 test_load(struct shtr* shtr)
    127 {
    128   const struct shtr_line l[] = {
    129     {0.000134, 2.672E-38, 0.0533, 0.410, 608.4727, 0.79, 0.000060, 1, 4},
    130     {0.000379, 1.055E-39, 0.0418, 0.329,1747.9686, 0.79, 0.000110, 1, 5},
    131     {0.000448, 5.560E-38, 0.0490, 0.364,1093.0269, 0.79, 0.000060, 1, 4},
    132     {0.000686, 1.633E-36, 0.0578, 0.394, 701.1162, 0.79, 0.000180, 1, 4},
    133     {0.000726, 6.613E-33, 0.0695, 0.428, 402.3295, 0.79, 0.000240, 1, 3}
    134   };
    135   const size_t nlines = sizeof(l) / sizeof(struct shtr_line);
    136 
    137   struct shtr_line_list* list = NULL;
    138   struct shtr_line_list_info info = SHTR_LINE_LIST_INFO_NULL;
    139   struct shtr_line_list_info info_ref = SHTR_LINE_LIST_INFO_NULL;
    140   struct shtr_line line = SHTR_LINE_NULL;
    141   struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL__;
    142   const char* filename = "test_lines.txt";
    143   FILE* fp = NULL;
    144   size_t i, n;
    145 
    146   CHK(fp = fopen(filename, "w+"));
    147   print_lines(fp, l, nlines);
    148   rewind(fp);
    149 
    150   args.file = fp;
    151   CHK(shtr_line_list_load(NULL, &args, &list) == RES_BAD_ARG);
    152   CHK(shtr_line_list_load(shtr, NULL, &list) == RES_BAD_ARG);
    153   CHK(shtr_line_list_load(shtr, &args, NULL) == RES_BAD_ARG);
    154   CHK(shtr_line_list_load(shtr, &args, &list) == RES_OK);
    155 
    156   CHK(shtr_line_list_get_info(NULL, &info) == RES_BAD_ARG);
    157   CHK(shtr_line_list_get_info(list, NULL) == RES_BAD_ARG);
    158   CHK(shtr_line_list_get_info(list, &info) == RES_OK);
    159 
    160   /* Setup the reference information, i.e., from the original list of lines */
    161   FOR_EACH(i, 0, nlines) {
    162     #define UPDATE_INFO(Name) { \
    163       info_ref.Name.range[0] = MMIN(l[i].Name, info_ref.Name.range[0]); \
    164       info_ref.Name.range[1] = MMAX(l[i].Name, info_ref.Name.range[1]); \
    165     } (void)0
    166     UPDATE_INFO(wavenumber);
    167     UPDATE_INFO(intensity);
    168     UPDATE_INFO(gamma_air);
    169     UPDATE_INFO(gamma_self);
    170     UPDATE_INFO(lower_state_energy);
    171     UPDATE_INFO(n_air);
    172     UPDATE_INFO(delta_air);
    173     #undef UPDATE_INFO
    174   }
    175 
    176   #define CHK_INFO(Name) { \
    177     CHK(eq_eps(info.Name.range[0], info_ref.Name.range[0], info.Name.err)); \
    178     CHK(eq_eps(info.Name.range[1], info_ref.Name.range[1], info.Name.err)); \
    179   } (void)0
    180   CHK_INFO(wavenumber);
    181   CHK_INFO(intensity);
    182   CHK_INFO(gamma_air);
    183   CHK_INFO(gamma_self);
    184   CHK_INFO(lower_state_energy);
    185   CHK_INFO(n_air);
    186   CHK_INFO(delta_air);
    187   #undef CHK_INFO
    188 
    189   CHK(shtr_line_list_get_size(NULL, &n) == RES_BAD_ARG);
    190   CHK(shtr_line_list_get_size(list, NULL) == RES_BAD_ARG);
    191   CHK(shtr_line_list_get_size(list, &n) == RES_OK);
    192   CHK(n == nlines);
    193 
    194   CHK(shtr_line_list_at(NULL, 0, &line) == RES_BAD_ARG);
    195   CHK(shtr_line_list_at(list, nlines, &line) == RES_BAD_ARG);
    196   CHK(shtr_line_list_at(list, 0, NULL) == RES_BAD_ARG);
    197   FOR_EACH(i, 0, n) {
    198     CHK(shtr_line_list_at(list, i, &line) == RES_OK);
    199     line_eq_eps(&line, l+i, &info);
    200   }
    201 
    202   CHK(shtr_line_list_ref_get(NULL) == RES_BAD_ARG);
    203   CHK(shtr_line_list_ref_get(list) == RES_OK);
    204   CHK(shtr_line_list_ref_put(NULL) == RES_BAD_ARG);
    205   CHK(shtr_line_list_ref_put(list) == RES_OK);
    206   CHK(shtr_line_list_ref_put(list) == RES_OK);
    207 
    208   CHK(fclose(fp) == 0);
    209 
    210   args.file = NULL;
    211   args.filename = filename;
    212   CHK(shtr_line_list_load(NULL, &args, &list) == RES_BAD_ARG);
    213   CHK(shtr_line_list_load(shtr, NULL, &list) == RES_BAD_ARG);
    214   CHK(shtr_line_list_load(shtr, &args, NULL) == RES_BAD_ARG);
    215   CHK(shtr_line_list_load(shtr, &args, &list) == RES_OK);
    216 
    217   CHK(shtr_line_list_get_info(list, &info) == RES_OK);
    218 
    219   CHK(shtr_line_list_get_size(list, &n) == RES_OK);
    220   CHK(n == nlines);
    221 
    222   FOR_EACH(i, 0, n) {
    223     CHK(shtr_line_list_at(list, i, &line) == RES_OK);
    224     line_eq_eps(&line, l+i, &info);
    225   }
    226 
    227   CHK(shtr_line_list_ref_put(list) == RES_OK);
    228 }
    229 
    230 static void
    231 test_line
    232   (struct shtr* shtr, const struct shtr_line* ln, const res_T res)
    233 {
    234   struct shtr_line_list* list = NULL;
    235   struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL__;
    236 
    237   CHK(args.file = tmpfile());
    238   print_lines(args.file, ln, 1);
    239   rewind(args.file);
    240   CHK(shtr_line_list_load(shtr, &args, &list) == res);
    241   CHK(fclose(args.file) == 0);
    242 
    243   if(res == RES_OK) CHK(shtr_line_list_ref_put(list) == RES_OK);
    244 }
    245 
    246 static void
    247 test_load_failures(struct shtr* shtr)
    248 {
    249   const struct shtr_line ln_ref = {
    250     0.000134, 2.672E-38, 0.0533, 0.410, 608.4727, 0.79, 0.000060, 1, 4,
    251   };
    252   struct shtr_line ln;
    253   struct shtr_line_list* list = NULL;
    254   struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL__;
    255 
    256   /* Check that the reference line is valid */
    257   test_line(shtr, &ln_ref, RES_OK);
    258 
    259   /* Invalid wavenumber */
    260   ln = ln_ref; ln.wavenumber = 0;
    261   test_line(shtr, &ln, RES_BAD_ARG);
    262 
    263   /* Invalid intensity */
    264   ln = ln_ref; ln.intensity = 0;
    265   test_line(shtr, &ln, RES_BAD_ARG);
    266 
    267   /* Invalid gamma air */
    268   ln = ln_ref; ln.gamma_air = -1;
    269   test_line(shtr, &ln, RES_BAD_ARG);
    270 
    271   /* Invalid gamma self */
    272   ln = ln_ref; ln.gamma_self = -1;
    273   test_line(shtr, &ln, RES_BAD_ARG);
    274 
    275   /* Unavailable lower state energy */
    276   ln = ln_ref; ln.lower_state_energy = -1;
    277   test_line(shtr, &ln, RES_OK);
    278 
    279   /* Invalid lower state energy */
    280   ln = ln_ref; ln.lower_state_energy = -2;
    281   test_line(shtr, &ln, RES_BAD_ARG);
    282 
    283   /* Invalid molecule id */
    284   ln = ln_ref; ln.molecule_id = -1;
    285   test_line(shtr, &ln, RES_BAD_ARG);
    286 
    287   /* Bad file formatting */
    288   ln = ln_ref; ln.molecule_id = 100;
    289   test_line(shtr, &ln, RES_BAD_ARG);
    290   ln = ln_ref; ln.isotope_id_local = 10;
    291   test_line(shtr, &ln, RES_BAD_ARG);
    292 
    293   /* Lines are not correctly sorted */
    294   CHK(args.file = tmpfile());
    295   ln = ln_ref;
    296   print_lines(args.file, &ln, 1);
    297   ln.wavenumber -= 1e-4;
    298   print_lines(args.file, &ln, 1);
    299   rewind(args.file);
    300   CHK(shtr_line_list_load(shtr, &args, &list) == RES_BAD_ARG);
    301   CHK(fclose(args.file) == 0);
    302 }
    303 
    304 static void
    305 check_line_list_equality
    306   (struct shtr_line_list* list1,
    307    struct shtr_line_list* list2)
    308 {
    309   size_t n1, n2;
    310   size_t iline, nlines;
    311   CHK(list1 && list2);
    312 
    313   CHK(shtr_line_list_get_size(list1, &n1) == RES_OK);
    314   CHK(shtr_line_list_get_size(list2, &n2) == RES_OK);
    315   CHK(n1 == n2);
    316   nlines = n1;
    317 
    318   FOR_EACH(iline, 0, nlines) {
    319     struct shtr_line lines1 = SHTR_LINE_NULL;
    320     struct shtr_line lines2 = SHTR_LINE_NULL;
    321 
    322     CHK(shtr_line_list_at(list1, iline, &lines1) == RES_OK);
    323     CHK(shtr_line_list_at(list2, iline, &lines2) == RES_OK);
    324 
    325     CHK(shtr_line_eq(&lines1, &lines2));
    326   }
    327 }
    328 
    329 static void
    330 test_serialization(struct shtr* shtr)
    331 {
    332   const struct shtr_line l[] = {
    333     {0.000134, 2.672E-38, 0.0533, 0.410, 608.4727, 0.79, 0.000060, 1, 4},
    334     {0.000379, 1.055E-39, 0.0418, 0.329,1747.9686, 0.79, 0.000110, 1, 5},
    335     {0.000448, 5.560E-38, 0.0490, 0.364,1093.0269, 0.79, 0.000060, 1, 4},
    336     {0.000686, 1.633E-36, 0.0578, 0.394, 701.1162, 0.79, 0.000180, 1, 4},
    337     {0.000726, 6.613E-33, 0.0695, 0.428, 402.3295, 0.79, 0.000240, 1, 3}
    338   };
    339   const size_t nlines = sizeof(l) / sizeof(struct shtr_line);
    340 
    341   struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL__;
    342   struct shtr_line_list* list1 = NULL;
    343   struct shtr_line_list* list2 = NULL;
    344 
    345   FILE* fp = NULL;
    346 
    347   CHK(args.file = tmpfile());
    348   print_lines(args.file, l, nlines);
    349   rewind(args.file);
    350 
    351   CHK(shtr_line_list_load(shtr, &args, &list1) == RES_OK);
    352   fclose(args.file);
    353 
    354   CHK(fp = tmpfile());
    355   CHK(shtr_line_list_write(NULL, fp) == RES_BAD_ARG);
    356   CHK(shtr_line_list_write(list1, NULL) == RES_BAD_ARG);
    357   CHK(shtr_line_list_write(list1, fp) == RES_OK);
    358   rewind(fp);
    359 
    360   CHK(shtr_line_list_create_from_stream(NULL, fp, &list2) == RES_BAD_ARG);
    361   CHK(shtr_line_list_create_from_stream(shtr, NULL, &list2) == RES_BAD_ARG);
    362   CHK(shtr_line_list_create_from_stream(shtr, fp, NULL) == RES_BAD_ARG);
    363   CHK(shtr_line_list_create_from_stream(shtr, fp, &list2) == RES_OK);
    364   fclose(fp);
    365 
    366   check_line_list_equality(list1, list2);
    367 
    368   CHK(shtr_line_list_ref_put(list1) == RES_OK);
    369   CHK(shtr_line_list_ref_put(list2) == RES_OK);
    370 }
    371 
    372 static void
    373 check_line(const struct shtr_line* ln)
    374 {
    375   /* Check NaN */
    376   CHK(ln->wavenumber == ln->wavenumber);
    377   CHK(ln->intensity == ln->intensity);
    378   CHK(ln->gamma_air == ln->gamma_air);
    379   CHK(ln->gamma_self == ln->gamma_self);
    380   CHK(ln->lower_state_energy == ln->lower_state_energy);
    381   CHK(ln->n_air == ln->n_air);
    382   CHK(ln->delta_air == ln->delta_air);
    383 
    384   CHK(ln->wavenumber > 0);
    385   CHK(ln->intensity > 0);
    386   CHK(ln->gamma_air >= 0);
    387   CHK(ln->gamma_self >= 0);
    388   CHK(ln->lower_state_energy >= 0);
    389   CHK(ln->molecule_id >= 0 && ln->molecule_id < 100);
    390   CHK(ln->isotope_id_local >= 0 && ln->isotope_id_local <= 9);
    391 }
    392 
    393 static void
    394 test_load_file(struct shtr* shtr, const char* path)
    395 {
    396   struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL__;
    397   struct shtr_line_list* list = NULL;
    398   size_t i, n;
    399   CHK(path);
    400   printf("Loading `%s'.\n", path);
    401   args.filename = path;
    402   CHK(shtr_line_list_load(shtr, &args, &list) == RES_OK);
    403   CHK(shtr_line_list_get_size(list, &n) == RES_OK);
    404   printf("  #lines: %lu\n", n);
    405 
    406   FOR_EACH(i, 0, n) {
    407     struct shtr_line line = SHTR_LINE_NULL;
    408     CHK(shtr_line_list_at(list, i, &line) == RES_OK);
    409     check_line(&line);
    410   }
    411 
    412   CHK(shtr_line_list_ref_put(list) == RES_OK);
    413 }
    414 
    415 int
    416 main(int argc, char** argv)
    417 {
    418   struct shtr_create_args args = SHTR_CREATE_ARGS_DEFAULT;
    419   struct shtr* shtr = NULL;
    420   int i;
    421   (void)argc, (void)argv;
    422 
    423   args.verbose = 1;
    424   CHK(shtr_create(&args, &shtr) == RES_OK);
    425 
    426   test_line_eq();
    427   test_load(shtr);
    428   test_load_failures(shtr);
    429   test_serialization(shtr);
    430 
    431   FOR_EACH(i, 1, argc) {
    432     char buf[64];
    433     struct time t0, t1;
    434 
    435     time_current(&t0);
    436     test_load_file(shtr, argv[i]);
    437     time_sub(&t0, time_current(&t1), &t0);
    438     time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
    439     printf("%s loaded in %s\n", argv[i], buf);
    440   }
    441 
    442   CHK(shtr_ref_put(shtr) == RES_OK);
    443   CHK(mem_allocated_size() == 0);
    444   return 0;
    445 }