sphin_spectral_property.c (10472B)
1 /* Copyright (C) 2024-2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2024-2026 Clermont Auvergne INP 3 * Copyright (C) 2024-2026 INSA Lyon 4 * Copyright (C) 2024-2026 Institut Mines Télécom Albi-Carmaux 5 * Copyright (C) 2024-2026 Institut National Polytechnique de Toulouse 6 * Copyright (C) 2024-2026 |Méso|Star> (contact@meso-star.com) 7 * Copyright (C) 2024-2026 PhotonLyX (info@photonlyx.com) 8 * Copyright (C) 2024-2026 Université de Lorraine 9 * Copyright (C) 2024-2026 Université Paul Sabatier 10 * Copyright (C) 2024-2026 Université Toulouse - Jean Jaurès 11 * 12 * This program is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU General Public License as published by 14 * the Free Software Foundation, either version 3 of the License, or 15 * (at your option) any later version. 16 * 17 * This program is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 * GNU General Public License for more details. 21 * 22 * You should have received a copy of the GNU General Public License 23 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 24 #define _POSIX_C_SOURCE 200112L /* for strtok_r support */ 25 26 #include "sphin.h" 27 #include "sphin_c.h" 28 #include "sphin_spectral_property.h" 29 30 #include <rsys/algorithm.h> 31 #include <rsys/cstr.h> 32 #include <rsys/mem_allocator.h> 33 #include <rsys/ref_count.h> 34 #include <rsys/text_reader.h> 35 36 #include <stdio.h> 37 #include <errno.h> 38 39 /******************************************************************************* 40 * Helper functions 41 ******************************************************************************/ 42 static res_T 43 spectral_property_create 44 (struct sphin* sphin, 45 const char* property_filename, 46 struct sphin_spectral_property** out_property) 47 { 48 struct sphin_spectral_property* property = NULL; 49 res_T res = RES_OK; 50 51 ASSERT(NULL != sphin); 52 ASSERT(NULL != out_property); 53 ASSERT(NULL != property_filename); 54 ASSERT('\0' != property_filename[0]); /* Name can't be empty */ 55 56 property = MEM_CALLOC(sphin->allocator, 1, 57 sizeof(struct sphin_spectral_property)); 58 if (NULL == property) { res = RES_MEM_ERR; goto error; } 59 60 /* Init geometry ref counter and init member variables */ 61 ref_init(&property->ref); 62 SPHIN(ref_get(sphin)); 63 property->sphin = sphin; 64 darray_double_init(sphin->allocator, &property->wavelengths); 65 darray_double_init(sphin->allocator, &property->values); 66 67 str_init(sphin->allocator, &property->property_filename); 68 res = str_set(&property->property_filename, property_filename); 69 if (RES_OK != res) { goto error; } 70 71 exit: 72 *out_property = property; 73 return res; 74 error: 75 if (NULL != property) { 76 SPHIN(spectral_property_ref_put(property)); 77 property = NULL; 78 } 79 goto exit; 80 } 81 82 static void 83 release_spectral_property 84 (ref_T* address) 85 { 86 struct sphin_spectral_property* property = NULL; 87 struct sphin* sphin = NULL; 88 89 ASSERT(NULL != address); 90 91 property = CONTAINER_OF(address, struct sphin_spectral_property, ref); 92 str_release(&property->property_filename); 93 darray_double_release(&property->wavelengths); 94 darray_double_release(&property->values); 95 sphin = property->sphin; 96 MEM_RM(sphin->allocator, property); 97 SPHIN(ref_put(sphin)); 98 } 99 100 /* As needed in search_lower_bound, see rsys/algorithm.h */ 101 static int 102 compare_wavelengths 103 (const void* key, 104 const void* element) 105 { 106 double target = *(const double*) key; 107 double array_element = *(const double*) element; 108 109 return (target > array_element) - (target < array_element); 110 } 111 112 /******************************************************************************* 113 * Local functions 114 ******************************************************************************/ 115 res_T 116 parse_spectral_property 117 (struct sphin* sphin, 118 char* filename, 119 struct sphin_spectral_property** out_property) 120 { 121 char* token = NULL; 122 char* token_ptr = NULL; 123 double wavelength = 0; 124 double wavelength_prev = 0; /* verify that the wavelengths are sorted */ 125 double value = 0; 126 127 FILE* stream = NULL; 128 struct sphin_spectral_property* property = NULL; 129 struct txtrdr* txtrdr = NULL; 130 struct str line; 131 res_T res = RES_OK; 132 133 ASSERT(NULL != sphin); 134 ASSERT(NULL != out_property); 135 ASSERT(NULL != filename); 136 ASSERT('\0' != filename[0]); /* filename can't be empty */ 137 138 str_init(sphin->allocator, &line); 139 140 stream = fopen(filename, "r"); 141 if (NULL == stream) { 142 ERROR(sphin, "Not possible to open file %s -- %s\n", 143 filename, strerror(errno)); 144 res = RES_IO_ERR; 145 goto error; 146 } 147 148 res = spectral_property_create(sphin, filename, &property); 149 if (RES_OK != res) { goto error; } 150 151 res = txtrdr_stream(sphin->allocator, stream, filename, '#', &txtrdr); 152 if (RES_OK != res) { goto error; } 153 154 res = txtrdr_read_line(txtrdr); 155 if (RES_OK != res) { goto error; } 156 157 while (NULL != txtrdr_get_line(txtrdr)) { 158 res = str_set(&line, txtrdr_get_cline(txtrdr)); 159 if (RES_OK != res) { goto error; } 160 161 /* Parse wavelength */ 162 token = strtok_r(str_get(&line), "\t ", &token_ptr); 163 if (NULL == token) { res = RES_BAD_ARG; goto error; } 164 if (NULL == token_ptr) { res = RES_BAD_ARG; goto error; } 165 166 res = cstr_to_double(token, &wavelength); 167 if (RES_OK != res) { 168 ERROR 169 (sphin, 170 "%s: %lu: could not read numerical value of wavelenght\n", 171 filename, txtrdr_get_line_num(txtrdr)); 172 goto error; } 173 174 if (wavelength_prev >= wavelength) { 175 res = RES_BAD_ARG; 176 ERROR 177 (sphin, 178 "%s: %lu: Wavelengths in spectral properties must be sorted\n", 179 filename, txtrdr_get_line_num(txtrdr)); 180 goto error; 181 } 182 183 res = darray_double_push_back(&property->wavelengths, &wavelength); 184 if (RES_OK != res) { goto error; } 185 186 res = cstr_to_double(token_ptr, &value); 187 if (RES_OK != res) { 188 ERROR 189 (sphin, 190 "%s: %lu: could not read numerical value of property\n", 191 filename, txtrdr_get_line_num(txtrdr)); 192 goto error; } 193 194 res = darray_double_push_back(&property->values, &value); 195 if (RES_OK != res) { goto error; } 196 197 wavelength_prev = wavelength; 198 199 res = txtrdr_read_line(txtrdr); 200 if (RES_OK != res) { goto error; } 201 } 202 203 exit: 204 str_release(&line); 205 if (NULL != stream) { 206 fclose(stream); 207 } 208 if (NULL != txtrdr) { 209 txtrdr_ref_put(txtrdr); 210 } 211 *out_property = property; 212 return res; 213 error: 214 if ( NULL != property) { 215 SPHIN(spectral_property_ref_put(property)); 216 property = NULL; 217 } 218 goto exit; 219 } 220 221 res_T 222 convert_wavelengths 223 (struct darray_double* wavelengths, 224 char* spectral_unit) 225 { 226 double scaling_factor = 1; 227 int invert = 0; 228 double* wls = NULL; 229 size_t i, wl_count = 0; 230 res_T res = RES_OK; 231 232 ASSERT(NULL != wavelengths); 233 ASSERT(NULL != spectral_unit); 234 235 if (0 == strcmp(spectral_unit, "nm")) {scaling_factor = 1e+00;} 236 else if (0 == strcmp(spectral_unit, "cm")) {scaling_factor = 1e07;} 237 else if (0 == strcmp(spectral_unit, "m")) {scaling_factor = 1e09;} 238 else if (0 == strcmp(spectral_unit, "cm^-1")) { 239 scaling_factor = 1e07; 240 invert = 1;} 241 else if (0 == strcmp(spectral_unit, "1/cm")) { 242 scaling_factor = 1e07; 243 invert = 1;} 244 else { res = RES_BAD_ARG; goto error; } 245 246 wl_count = darray_double_size_get(wavelengths); 247 wls = darray_double_data_get(wavelengths); 248 249 FOR_EACH(i, 0, wl_count) { 250 if (invert) { 251 if (0.0 == wls[i]) { res = RES_BAD_ARG; goto error; } 252 wls[i] = scaling_factor / wls[i]; 253 } 254 else {wls[i] = wls[i] * scaling_factor;} 255 } 256 257 exit: 258 return res; 259 error: 260 goto exit; 261 } 262 263 /******************************************************************************* 264 * Exported functions 265 ******************************************************************************/ 266 res_T 267 sphin_spectral_property_ref_get 268 (struct sphin_spectral_property* property) 269 { 270 if (NULL == property) { 271 return RES_BAD_ARG; 272 } 273 ref_get(&property->ref); 274 return RES_OK; 275 } 276 277 res_T 278 sphin_spectral_property_ref_put 279 (struct sphin_spectral_property* property) 280 { 281 if (NULL == property) { 282 return RES_BAD_ARG; 283 } 284 ref_put(&property->ref, release_spectral_property); 285 return RES_OK; 286 } 287 288 res_T 289 sphin_spectral_property_get_desc 290 (const struct sphin_spectral_property* property, 291 struct sphin_spectral_property_descriptor* desc) 292 { 293 if (NULL == property || NULL == desc) { 294 return RES_BAD_ARG; 295 } 296 297 desc->wavelengths = darray_double_data_get 298 ((struct darray_double*)&property->wavelengths); 299 desc->values = darray_double_data_get 300 ((struct darray_double*)&property->values); 301 desc->filename = (char*)str_get((struct str*)&property->property_filename); 302 303 desc->data_count = darray_double_size_get(&property->wavelengths); 304 305 return RES_OK; 306 } 307 308 res_T 309 sphin_spectral_property_interpolate_at_wavelength 310 (const struct sphin_spectral_property* property, 311 double wavelength, 312 enum sphin_interpolation_type interpolation_type, 313 double* value) 314 { 315 double wl_min, wl_max; /* Min and max wl values in the whole table */ 316 double wl_lower, wl_upper; /* Wl values used in interpolation */ 317 double v_lower, v_upper; /* Property values used in interpolation */ 318 const double* wavelengths = NULL; 319 const double* values = NULL; 320 const double* upper = NULL; 321 322 size_t data_count, upper_index, lower_index; 323 324 if (NULL == property || NULL == value) { 325 return RES_BAD_ARG; 326 } 327 328 if ((unsigned) interpolation_type >= SPHIN_INTERPOLATION_NONE__) { 329 return RES_BAD_ARG; 330 } 331 332 data_count = darray_double_size_get(&property->wavelengths); 333 wavelengths = darray_double_cdata_get(&property->wavelengths); 334 values = darray_double_cdata_get(&property->values); 335 wl_min = wavelengths[0]; 336 wl_max = wavelengths[data_count - 1]; 337 338 if (wavelength < wl_min || wavelength > wl_max) { 339 return RES_BAD_ARG; 340 } 341 342 upper = search_lower_bound 343 (&wavelength, wavelengths, data_count, sizeof(double), compare_wavelengths); 344 345 if (upper == wavelengths) { 346 /* Exact match is the first element of the wavelengths */ 347 *value = values[0]; 348 349 return RES_OK; 350 } 351 upper_index = (size_t)(upper - wavelengths); 352 lower_index = upper_index - 1; 353 354 wl_lower = wavelengths[lower_index]; 355 wl_upper = wavelengths[upper_index]; 356 v_lower = values[lower_index]; 357 v_upper = values[upper_index]; 358 359 *value = v_lower + (wavelength - wl_lower) * 360 (v_upper - v_lower) / (wl_upper - wl_lower); 361 362 return RES_OK; 363 }