test_sln_thermo_props.c (11622B)
1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2026 Université de Lorraine 3 * Copyright (C) 2022 Centre National de la Recherche Scientifique 4 * Copyright (C) 2022 Université Paul Sabatier 5 * 6 * This file is part of Star-Line. 7 * 8 * This program is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 #include <test_sln_lines.h> 22 23 #include "sln.h" 24 25 #include <rsys/mem_allocator.h> 26 #include <rsys/rsys.h> 27 28 struct thermo_props { 29 double xH2O; 30 double xCO2; 31 double xO3; 32 double pressure; /*[atm]*/ 33 double temperature; /*[K]*/ 34 }; 35 36 static const struct thermo_props thermo_props1 = {0.15, 0.10, 0.05, 10, 600}; 37 static const struct thermo_props thermo_props2 = {0.10, 0.15, 0.25, 20, 450}; 38 39 /******************************************************************************* 40 * Helper functions 41 ******************************************************************************/ 42 static struct shtr_isotope_metadata* 43 setup_isotopes(struct shtr* shtr) 44 { 45 struct shtr_isotope_metadata* metadata = NULL; 46 FILE* fp = NULL; 47 48 CHK(fp = tmpfile()); 49 fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n"); 50 write_shtr_molecule(fp, &g_H2O); 51 write_shtr_molecule(fp, &g_CO2); 52 write_shtr_molecule(fp, &g_O3); 53 rewind(fp); 54 55 CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &metadata) == RES_OK); 56 57 CHK(fclose(fp) == 0); 58 59 return metadata; 60 } 61 62 static struct shtr_line_list* 63 setup_lines(struct shtr* shtr) 64 { 65 struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL; 66 struct shtr_line_list* lines = NULL; 67 FILE* fp = NULL; 68 69 CHK(fp = tmpfile()); 70 write_shtr_lines(fp, g_lines, g_nlines); 71 rewind(fp); 72 73 args.filename = "stream"; 74 args.file = fp; 75 CHK(shtr_line_list_load(shtr, &args, &lines) == RES_OK); 76 77 CHK(fclose(fp) == 0); 78 79 return lines; 80 } 81 82 static struct sln_tree* 83 create_tree 84 (struct sln_device* sln, 85 struct shtr_isotope_metadata* mdata, 86 struct shtr_line_list* lines, 87 const struct thermo_props* props) 88 { 89 struct sln_tree_create_args args = SLN_TREE_CREATE_ARGS_DEFAULT; 90 struct sln_tree* tree = NULL; 91 92 args.metadata = mdata; 93 args.lines = lines; 94 95 args.molecules[SHTR_H2O].concentration = props->xH2O; 96 args.molecules[SHTR_H2O].cutoff = 25; /* [cm^-1] */ 97 args.molecules[SHTR_CO2].concentration = props->xCO2; 98 args.molecules[SHTR_CO2].cutoff = 50; /* [cm^-1] */ 99 args.molecules[SHTR_O3].concentration = props->xO3; 100 args.molecules[SHTR_O3].cutoff = 25; /* [cm^-1] */ 101 102 args.pressure = props->pressure; /*[atm]*/ 103 args.temperature = props->temperature; /*[K]*/ 104 105 CHK(sln_tree_create(sln, &args, &tree) == RES_OK); 106 return tree; 107 } 108 109 static INLINE double /* in [0,1[ */ 110 rand_canonic(void) 111 { 112 return (double)rand() / (double)((long)RAND_MAX+1); 113 } 114 115 static double /* [cm^-1] */ 116 line_sample_nu(const struct sln_tree* tree, const size_t iline) 117 { 118 struct sln_line line = SLN_LINE_NULL; 119 double nu_range[2] = {0,0}; /* [cm^-1] */ 120 double nu = 0; /* [cm^-1] */ 121 122 CHK(sln_tree_get_line(tree, iline, NULL, &line) == RES_OK); 123 nu_range[0] = line.wavenumber - 20; 124 nu_range[1] = line.wavenumber + 20; 125 nu = nu_range[0] + rand_canonic() * (nu_range[1] - nu_range[0]); 126 return nu; 127 } 128 129 static double /* [cm^-1] */ 130 node_sample_nu(const struct sln_tree* tree, const struct sln_node* node) 131 { 132 struct sln_line line = SLN_LINE_NULL; 133 struct sln_node_desc desc = SLN_NODE_DESC_NULL; 134 double nu_range[2] = {0,0}; /* [cm^-1] */ 135 double nu = 0; /* [cm^-1] */ 136 137 CHK(sln_node_get_desc(tree, node, &desc) == RES_OK); 138 139 CHK(sln_tree_get_line(tree, desc.ilines[0], NULL, &line) == RES_OK); 140 nu_range[0] = line.wavenumber - 10; 141 CHK(sln_tree_get_line(tree, desc.ilines[1], NULL, &line) == RES_OK); 142 nu_range[1] = line.wavenumber + 10; 143 144 nu = nu_range[0] + rand_canonic() * (nu_range[1] - nu_range[0]); 145 return nu; 146 } 147 148 static size_t 149 node_sample_line(const struct sln_tree* tree, const struct sln_node* node) 150 { 151 struct sln_node_desc desc = SLN_NODE_DESC_NULL; 152 const double r = rand_canonic(); 153 size_t iline = 0; 154 155 CHK(sln_node_get_desc(tree, node, &desc) == RES_OK); 156 iline = desc.ilines[0]; 157 iline += (size_t)(r * (double)(desc.ilines[1] - desc.ilines[0] + 1)); 158 return iline; 159 } 160 161 /* Check that, even belonging to 2 trees built from different thermodynamic 162 * properties, a line has the _exact_ same value when queried with the same 163 * thermodynamic properties. */ 164 static void 165 cmp_lines_values 166 (const struct sln_tree* tree1, 167 const struct sln_tree* tree2, 168 const size_t iline, 169 const double nu/*[cm^-1]*/) 170 { 171 struct sln_line line1 = SLN_LINE_NULL; 172 struct sln_line line2 = SLN_LINE_NULL; 173 struct sln_thermo_props props = SLN_THERMO_PROPS_NULL; 174 double ka1 = 0; 175 double ka2 = 0; 176 177 CHK(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK); 178 CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK); 179 ka1 = sln_line_eval(tree1, &line1, nu); 180 ka2 = sln_line_eval(tree2, &line2, nu); 181 CHK(ka1 != ka2); 182 183 props.concentrations[SHTR_H2O] = thermo_props1.xH2O; 184 props.concentrations[SHTR_CO2] = thermo_props1.xCO2; 185 props.concentrations[SHTR_O3] = thermo_props1.xO3; 186 props.pressure = thermo_props1.pressure; 187 props.temperature = thermo_props1.temperature; 188 CHK(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK); 189 CHK(sln_tree_get_line(tree2, iline, &props, &line2) == RES_OK); 190 ka1 = sln_line_eval(tree1, &line1, nu); 191 ka2 = sln_line_eval(tree2, &line2, nu); 192 CHK(ka1 == ka2); 193 194 props.concentrations[SHTR_H2O] = thermo_props2.xH2O; 195 props.concentrations[SHTR_CO2] = thermo_props2.xCO2; 196 props.concentrations[SHTR_O3] = thermo_props2.xO3; 197 props.pressure = thermo_props2.pressure; 198 props.temperature = thermo_props2.temperature; 199 CHK(sln_tree_get_line(tree1, iline, &props, &line1) == RES_OK); 200 CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK); 201 ka1 = sln_line_eval(tree1, &line1, nu); 202 ka2 = sln_line_eval(tree2, &line2, nu); 203 CHK(ka1 == ka2); 204 205 props.concentrations[SHTR_H2O] = thermo_props1.xH2O; 206 props.concentrations[SHTR_CO2] = thermo_props2.xCO2; 207 props.concentrations[SHTR_O3] = thermo_props1.xO3; 208 props.pressure = thermo_props2.pressure; 209 props.temperature = thermo_props1.temperature; 210 CHK(sln_tree_get_line(tree1, iline, &props, &line1) == RES_OK); 211 CHK(sln_tree_get_line(tree2, iline, &props, &line2) == RES_OK); 212 ka1 = sln_line_eval(tree1, &line1, nu); 213 ka2 = sln_line_eval(tree2, &line2, nu); 214 CHK(ka1 == ka2); 215 } 216 217 static void 218 cmp_lines 219 (const struct sln_tree* tree1, 220 const struct sln_tree* tree2, 221 const size_t iline) 222 { 223 const size_t N = 50; 224 size_t i = 0; 225 226 FOR_EACH(i, 0, N) { 227 const double nu = line_sample_nu(tree1, iline); 228 cmp_lines_values(tree1, tree2, iline, nu); 229 } 230 } 231 232 static void 233 cmp_nodes_lines 234 (const struct sln_tree* tree1, 235 const struct sln_node* node1, 236 const struct sln_tree* tree2, 237 const struct sln_node* node2) 238 { 239 const size_t N = 50; 240 size_t i = 0; 241 (void)node2; 242 243 FOR_EACH(i, 0, N) { 244 const size_t iline = node_sample_line(tree1, node1); 245 cmp_lines(tree1, tree2, iline); 246 } 247 } 248 249 /* Check that, even belonging to 2 trees built from different thermodynamic 250 * properties, a node has the _exact_ same value when evaluated with the same 251 * thermodynamic properties. */ 252 static void 253 cmp_nodes_values 254 (const struct sln_tree* tree1, 255 const struct sln_node* node1, 256 const struct sln_tree* tree2, 257 const struct sln_node* node2, 258 const double nu/*[cm^-1]*/) 259 { 260 struct sln_thermo_props props = SLN_THERMO_PROPS_NULL; 261 double ka1 = 0; 262 double ka2 = 0; 263 264 ka1 = sln_node_eval(tree1, node1, NULL, nu); 265 ka2 = sln_node_eval(tree2, node2, NULL, nu); 266 CHK(ka1 != ka2); 267 268 props.concentrations[SHTR_H2O] = thermo_props1.xH2O; 269 props.concentrations[SHTR_CO2] = thermo_props1.xCO2; 270 props.concentrations[SHTR_O3] = thermo_props1.xO3; 271 props.pressure = thermo_props1.pressure; 272 props.temperature = thermo_props1.temperature; 273 ka1 = sln_node_eval(tree1, node1, NULL, nu); 274 ka2 = sln_node_eval(tree2, node2, &props, nu); 275 CHK(ka1 == ka2); 276 277 props.concentrations[SHTR_H2O] = thermo_props2.xH2O; 278 props.concentrations[SHTR_CO2] = thermo_props2.xCO2; 279 props.concentrations[SHTR_O3] = thermo_props2.xO3; 280 props.pressure = thermo_props2.pressure; 281 props.temperature = thermo_props2.temperature; 282 ka1 = sln_node_eval(tree1, node1, &props, nu); 283 ka2 = sln_node_eval(tree2, node2, NULL, nu); 284 CHK(ka1 == ka2); 285 286 props.concentrations[SHTR_H2O] = thermo_props1.xH2O; 287 props.concentrations[SHTR_CO2] = thermo_props2.xCO2; 288 props.concentrations[SHTR_O3] = thermo_props1.xO3; 289 props.pressure = thermo_props2.pressure; 290 props.temperature = thermo_props1.temperature; 291 ka1 = sln_node_eval(tree1, node1, &props, nu); 292 ka2 = sln_node_eval(tree2, node2, &props, nu); 293 CHK(ka1 == ka2); 294 } 295 296 static void 297 cmp_nodes 298 (const struct sln_tree* tree1, 299 const struct sln_node* node1, 300 const struct sln_tree* tree2, 301 const struct sln_node* node2) 302 { 303 const size_t N = 50; 304 size_t i = 0; 305 306 FOR_EACH(i, 0, N) { 307 const double nu = node_sample_nu(tree1, node1); 308 cmp_nodes_values(tree1, node1, tree2, node2, nu); 309 } 310 311 cmp_nodes_lines(tree1, node1, tree2, node2); 312 } 313 314 static void 315 cmp_trees(const struct sln_tree* tree1, const struct sln_tree* tree2) 316 { 317 const struct sln_node* node1 = NULL; 318 const struct sln_node* node2 = NULL; 319 struct sln_tree_desc desc = SLN_TREE_DESC_NULL; 320 321 node1 = sln_tree_get_root(tree1); 322 node2 = sln_tree_get_root(tree2); 323 324 CHK(sln_tree_get_desc(tree1, &desc) == RES_OK); 325 CHK(desc.arity == 2); /* Assume that the arity of the tree is 2 */ 326 327 for(;;) { 328 unsigned ichild = 0; 329 330 cmp_nodes(tree1, node1, tree2, node2); 331 if(sln_node_is_leaf(node1)) break; 332 333 /* Randomly choose one node child */ 334 ichild = rand_canonic() < 0.5 ? 0 : 1; 335 node1 = sln_node_get_child(tree1, node1, ichild); 336 node2 = sln_node_get_child(tree2, node2, ichild); 337 } 338 } 339 340 /******************************************************************************* 341 * The test 342 ******************************************************************************/ 343 int 344 main(void) 345 { 346 struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT; 347 struct shtr* shtr = NULL; 348 struct shtr_isotope_metadata* mdata = NULL; 349 struct shtr_line_list* lines = NULL; 350 351 struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT; 352 struct sln_device* sln = NULL; 353 struct sln_tree* tree1 = NULL; 354 struct sln_tree* tree2 = NULL; 355 356 shtr_args.verbose = 1; 357 CHK(shtr_create(&shtr_args, &shtr) == RES_OK); 358 sln_args.verbose = 1; 359 CHK(sln_device_create(&sln_args, &sln) == RES_OK); 360 361 mdata = setup_isotopes(shtr); 362 lines = setup_lines(shtr); 363 364 tree1 = create_tree(sln, mdata, lines, &thermo_props1); 365 tree2 = create_tree(sln, mdata, lines, &thermo_props2); 366 367 cmp_trees(tree1, tree2); 368 369 CHK(shtr_ref_put(shtr) == RES_OK); 370 CHK(shtr_line_list_ref_put(lines) == RES_OK); 371 CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); 372 373 CHK(sln_device_ref_put(sln) == RES_OK); 374 CHK(sln_tree_ref_put(tree1) == RES_OK); 375 CHK(sln_tree_ref_put(tree2) == RES_OK); 376 377 CHK(mem_allocated_size() == 0); 378 return 0; 379 }