sln.h (16336B)
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 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 #ifndef SLN_H 20 #define SLN_H 21 22 #include <star/shtr.h> 23 #include <rsys/rsys.h> 24 25 #include <float.h> 26 #include <math.h> 27 28 /* Library symbol management */ 29 #if defined(SLN_SHARED_BUILD) /* Build shared library */ 30 #define SLN_API extern EXPORT_SYM 31 #elif defined(SLN_STATIC) /* Use/build static library */ 32 #define SLN_API extern LOCAL_SYM 33 #else 34 #define SLN_API extern IMPORT_SYM 35 #endif 36 37 /* Helper macro that asserts if the invocation of the sln function `Func' 38 * returns an error. One should use this macro on sln calls for which no 39 * explicit error checking is performed */ 40 #ifndef NDEBUG 41 #define SLN(Func) ASSERT(sln_ ## Func == RES_OK) 42 #else 43 #define SLN(Func) sln_ ## Func 44 #endif 45 46 #define SLN_TREE_DEPTH_MAX 64 /* Maximum depth of a tree */ 47 #define SLN_TREE_ARITY_MAX 256 /* Maximum arity of a tree */ 48 #define SLN_LEAF_NLINES_MAX 16384 /* Maximum number of lines per leaf */ 49 50 /* Forward declaration of external data structures */ 51 struct logger; 52 struct mem_allocator; 53 struct shtr; 54 struct shtr_line; 55 struct shtr_isotope_metadata; 56 struct shtr_line_list; 57 58 enum sln_mesh_type { 59 SLN_MESH_FIT, /* Fit the spectrum */ 60 SLN_MESH_UPPER, /* Upper limit of the spectrum */ 61 SLN_MESH_TYPES_COUNT__ 62 }; 63 64 enum sln_line_profile { 65 SLN_LINE_PROFILE_VOIGT, 66 SLN_LINE_PROFILES_COUNT__ 67 }; 68 69 struct sln_device_create_args { 70 struct logger* logger; /* May be NULL <=> default logger */ 71 struct mem_allocator* allocator; /* NULL <=> use default allocator */ 72 int verbose; /* Verbosity level */ 73 }; 74 #define SLN_DEVICE_CREATE_ARGS_DEFAULT__ {NULL,NULL,0} 75 static const struct sln_device_create_args SLN_DEVICE_CREATE_ARGS_DEFAULT = 76 SLN_DEVICE_CREATE_ARGS_DEFAULT__; 77 78 struct sln_isotope { 79 double abundance; /* in [0, 1] */ 80 int id; /* Identifier of the isotope */ 81 }; 82 83 struct sln_molecule { 84 struct sln_isotope isotopes[SHTR_MAX_ISOTOPE_COUNT]; 85 double concentration; 86 double cutoff; /* [cm^-1] */ 87 int non_default_isotope_abundances; 88 }; 89 #define SLN_MOLECULE_NULL__ {{{0}},0,0,0} 90 static const struct sln_molecule SLN_MOLECULE_NULL = SLN_MOLECULE_NULL__; 91 92 struct sln_tree_create_args { 93 /* Isotope metadata and list of spectral lines */ 94 struct shtr_isotope_metadata* metadata; 95 struct shtr_line_list* lines; 96 97 enum sln_line_profile line_profile; 98 /* Mixture description */ 99 struct sln_molecule molecules[SHTR_MAX_MOLECULE_COUNT]; 100 101 /* Thermo dynamic properties */ 102 double pressure; /* [atm] */ 103 double temperature; /* [K] */ 104 105 /* Hint on the number of vertices around the line center */ 106 size_t nvertices_hint; 107 108 /* Relative error used to simplify the spectrum mesh. The larger it is, the 109 * coarser the mesh */ 110 double mesh_decimation_err; /* > 0 */ 111 enum sln_mesh_type mesh_type; /* Type of mesh to generate */ 112 113 /* Maximum number of children per node */ 114 unsigned arity; 115 116 /* Maximum number of lines per leaf */ 117 unsigned leaf_nlines; 118 119 /* When this option is enabled, the polylines of internal nodes are 120 * constructed by merging their children's polylines in pairs (and then 121 * simplifying the result), and repeating the process until only a single 122 * polyline remains, which becomes the internal node's polyline. 123 * 124 * If this option is disabled, all child polylines are merged in a single step 125 * before being simplified. 126 * 127 * Enabling this option only makes sense for trees with an arity greater than 128 * two. For a binary tree, both methods should produce exactly the same tree, 129 * down to the bit */ 130 int collapse_polylines; 131 132 /* Advice on the number of threads to use */ 133 unsigned nthreads_hint; 134 }; 135 #define SLN_TREE_CREATE_ARGS_DEFAULT__ { \ 136 NULL, /* metadata */ \ 137 NULL, /* line list */ \ 138 SLN_LINE_PROFILE_VOIGT, /* Profile */ \ 139 {SLN_MOLECULE_NULL__}, /* Molecules */ \ 140 0, /* Pressure [atm] */ \ 141 0, /* Temperature [K] */ \ 142 16, /* #vertices hint */ \ 143 0.01f, /* Mesh decimation error */ \ 144 SLN_MESH_UPPER, /* Mesh type */ \ 145 2, /* Arity */ \ 146 1, /* Number of lines per leaf */ \ 147 0, /* Collapse polylines */ \ 148 (unsigned)(-1), /* #threads hint */ \ 149 } 150 static const struct sln_tree_create_args SLN_TREE_CREATE_ARGS_DEFAULT = 151 SLN_TREE_CREATE_ARGS_DEFAULT__; 152 153 struct sln_tree_read_args { 154 /* Metadata and list of spectral lines from which the tree was constructed */ 155 struct shtr_isotope_metadata* metadata; 156 struct shtr_line_list* lines; 157 158 /* Name of the file to read or of the provided stream. 159 * NULL <=> uses a default name for the stream to be read, which must 160 * therefore be defined. */ 161 const char* filename; /* Name of the file to read */ 162 FILE* file; /* Stream from where data are read. NULL <=> read from file */ 163 164 /* Verify that the digital signature of the input lines matches the one stored 165 * in the tree. In other words, ensure that this list of lines is indeed the 166 * one used to construct the tree. An error is returned if the signatures do 167 * not match. 168 * 169 * Although it is always advisable to verify that the data matches what is 170 * expected, calculating the signatures of the lines can be time-consuming. 171 * Therefore, a user who is _certain_ that the data matches can disable this 172 * verification */ 173 int disable_line_hash_check; 174 }; 175 #define SLN_TREE_READ_ARGS_NULL__ {NULL,NULL,NULL,NULL,0} 176 static const struct sln_tree_read_args SLN_TREE_READ_ARGS_NULL = 177 SLN_TREE_READ_ARGS_NULL__; 178 179 struct sln_tree_write_args { 180 /* Name of the file in which the tree is serialized. 181 * NULL <=> uses a default name for the stream to be written, which must 182 * therefore be defined. */ 183 const char* filename; /* Name of the file to read */ 184 185 /* Stream where data is written. 186 * NULL <=> write to the file defined by "filename" */ 187 FILE* file; 188 }; 189 #define SLN_TREE_WRITE_ARGS_NULL__ {NULL,NULL} 190 static const struct sln_tree_write_args SLN_TREE_WRITE_ARGS_NULL = 191 SLN_TREE_WRITE_ARGS_NULL__; 192 193 struct sln_tree_desc { 194 double mesh_decimation_err; 195 enum sln_mesh_type mesh_type; 196 enum sln_line_profile line_profile; 197 198 double pressure; /* [atm] */ 199 double temperature; /* [K] */ 200 201 unsigned depth; /* #edges from the root to the deepest leaf */ 202 size_t nlines; 203 size_t nvertices; 204 size_t nnodes; 205 unsigned arity; 206 unsigned leaf_nlines; 207 }; 208 #define SLN_TREE_DESC_NULL__ { \ 209 0,SLN_MESH_TYPES_COUNT__,SLN_LINE_PROFILES_COUNT__,0,0,0,0,0,0,0,0 \ 210 } 211 static const struct sln_tree_desc SLN_TREE_DESC_NULL = SLN_TREE_DESC_NULL__; 212 213 struct sln_node_desc { 214 /* Range of lines belonging to the node. The endpoints are included */ 215 size_t ilines[2]; 216 size_t nvertices; 217 unsigned nchildren; 218 }; 219 #define SLN_NODE_DESC_NULL__ {{0,0},0,0} 220 static const struct sln_node_desc SLN_NODE_DESC_NULL = SLN_NODE_DESC_NULL__; 221 222 struct sln_vertex { /* 8 Bytes */ 223 float wavenumber; /* in cm^-1 */ 224 float ka; 225 }; 226 #define SLN_VERTEX_NULL__ {0,0} 227 static const struct sln_vertex SLN_VERTEX_NULL = SLN_VERTEX_NULL__; 228 229 struct sln_mesh { 230 const struct sln_vertex* vertices; 231 size_t nvertices; 232 }; 233 #define SLN_MESH_NULL__ {NULL,0} 234 static const struct sln_mesh SLN_MESH_NULL = SLN_MESH_NULL__; 235 236 struct sln_mixture_load_args { 237 const char* filename; /* Name of the file to load or of the provided stream */ 238 FILE* file; /* Stream from where data are loaded. NULL <=> load from file */ 239 240 /* Metadata from which the mix is defined */ 241 struct shtr_isotope_metadata* molparam; 242 }; 243 #define SLN_MIXTURE_LOAD_ARGS_NULL__ {NULL,NULL,NULL} 244 static const struct sln_mixture_load_args SLN_MIXTURE_LOAD_ARGS_NULL = 245 SLN_MIXTURE_LOAD_ARGS_NULL__; 246 247 struct sln_line { 248 double wavenumber; /* Line center wrt pressure in cm^-1 */ 249 double profile_factor; /* m^-1.cm^-1 (1e2*density*intensity) */ 250 double gamma_d; /* Doppler half width */ 251 double gamma_l; /* Lorentz half width */ 252 enum shtr_molecule_id molecule_id; 253 }; 254 #define SLN_LINE_NULL__ {0,0,0,0,SHTR_MOLECULE_ID_NULL} 255 static const struct sln_line SLN_LINE_NULL = SLN_LINE_NULL__; 256 257 /* External data structure */ 258 struct ssp_rng; 259 260 /* Forward declarations of opaque data structures */ 261 struct sln_device; 262 struct sln_mixture; 263 struct sln_node; 264 struct sln_tree; 265 266 BEGIN_DECLS 267 268 /******************************************************************************* 269 * Device API 270 ******************************************************************************/ 271 SLN_API res_T 272 sln_device_create 273 (const struct sln_device_create_args* args, 274 struct sln_device** sln); 275 276 SLN_API res_T 277 sln_device_ref_get 278 (struct sln_device* sln); 279 280 SLN_API res_T 281 sln_device_ref_put 282 (struct sln_device* sln); 283 284 285 /******************************************************************************* 286 * Mixture API 287 ******************************************************************************/ 288 SLN_API res_T 289 sln_mixture_load 290 (struct sln_device* dev, 291 const struct sln_mixture_load_args* args, 292 struct sln_mixture** mixture); 293 294 SLN_API res_T 295 sln_mixture_ref_get 296 (struct sln_mixture* mixture); 297 298 SLN_API res_T 299 sln_mixture_ref_put 300 (struct sln_mixture* mixture); 301 302 SLN_API int 303 sln_mixture_get_molecule_count 304 (const struct sln_mixture* mixture); 305 306 SLN_API enum shtr_molecule_id 307 sln_mixture_get_molecule_id 308 (const struct sln_mixture* mixture, 309 const int index); 310 311 SLN_API res_T 312 sln_mixture_get_molecule 313 (const struct sln_mixture* mixture, 314 const int index, 315 struct sln_molecule* molecule); 316 317 /******************************************************************************* 318 * Tree API 319 ******************************************************************************/ 320 SLN_API res_T 321 sln_tree_create 322 (struct sln_device* dev, 323 const struct sln_tree_create_args* args, 324 struct sln_tree** tree); 325 326 /* Read a tree serialized with the "sln_tree_write" function */ 327 SLN_API res_T 328 sln_tree_read 329 (struct sln_device* sln, 330 const struct sln_tree_read_args* args, 331 struct sln_tree** tree); 332 333 SLN_API res_T 334 sln_tree_ref_get 335 (struct sln_tree* tree); 336 337 SLN_API res_T 338 sln_tree_ref_put 339 (struct sln_tree* tree); 340 341 SLN_API res_T 342 sln_tree_get_desc 343 (const struct sln_tree* tree, 344 struct sln_tree_desc* desc); 345 346 SLN_API const struct sln_node* /* NULL <=> No node */ 347 sln_tree_get_root 348 (const struct sln_tree* tree); 349 350 SLN_API res_T 351 sln_tree_get_line 352 (const struct sln_tree* tree, 353 const size_t iline, 354 struct sln_line* line); 355 356 SLN_API res_T 357 sln_tree_write 358 (const struct sln_tree* tree, 359 const struct sln_tree_write_args* args); 360 361 /******************************************************************************* 362 * Node API 363 ******************************************************************************/ 364 SLN_API int 365 sln_node_is_leaf 366 (const struct sln_node* node); 367 368 SLN_API unsigned 369 sln_node_get_child_count 370 (const struct sln_tree* tree, 371 const struct sln_node* node); 372 373 /* The node must not be a leaf */ 374 SLN_API const struct sln_node* 375 sln_node_get_child 376 (const struct sln_tree* tree, 377 const struct sln_node* node, 378 const unsigned ichild); /* 0 or #children */ 379 380 SLN_API double 381 sln_node_eval 382 (const struct sln_tree* tree, 383 const struct sln_node* node, 384 const double wavenumber); /* In cm^-1 */ 385 386 SLN_API res_T 387 sln_node_get_desc 388 (const struct sln_tree* tree, 389 const struct sln_node* node, 390 struct sln_node_desc* desc); 391 392 SLN_API res_T 393 sln_node_get_mesh 394 (const struct sln_tree* tree, 395 const struct sln_node* node, 396 struct sln_mesh* mesh); 397 398 /* Sample a leaf based on its importance, that is, its contribution to the 399 * node's absorption spectrum at the given wave number */ 400 SLN_API const struct sln_node* /* NULL <=> an error occurs */ 401 sln_node_sample_leaf 402 (const struct sln_tree* tree, 403 const struct sln_node* node, 404 const double nu, /* [cm^-1] */ 405 struct ssp_rng* rng, 406 double* proba); /* May be NULL */ 407 408 /******************************************************************************* 409 * Miscellaneous 410 ******************************************************************************/ 411 SLN_API double 412 sln_line_eval 413 (const struct sln_tree* tree, 414 const struct sln_line* line, 415 const double wavenumber); /* In cm^-1 */ 416 417 SLN_API double 418 sln_mesh_eval 419 (const struct sln_mesh* mesh, 420 const double wavenumber); /* In cm^-1 */ 421 422 /******************************************************************************* 423 * Helper functions 424 ******************************************************************************/ 425 /* Purpose: to calculate the Faddeeva function with relative error less than 426 * 10^(-4). 427 * 428 * Inputs: x and y, parameters for the Voigt function : 429 * - x is defined as x=(nu-nu_c)/gamma_D*sqrt(ln(2)) with nu the current 430 * wavenumber, nu_c the wavenumber at line center, gamma_D the Doppler 431 * linewidth. 432 * - y is defined as y=gamma_L/gamma_D*sqrt(ln(2)) with gamma_L the Lorentz 433 * linewith and gamma_D the Doppler linewidth 434 * 435 * Output: k, the Voigt function; it has to be multiplied by 436 * sqrt(ln(2)/pi)*1/gamma_D so that the result may be interpretable in terms of 437 * line profile. 438 * 439 * TODO check the copyright */ 440 SLN_API double 441 sln_faddeeva 442 (const double x, 443 const double y); 444 445 static INLINE double 446 sln_compute_line_half_width_doppler 447 (const double nu, /* Line center wrt pressure in cm^-1 */ /* TODO check this */ 448 const double molar_mass, /* In kg.mol^-1 */ 449 const double temperature) /* In K */ 450 { 451 /* kb = 1.3806e-23 452 * Na = 6.02214076e23 453 * c = 299792458 454 * sqrt(2*log(2)*kb*Na)/c */ 455 const double sqrt_two_ln2_kb_Na_over_c = 1.1324431552553545042e-08; 456 const double gamma_d = nu * sqrt_two_ln2_kb_Na_over_c * sqrt(temperature/molar_mass); 457 ASSERT(temperature >= 0 && molar_mass > 0); 458 return gamma_d; 459 } 460 461 static INLINE double 462 sln_compute_line_half_width_lorentz 463 (const double gamma_air, /* Air broadening half width [cm^-1.atm^-1] */ 464 const double gamma_self, /* Air broadening half width [cm^-1.atm^-1] */ 465 const double temperature, /* [K] */ 466 const double pressure, /* [atm^-1] */ 467 const double n_air, 468 const double concentration) 469 { 470 const double TREF=296; /* Ref temperature [K] for HITRAN/HITEMP database */ 471 const double Ps = pressure * concentration; 472 const double n_self = n_air; /* In HITRAN n_air == n_self */ 473 const double gamma_l = 474 pow(TREF/temperature, n_air) * (pressure - Ps) * gamma_air 475 + pow(TREF/temperature, n_self) * Ps * gamma_self; 476 ASSERT(gamma_air > 0 && gamma_self > 0); 477 ASSERT(pressure > 0 && concentration >= 0 && concentration <= 1); 478 479 return gamma_l; 480 } 481 482 static INLINE double 483 sln_compute_voigt_profile 484 (const double wavenumber, /* In cm^-1 */ 485 const double nu, /* Line center in cm^-1 */ 486 const double gamma_d, /* Doppler line half width in cm^-1 */ 487 const double gamma_l) /* Lorentz line half width in cm^-1 */ 488 { 489 /* Constants */ 490 const double sqrt_ln2 = 0.83255461115769768821; /* sqrt(log(2)) */ 491 const double sqrt_ln2_over_pi = 0.46971863934982566180; /* sqrt(log(2)/M_PI) */ 492 const double sqrt_ln2_over_gamma_d = sqrt_ln2 / gamma_d; 493 494 const double x = (wavenumber - nu) * sqrt_ln2_over_gamma_d; 495 const double y = gamma_l * sqrt_ln2_over_gamma_d; 496 const double k = sln_faddeeva(x, y); 497 return k*sqrt_ln2_over_pi/gamma_d; 498 } 499 500 static INLINE const char* 501 sln_mesh_type_cstr(const enum sln_mesh_type type) 502 { 503 const char* cstr = NULL; 504 505 switch(type) { 506 case SLN_MESH_FIT: cstr = "fit"; break; 507 case SLN_MESH_UPPER: cstr = "upper"; break; 508 default: FATAL("Unreachable code\n"); break; 509 } 510 return cstr; 511 } 512 513 static INLINE const char* 514 sln_line_profile_cstr(const enum sln_line_profile profile) 515 { 516 const char* cstr = NULL; 517 518 switch(profile) { 519 case SLN_LINE_PROFILE_VOIGT: cstr = "voigt"; break; 520 default: FATAL("Unreachable code\n"); break; 521 } 522 return cstr; 523 } 524 525 END_DECLS 526 527 #endif /* SLN_H */