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