sln.h (13582B)
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 49 /* Forward declaration of external data structures */ 50 struct logger; 51 struct mem_allocator; 52 struct shtr; 53 struct shtr_line; 54 struct shtr_isotope_metadata; 55 struct shtr_line_list; 56 57 enum sln_mesh_type { 58 SLN_MESH_FIT, /* Fit the spectrum */ 59 SLN_MESH_UPPER, /* Upper limit of the spectrum */ 60 SLN_MESH_TYPES_COUNT__ 61 }; 62 63 enum sln_line_profile { 64 SLN_LINE_PROFILE_VOIGT, 65 SLN_LINE_PROFILES_COUNT__ 66 }; 67 68 struct sln_device_create_args { 69 struct logger* logger; /* May be NULL <=> default logger */ 70 struct mem_allocator* allocator; /* NULL <=> use default allocator */ 71 int verbose; /* Verbosity level */ 72 }; 73 #define SLN_DEVICE_CREATE_ARGS_DEFAULT__ {NULL,NULL,0} 74 static const struct sln_device_create_args SLN_DEVICE_CREATE_ARGS_DEFAULT = 75 SLN_DEVICE_CREATE_ARGS_DEFAULT__; 76 77 struct sln_isotope { 78 double abundance; /* in [0, 1] */ 79 int id; /* Identifier of the isotope */ 80 }; 81 82 struct sln_molecule { 83 struct sln_isotope isotopes[SHTR_MAX_ISOTOPE_COUNT]; 84 double concentration; 85 double cutoff; /* [cm^-1] */ 86 int non_default_isotope_abundances; 87 }; 88 #define SLN_MOLECULE_NULL__ {{{0}},0,0,0} 89 static const struct sln_molecule SLN_MOLECULE_NULL = SLN_MOLECULE_NULL__; 90 91 struct sln_tree_create_args { 92 /* Isotope metadata and list of spectral lines */ 93 struct shtr_isotope_metadata* metadata; 94 struct shtr_line_list* lines; 95 96 enum sln_line_profile line_profile; 97 /* Mixture description */ 98 struct sln_molecule molecules[SHTR_MAX_MOLECULE_COUNT]; 99 100 /* Thermo dynamic properties */ 101 double pressure; /* [atm] */ 102 double temperature; /* [K] */ 103 104 /* Hint on the number of vertices around the line center */ 105 size_t nvertices_hint; 106 107 /* Relative error used to simplify the spectrum mesh. The larger it is, the 108 * coarser the mesh */ 109 double mesh_decimation_err; /* > 0 */ 110 enum sln_mesh_type mesh_type; /* Type of mesh to generate */ 111 112 /* Maximum number of children per node */ 113 unsigned arity; 114 }; 115 #define SLN_TREE_CREATE_ARGS_DEFAULT__ { \ 116 NULL, /* metadata */ \ 117 NULL, /* line list */ \ 118 SLN_LINE_PROFILE_VOIGT, /* Profile */ \ 119 {SLN_MOLECULE_NULL__}, /* Molecules */ \ 120 0, /* Pressure [atm] */ \ 121 0, /* Temperature [K] */ \ 122 16, /* #vertices hint */ \ 123 0.01f, /* Mesh decimation error */ \ 124 SLN_MESH_UPPER, /* Mesh type */ \ 125 2 /* Arity */ \ 126 } 127 static const struct sln_tree_create_args SLN_TREE_CREATE_ARGS_DEFAULT = 128 SLN_TREE_CREATE_ARGS_DEFAULT__; 129 130 struct sln_tree_read_args { 131 /* Metadata and list of spectral lines from which the tree was constructed */ 132 struct shtr_isotope_metadata* metadata; 133 struct shtr_line_list* lines; 134 135 /* Name of the file to read or of the provided stream. 136 * NULL <=> uses a default name for the stream to be read, which must 137 * therefore be defined. */ 138 const char* filename; /* Name of the file to read */ 139 FILE* file; /* Stream from where data are read. NULL <=> read from file */ 140 }; 141 #define SLN_TREE_READ_ARGS_NULL__ {NULL,NULL,NULL,NULL} 142 static const struct sln_tree_read_args SLN_TREE_READ_ARGS_NULL = 143 SLN_TREE_READ_ARGS_NULL__; 144 145 struct sln_tree_write_args { 146 /* Name of the file in which the tree is serialized. 147 * NULL <=> uses a default name for the stream to be written, which must 148 * therefore be defined. */ 149 const char* filename; /* Name of the file to read */ 150 151 /* Stream where data is written. 152 * NULL <=> write to the file defined by "filename" */ 153 FILE* file; 154 }; 155 #define SLN_TREE_WRITE_ARGS_NULL__ {NULL,NULL} 156 static const struct sln_tree_write_args SLN_TREE_WRITE_ARGS_NULL = 157 SLN_TREE_WRITE_ARGS_NULL__; 158 159 struct sln_tree_desc { 160 size_t max_nlines_per_leaf; 161 double mesh_decimation_err; 162 enum sln_mesh_type mesh_type; 163 enum sln_line_profile line_profile; 164 165 double pressure; /* [atm] */ 166 double temperature; /* [K] */ 167 168 unsigned depth; /* #edges from the root to the deepest leaf */ 169 size_t nlines; 170 size_t nvertices; 171 size_t nnodes; 172 unsigned arity; 173 }; 174 #define SLN_TREE_DESC_NULL__ { \ 175 0, 0, SLN_MESH_TYPES_COUNT__, SLN_LINE_PROFILES_COUNT__, 0, 0, 0, 0, 0, 0, 0 \ 176 } 177 static const struct sln_tree_desc SLN_TREE_DESC_NULL = SLN_TREE_DESC_NULL__; 178 179 struct sln_node_desc { 180 size_t nlines; 181 size_t nvertices; 182 unsigned nchildren; 183 }; 184 #define SLN_NODE_DESC_NULL__ {0,0,0} 185 static const struct sln_node_desc SLN_NODE_DESC_NULL = SLN_NODE_DESC_NULL__; 186 187 struct sln_vertex { /* 8 Bytes */ 188 float wavenumber; /* in cm^-1 */ 189 float ka; 190 }; 191 #define SLN_VERTEX_NULL__ {0,0} 192 static const struct sln_vertex SLN_VERTEX_NULL = SLN_VERTEX_NULL__; 193 194 struct sln_mesh { 195 const struct sln_vertex* vertices; 196 size_t nvertices; 197 }; 198 #define SLN_MESH_NULL__ {NULL,0} 199 static const struct sln_mesh SLN_MESH_NULL = SLN_MESH_NULL__; 200 201 struct sln_mixture_load_args { 202 const char* filename; /* Name of the file to load or of the provided stream */ 203 FILE* file; /* Stream from where data are loaded. NULL <=> load from file */ 204 205 /* Metadata from which the mix is defined */ 206 struct shtr_isotope_metadata* molparam; 207 }; 208 #define SLN_MIXTURE_LOAD_ARGS_NULL__ {NULL,NULL,NULL} 209 static const struct sln_mixture_load_args SLN_MIXTURE_LOAD_ARGS_NULL = 210 SLN_MIXTURE_LOAD_ARGS_NULL__; 211 212 /* Forward declarations of opaque data structures */ 213 struct sln_device; 214 struct sln_mixture; 215 struct sln_node; 216 struct sln_tree; 217 218 BEGIN_DECLS 219 220 /******************************************************************************* 221 * Device API 222 ******************************************************************************/ 223 SLN_API res_T 224 sln_device_create 225 (const struct sln_device_create_args* args, 226 struct sln_device** sln); 227 228 SLN_API res_T 229 sln_device_ref_get 230 (struct sln_device* sln); 231 232 SLN_API res_T 233 sln_device_ref_put 234 (struct sln_device* sln); 235 236 237 /******************************************************************************* 238 * Mixture API 239 ******************************************************************************/ 240 SLN_API res_T 241 sln_mixture_load 242 (struct sln_device* dev, 243 const struct sln_mixture_load_args* args, 244 struct sln_mixture** mixture); 245 246 SLN_API res_T 247 sln_mixture_ref_get 248 (struct sln_mixture* mixture); 249 250 SLN_API res_T 251 sln_mixture_ref_put 252 (struct sln_mixture* mixture); 253 254 SLN_API int 255 sln_mixture_get_molecule_count 256 (const struct sln_mixture* mixture); 257 258 SLN_API enum shtr_molecule_id 259 sln_mixture_get_molecule_id 260 (const struct sln_mixture* mixture, 261 const int index); 262 263 SLN_API res_T 264 sln_mixture_get_molecule 265 (const struct sln_mixture* mixture, 266 const int index, 267 struct sln_molecule* molecule); 268 269 /******************************************************************************* 270 * Tree API 271 ******************************************************************************/ 272 SLN_API res_T 273 sln_tree_create 274 (struct sln_device* dev, 275 const struct sln_tree_create_args* args, 276 struct sln_tree** tree); 277 278 /* Read a tree serialized with the "sln_tree_write" function */ 279 SLN_API res_T 280 sln_tree_read 281 (struct sln_device* sln, 282 const struct sln_tree_read_args* args, 283 struct sln_tree** tree); 284 285 SLN_API res_T 286 sln_tree_ref_get 287 (struct sln_tree* tree); 288 289 SLN_API res_T 290 sln_tree_ref_put 291 (struct sln_tree* tree); 292 293 SLN_API res_T 294 sln_tree_get_desc 295 (const struct sln_tree* tree, 296 struct sln_tree_desc* desc); 297 298 SLN_API const struct sln_node* /* NULL <=> No node */ 299 sln_tree_get_root 300 (const struct sln_tree* tree); 301 302 SLN_API int 303 sln_node_is_leaf 304 (const struct sln_node* node); 305 306 SLN_API unsigned 307 sln_node_get_child_count 308 (const struct sln_tree* tree, 309 const struct sln_node* node); 310 311 /* The node must not be a leaf */ 312 SLN_API const struct sln_node* 313 sln_node_get_child 314 (const struct sln_tree* tree, 315 const struct sln_node* node, 316 const unsigned ichild); /* 0 or #children */ 317 318 SLN_API size_t 319 sln_node_get_line_count 320 (const struct sln_node* node); 321 322 SLN_API res_T 323 sln_node_get_line 324 (const struct sln_tree* tree, 325 const struct sln_node* node, 326 const size_t iline, 327 struct shtr_line* line); 328 329 SLN_API res_T 330 sln_node_get_mesh 331 (const struct sln_tree* tree, 332 const struct sln_node* node, 333 struct sln_mesh* mesh); 334 335 SLN_API double 336 sln_node_eval 337 (const struct sln_tree* tree, 338 const struct sln_node* node, 339 const double wavenumber); /* In cm^-1 */ 340 341 SLN_API res_T 342 sln_node_get_desc 343 (const struct sln_tree* tree, 344 const struct sln_node* node, 345 struct sln_node_desc* desc); 346 347 SLN_API double 348 sln_mesh_eval 349 (const struct sln_mesh* mesh, 350 const double wavenumber); /* In cm^-1 */ 351 352 SLN_API res_T 353 sln_tree_write 354 (const struct sln_tree* tree, 355 const struct sln_tree_write_args* args); 356 357 /******************************************************************************* 358 * Helper functions 359 ******************************************************************************/ 360 /* Purpose: to calculate the Faddeeva function with relative error less than 361 * 10^(-4). 362 * 363 * Inputs: x and y, parameters for the Voigt function : 364 * - x is defined as x=(nu-nu_c)/gamma_D*sqrt(ln(2)) with nu the current 365 * wavenumber, nu_c the wavenumber at line center, gamma_D the Doppler 366 * linewidth. 367 * - y is defined as y=gamma_L/gamma_D*sqrt(ln(2)) with gamma_L the Lorentz 368 * linewith and gamma_D the Doppler linewidth 369 * 370 * Output: k, the Voigt function; it has to be multiplied by 371 * sqrt(ln(2)/pi)*1/gamma_D so that the result may be interpretable in terms of 372 * line profile. 373 * 374 * TODO check the copyright */ 375 SLN_API double 376 sln_faddeeva 377 (const double x, 378 const double y); 379 380 static INLINE double 381 sln_compute_line_half_width_doppler 382 (const double nu, /* Line center wrt pressure in cm^-1 */ /* TODO check this */ 383 const double molar_mass, /* In kg.mol^-1 */ 384 const double temperature) /* In K */ 385 { 386 /* kb = 1.3806e-23 387 * Na = 6.02214076e23 388 * c = 299792458 389 * sqrt(2*log(2)*kb*Na)/c */ 390 const double sqrt_two_ln2_kb_Na_over_c = 1.1324431552553545042e-08; 391 const double gamma_d = nu * sqrt_two_ln2_kb_Na_over_c * sqrt(temperature/molar_mass); 392 ASSERT(temperature >= 0 && molar_mass > 0); 393 return gamma_d; 394 } 395 396 static INLINE double 397 sln_compute_line_half_width_lorentz 398 (const double gamma_air, /* Air broadening half width [cm^-1.atm^-1] */ 399 const double gamma_self, /* Air broadening half width [cm^-1.atm^-1] */ 400 const double temperature, /* [K] */ 401 const double pressure, /* [atm^-1] */ 402 const double n_air, 403 const double concentration) 404 { 405 const double TREF=296; /* Ref temperature [K] for HITRAN/HITEMP database */ 406 const double Ps = pressure * concentration; 407 const double n_self = n_air; /* In HITRAN n_air == n_self */ 408 const double gamma_l = 409 pow(TREF/temperature, n_air) * (pressure - Ps) * gamma_air 410 + pow(TREF/temperature, n_self) * Ps * gamma_self; 411 ASSERT(gamma_air > 0 && gamma_self > 0); 412 ASSERT(pressure > 0 && concentration >= 0 && concentration <= 1); 413 414 return gamma_l; 415 } 416 417 static INLINE double 418 sln_compute_voigt_profile 419 (const double wavenumber, /* In cm^-1 */ 420 const double nu, /* Line center in cm^-1 */ 421 const double gamma_d, /* Doppler line half width in cm^-1 */ 422 const double gamma_l) /* Lorentz line half width in cm^-1 */ 423 { 424 /* Constants */ 425 const double sqrt_ln2 = 0.83255461115769768821; /* sqrt(log(2)) */ 426 const double sqrt_ln2_over_pi = 0.46971863934982566180; /* sqrt(log(2)/M_PI) */ 427 const double sqrt_ln2_over_gamma_d = sqrt_ln2 / gamma_d; 428 429 const double x = (wavenumber - nu) * sqrt_ln2_over_gamma_d; 430 const double y = gamma_l * sqrt_ln2_over_gamma_d; 431 const double k = sln_faddeeva(x, y); 432 return k*sqrt_ln2_over_pi/gamma_d; 433 } 434 435 static INLINE const char* 436 sln_mesh_type_cstr(const enum sln_mesh_type type) 437 { 438 const char* cstr = NULL; 439 440 switch(type) { 441 case SLN_MESH_FIT: cstr = "fit"; break; 442 case SLN_MESH_UPPER: cstr = "upper"; break; 443 default: FATAL("Unreachable code\n"); break; 444 } 445 return cstr; 446 } 447 448 static INLINE const char* 449 sln_line_profile_cstr(const enum sln_line_profile profile) 450 { 451 const char* cstr = NULL; 452 453 switch(profile) { 454 case SLN_LINE_PROFILE_VOIGT: cstr = "voigt"; break; 455 default: FATAL("Unreachable code\n"); break; 456 } 457 return cstr; 458 } 459 460 END_DECLS 461 462 #endif /* SLN_H */