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