sln.h (9361B)
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 77 struct sln_molecule { 78 double isotope_abundances[SLN_MAX_ISOTOPES_COUNT]; 79 double concentration; 80 double cutoff; /* [cm^-1] */ 81 int non_default_isotope_abundances; 82 }; 83 #define SLN_MOLECULE_NULL__ {{0},0,0,0} 84 static const struct sln_molecule SLN_MOLECULE_NULL = SLN_MOLECULE_NULL__; 85 86 struct sln_tree_create_args { 87 /* Isotope metadata and lise of spectral lines */ 88 struct shtr_isotope_metadata* metadata; 89 struct shtr_line_list* lines; 90 91 enum sln_line_profile line_profile; 92 /* Mixture description */ 93 struct sln_molecule molecules[SHTR_MAX_MOLECULE_COUNT]; 94 95 /* Thermo dynamic properties */ 96 double pressure; /* [atm] */ 97 double temperature; /* [K] */ 98 99 /* Hint on the number of vertices around the line center */ 100 size_t nvertices_hint; 101 102 /* Relative error used to simplify the spectrum mesh. The larger it is, the 103 * coarser the mesh */ 104 float mesh_decimation_err; /* > 0 */ 105 enum sln_mesh_type mesh_type; /* Type of mesh to generate */ 106 }; 107 #define SLN_TREE_CREATE_ARGS_DEFAULT__ { \ 108 NULL, /* metadata */ \ 109 NULL, /* line list */ \ 110 SLN_LINE_PROFILE_VOIGT, /* Profile */ \ 111 {SLN_MOLECULE_NULL__}, /* Molecules */ \ 112 0, /* Pressure [atm] */ \ 113 0, /* Temperature [K] */ \ 114 16, /* #vertices hint */ \ 115 0.01f, /* Mesh decimation error */ \ 116 SLN_MESH_UPPER, /* Mesh type */ \ 117 } 118 static const struct sln_tree_create_args SLN_TREE_CREATE_ARGS_DEFAULT = 119 SLN_TREE_CREATE_ARGS_DEFAULT__; 120 121 struct sln_tree_desc { 122 size_t max_nlines_per_leaf; 123 float mesh_decimation_err; 124 enum sln_mesh_type mesh_type; 125 enum sln_line_profile line_profile; 126 }; 127 #define SLN_TREE_DESC_NULL__ { \ 128 0, 0, SLN_MESH_TYPES_COUNT__, SLN_LINE_PROFILES_COUNT__ \ 129 } 130 static const struct sln_tree_desc SLN_TREE_DESC_NULL = SLN_TREE_DESC_NULL__; 131 132 struct sln_vertex { /* 8 Bytes */ 133 float wavenumber; /* in cm^-1 */ 134 float ka; 135 }; 136 #define SLN_VERTEX_NULL__ {0,0} 137 static const struct sln_vertex SLN_VERTEX_NULL = SLN_VERTEX_NULL__; 138 139 struct sln_mesh { 140 const struct sln_vertex* vertices; 141 size_t nvertices; 142 }; 143 #define SLN_MESH_NULL__ {NULL,0} 144 static const struct sln_mesh SLN_MESH_NULL = SLN_MESH_NULL__; 145 146 /* Forward declarations of opaque data structures */ 147 struct sln_device; 148 struct sln_node; 149 struct sln_tree; 150 151 /******************************************************************************* 152 * Device API 153 ******************************************************************************/ 154 SLN_API res_T 155 sln_device_create 156 (const struct sln_device_create_args* args, 157 struct sln_device** sln); 158 159 SLN_API res_T 160 sln_device_ref_get 161 (struct sln_device* sln); 162 163 SLN_API res_T 164 sln_device_ref_put 165 (struct sln_device* sln); 166 167 /******************************************************************************* 168 * Tree API 169 ******************************************************************************/ 170 SLN_API res_T 171 sln_tree_create 172 (struct sln_device* dev, 173 const struct sln_tree_create_args* args, 174 struct sln_tree** tree); 175 176 /* Load a tree serialized with the "sln_tree_write" function */ 177 SLN_API res_T 178 sln_tree_create_from_stream 179 (struct sln_device* sln, 180 struct shtr* shtr, 181 FILE* stream, 182 struct sln_tree** tree); 183 184 SLN_API res_T 185 sln_tree_ref_get 186 (struct sln_tree* tree); 187 188 SLN_API res_T 189 sln_tree_ref_put 190 (struct sln_tree* tree); 191 192 SLN_API res_T 193 sln_tree_get_desc 194 (const struct sln_tree* tree, 195 struct sln_tree_desc* desc); 196 197 SLN_API const struct sln_node* 198 sln_tree_get_root 199 (const struct sln_tree* tree); 200 201 SLN_API int 202 sln_node_is_leaf 203 (const struct sln_node* node); 204 205 /* Return NULL if the node is a leaf */ 206 SLN_API const struct sln_node* 207 sln_node_get_child 208 (const struct sln_node* node, 209 const unsigned ichild); /* 0 or 1 */ 210 211 SLN_API size_t 212 sln_node_get_lines_count 213 (const struct sln_node* node); 214 215 SLN_API res_T 216 sln_node_get_line 217 (const struct sln_tree* tree, 218 const struct sln_node* node, 219 const size_t iline, 220 struct shtr_line* line); 221 222 SLN_API res_T 223 sln_node_get_mesh 224 (const struct sln_tree* tree, 225 const struct sln_node* node, 226 struct sln_mesh* mesh); 227 228 SLN_API double 229 sln_node_eval 230 (const struct sln_tree* tree, 231 const struct sln_node* node, 232 const double wavenumber); /* In cm^-1 */ 233 234 SLN_API double 235 sln_mesh_eval 236 (const struct sln_mesh* mesh, 237 const double wavenumber); /* In cm^-1 */ 238 239 SLN_API res_T 240 sln_tree_write 241 (const struct sln_tree* tree, 242 FILE* stream); 243 244 /******************************************************************************* 245 * Helper functions 246 ******************************************************************************/ 247 /* Purpose: to calculate the Faddeeva function with relative error less than 248 * 10^(-4). 249 * 250 * Inputs: x and y, parameters for the Voigt function : 251 * - x is defined as x=(nu-nu_c)/gamma_D*sqrt(ln(2)) with nu the current 252 * wavenumber, nu_c the wavenumber at line center, gamma_D the Doppler 253 * linewidth. 254 * - y is defined as y=gamma_L/gamma_D*sqrt(ln(2)) with gamma_L the Lorentz 255 * linewith and gamma_D the Doppler linewidth 256 * 257 * Output: k, the Voigt function; it has to be multiplied by 258 * sqrt(ln(2)/pi)*1/gamma_D so that the result may be interpretable in terms of 259 * line profile. 260 * 261 * TODO check the copyright */ 262 SLN_API double 263 sln_faddeeva 264 (const double x, 265 const double y); 266 267 static INLINE double 268 sln_compute_line_half_width_doppler 269 (const double nu, /* Line center wrt pressure in cm^-1 */ /* TODO check this */ 270 const double molar_mass, /* In kg.mol^-1 */ 271 const double temperature) /* In K */ 272 { 273 /* kb = 1.3806e-23 274 * Na = 6.02214076e23 275 * c = 299792458 276 * sqrt(2*log(2)*kb*Na)/c */ 277 const double sqrt_two_ln2_kb_Na_over_c = 1.1324431552553545042e-08; 278 const double gamma_d = nu * sqrt_two_ln2_kb_Na_over_c * sqrt(temperature/molar_mass); 279 ASSERT(nu > 0 && temperature >= 0 && molar_mass > 0); 280 return gamma_d; 281 } 282 283 static INLINE double 284 sln_compute_line_half_width_lorentz 285 (const double gamma_air, /* Air broadening half width [cm^-1.atm^-1] */ 286 const double gamma_self, /* Air broadening half width [cm^-1.atm^-1] */ 287 const double pressure, /* [atm^-1] */ 288 const double concentration) 289 { 290 const double Ps = pressure * concentration; 291 const double gamma_l = (pressure - Ps) * gamma_air + Ps * gamma_self; 292 ASSERT(gamma_air > 0 && gamma_self > 0); 293 ASSERT(pressure > 0 && concentration >= 0 && concentration <= 1); 294 return gamma_l; 295 } 296 297 static INLINE double 298 sln_compute_voigt_profile 299 (const double wavenumber, /* In cm^-1 */ 300 const double nu, /* Line center in cm^-1 */ 301 const double gamma_d, /* Doppler line half width in cm^-1 */ 302 const double gamma_l) /* Lorentz line half width in cm^-1 */ 303 { 304 /* Constants */ 305 const double sqrt_ln2 = 0.83255461115769768821; /* sqrt(log(2)) */ 306 const double sqrt_ln2_over_pi = 0.46971863934982566180; /* sqrt(log(2)/M_PI) */ 307 const double sqrt_ln2_over_gamma_d = sqrt_ln2 / gamma_d; 308 309 const double x = (wavenumber - nu) * sqrt_ln2_over_gamma_d; 310 const double y = gamma_l * sqrt_ln2_over_gamma_d; 311 const double k = sln_faddeeva(x, y); 312 return k*sqrt_ln2_over_pi/gamma_d; 313 } 314 315 #endif /* SLN_H */