schiff.c (14429B)
1 /* Copyright (C) 2015, 2016, 2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2026 Clermont Auvergne INP 3 * Copyright (C) 2026 Institut Mines Télécom Albi-Carmaux 4 * Copyright (C) 2017, 2019-2021, 2026 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2026 Université de Lorraine 6 * Copyright (C) 2026 Université de Toulouse 7 * 8 * This program is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 #include "schiff_args.h" 22 #include "schiff_geometry.h" 23 #include "schiff_optical_properties.h" 24 25 #include <rsys/clock_time.h> 26 #include <rsys/stretchy_array.h> 27 28 #include <star/s3d.h> 29 #include <star/sschiff.h> 30 #include <star/ssp.h> 31 32 /******************************************************************************* 33 * Helper functions 34 ******************************************************************************/ 35 static res_T 36 dump_geometries 37 (const struct schiff_args* args, 38 struct s3d_device* s3d, 39 struct sschiff_geometry_distribution* distrib, 40 struct ssp_rng* rng, 41 FILE* stream) 42 { 43 struct s3d_scene* scn = NULL; 44 struct s3d_shape* shape = NULL; 45 unsigned i; 46 res_T res = RES_OK; 47 ASSERT(args && distrib && rng && stream); 48 49 res = s3d_scene_create(s3d, &scn); 50 if(res != RES_OK) { 51 fprintf(stderr, "Couldn't create the Star-3D scene.\n"); 52 goto error; 53 } 54 res = s3d_shape_create_mesh(s3d, &shape); 55 if(res != RES_OK) { 56 fprintf(stderr, "Couldn't create the Star-3D shape.\n"); 57 goto error; 58 } 59 res = s3d_scene_attach_shape(scn, shape); 60 if(res != RES_OK) { 61 fprintf(stderr, "Couldn't attach the Star-3D shape to the Star-3D scene.\n"); 62 goto error; 63 } 64 65 FOR_EACH(i, 0, args->ngeoms_dump) { 66 unsigned ivert, nverts; 67 unsigned itri, ntris; 68 double volume_scaling; 69 double dist_scaling; 70 void* shape_data = NULL; 71 res = distrib->sample(rng, &shape_data, shape, distrib->context); 72 if(res != RES_OK) { 73 fprintf(stderr, "Couldn't sample the micro organism geometry.\n"); 74 goto error; 75 } 76 res = distrib->sample_volume_scaling 77 (rng, shape_data, &volume_scaling, distrib->context); 78 if(res != RES_OK) { 79 fprintf(stderr, "Couldn't sample the volume scaling of the geometry.\n"); 80 goto error; 81 } 82 83 dist_scaling = pow(volume_scaling, 1.0/3.0); 84 85 fprintf(stream, "g schiff_geometry_%u\n", i); 86 87 /* Dump vertex position */ 88 fprintf(stream, "# List of vertex positions\n"); 89 S3D(mesh_get_vertices_count(shape, &nverts)); 90 FOR_EACH(ivert, 0, nverts) { 91 struct s3d_attrib attr; 92 S3D(mesh_get_vertex_attrib(shape, ivert, S3D_POSITION, &attr)); 93 f3_mulf(attr.value, attr.value, (float)dist_scaling); 94 fprintf(stream, "v %f %f %f\n", SPLIT3(attr.value)); 95 } 96 97 /* Dump triangle indices */ 98 fprintf(stream, "# Vertex indices of the triangular faces\n"); 99 S3D(mesh_get_triangles_count(shape, &ntris)); 100 FOR_EACH(itri, 0, ntris) { 101 unsigned ids[3]; 102 S3D(mesh_get_triangle_indices(shape, itri, ids)); 103 fprintf(stream, "f %u %u %u\n", ids[0]+1, ids[1]+1, ids[2]+1); 104 } 105 } 106 107 exit: 108 if(scn) S3D(scene_ref_put(scn)); 109 if(shape) S3D(shape_ref_put(shape)); 110 return res; 111 error: 112 goto exit; 113 } 114 115 static void 116 write_cross_sections 117 (FILE* stream, 118 struct sschiff_estimator* estimator, 119 const double* wlens, /* List of wavelengths in microns */ 120 const size_t nwlens) /* #wavelengths */ 121 { 122 size_t iwlen; 123 ASSERT(stream && estimator && wlens && nwlens); 124 125 FOR_EACH(iwlen, 0, nwlens) { 126 const struct sschiff_state* val; 127 struct sschiff_cross_section cross_section; 128 129 SSCHIFF(estimator_get_cross_section(estimator, iwlen, &cross_section)); 130 fprintf(stream, "%g ", wlens[iwlen]); 131 132 val = &cross_section.extinction; 133 fprintf(stream, "%g %g ", val->E, val->SE); 134 val = &cross_section.absorption; 135 fprintf(stream, "%g %g ", val->E, val->SE); 136 val = &cross_section.scattering; 137 fprintf(stream, "%g %g ", val->E, val->SE); 138 val = &cross_section.average_projected_area; 139 fprintf(stream, "%g %g\n", val->E, val->SE); 140 } 141 fprintf(stream, "\n"); 142 } 143 144 static void 145 write_phase_function_descriptors 146 (FILE* stream, 147 struct sschiff_estimator* estimator, 148 const double* wlens, /* List of wavelengths in microns */ 149 const size_t nwlens, /* #wavelenghts */ 150 const size_t nangles_inv) /* Number of inversed cumulative entries */ 151 { 152 size_t iwlen; 153 const double* angles; 154 size_t nangles; 155 res_T res = RES_OK; 156 ASSERT(stream && estimator && wlens && nwlens && nangles_inv > 1); 157 158 SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles)); 159 160 FOR_EACH(iwlen, 0, nwlens) { 161 size_t iangle; 162 struct sschiff_state cdf, pf; 163 double n; 164 165 res = sschiff_estimator_get_limit_scattering_angle_index 166 (estimator, iwlen, &iangle); 167 168 if(res == RES_BAD_OP) { 169 /* No valid limit angle => no phase function. Print -1 for all the data */ 170 fprintf(stderr, 171 "The phase function couldn't be computed.\n" 172 "wavelength = %g microns\n", wlens[iwlen]); 173 fprintf(stream, "-1 -1 -1 -1 -1 -1 -1 -1 -1\n"); 174 175 } else { 176 /* Retrieve the phase function descriptor data */ 177 ASSERT(res == RES_OK); 178 SSCHIFF(estimator_get_wide_scattering_angle_model_parameter 179 (estimator, iwlen, &n)); 180 SSCHIFF(estimator_get_differential_cross_section 181 (estimator, iwlen, iangle, &pf)); 182 SSCHIFF(estimator_get_differential_cross_section_cumulative 183 (estimator, iwlen, iangle, &cdf)); 184 185 /* Write the phase function descriptor */ 186 fprintf(stream, "%g %g %g %g %g %g %g %lu %lu\n", 187 wlens[iwlen], angles[iangle], pf.E, pf.SE, cdf.E, cdf.SE, n, 188 (unsigned long)nangles, 189 (unsigned long)nangles_inv); 190 } 191 } 192 fprintf(stream, "\n"); /* Empty line */ 193 } 194 195 static void 196 write_phase_functions 197 (FILE* stream, 198 const struct sschiff_estimator* estimator, 199 const double* wlens, 200 const size_t nwlens) 201 { 202 size_t iwlen; 203 const double* angles; 204 const struct sschiff_state* phase_funcs; 205 size_t iangle, nangles; 206 ASSERT(stream && estimator && wlens && nwlens); 207 (void)wlens; 208 209 SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles)); 210 FOR_EACH(iwlen, 0, nwlens) { 211 const res_T res = sschiff_estimator_get_phase_function 212 (estimator, iwlen, &phase_funcs); 213 if(res == RES_BAD_OP) { 214 /* The phase function couldn't be computed. Write default data. */ 215 FOR_EACH(iangle, 0, nangles) { 216 fprintf(stream, "-1 -1 -1\n"); 217 } 218 } else { 219 ASSERT(res == RES_OK); 220 FOR_EACH(iangle, 0, nangles) { 221 fprintf(stream, "%g %g %g\n", 222 angles[iangle], 223 phase_funcs[iangle].E, 224 phase_funcs[iangle].SE); 225 } 226 } 227 fprintf(stream, "\n"); 228 } 229 } 230 231 static void 232 write_phase_functions_cum 233 (FILE* stream, 234 struct sschiff_estimator* estimator, 235 const double* wlens, 236 const size_t nwlens) 237 { 238 size_t iwlen; 239 const double* angles; 240 const struct sschiff_state* phase_funcs_cum; 241 size_t iangle, nangles; 242 ASSERT(stream && estimator && wlens && nwlens); 243 (void)wlens; 244 245 SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles)); 246 FOR_EACH(iwlen, 0, nwlens) { 247 const res_T res = sschiff_estimator_get_phase_function_cumulative 248 (estimator, iwlen, &phase_funcs_cum); 249 if(res == RES_BAD_OP) { 250 /* The cumulative phase func couldn't be computed. Write default data */ 251 FOR_EACH(iangle, 0, nangles) { 252 fprintf(stream, "-1 -1 -1\n"); 253 } 254 } else { 255 ASSERT(res == RES_OK); 256 FOR_EACH(iangle, 0, nangles) { 257 fprintf(stream, "%g %g %g\n", 258 angles[iangle], 259 phase_funcs_cum[iangle].E, 260 phase_funcs_cum[iangle].SE); 261 } 262 } 263 fprintf(stream, "\n"); 264 } 265 } 266 267 static void 268 write_phase_functions_invcum 269 (FILE* stream, 270 struct sschiff_estimator* estimator, 271 const double* wlens, 272 const size_t nwlens, 273 const size_t nangles_inv) 274 { 275 double* invcum = NULL; 276 size_t iangle; 277 size_t iwlen; 278 double step; 279 ASSERT(stream && estimator && wlens && nwlens && nangles_inv > 1); 280 281 step = 1.0 / (double)(nangles_inv - 1); 282 283 if(!sa_add(invcum, nangles_inv)) { 284 fprintf(stderr, 285 "Couldn't allocate the list of inverse cumulative phase function angles.\n"); 286 287 /* Print the default -1 value */ 288 FOR_EACH(iwlen, 0, nwlens) { 289 FOR_EACH(iangle, 0, nangles_inv) { 290 fprintf(stream, "-1 -1\n"); 291 } 292 fprintf(stream, "\n"); 293 } 294 } else { 295 FOR_EACH(iwlen, 0, nwlens) { 296 const res_T res = sschiff_estimator_inverse_cumulative_phase_function 297 (estimator, iwlen, invcum, nangles_inv); 298 if(res != RES_OK) { 299 /* The inverse cumulative cannot be computed. Write -1 */ 300 fprintf(stderr, 301 "Couldn't inverse the cumulative phase function.\n" 302 "wavelength = %g microns\n", wlens[iwlen]); 303 FOR_EACH(iangle, 0, nangles_inv) { 304 fprintf(stream, "-1 -1\n"); 305 } 306 } else { 307 /* Write the inverse cumulative values */ 308 FOR_EACH(iangle, 0, nangles_inv-1) { 309 fprintf(stream, "%g %g\n", (double)iangle*step, invcum[iangle]); 310 } 311 /* Handle precision issues */ 312 fprintf(stream, "1 %g\n", invcum[nangles_inv-1]); 313 } 314 fprintf(stream, "\n"); /* Empty line */ 315 } 316 } 317 318 if(invcum) 319 sa_release(invcum); 320 } 321 322 static res_T 323 run_integration 324 (const struct schiff_args* args, 325 struct s3d_device* s3d, 326 struct sschiff_geometry_distribution* distrib, 327 struct ssp_rng* rng, 328 FILE* stream) 329 { 330 char dump[128]; 331 struct time t0, t1; 332 struct sschiff_device* sschiff = NULL; 333 struct sschiff_estimator* estimator = NULL; 334 size_t nwlens; 335 res_T res = RES_OK; 336 ASSERT(args && sa_size(args->wavelengths) && rng && stream); 337 338 res = sschiff_device_create(NULL, NULL, args->nthreads, 1, s3d, &sschiff); 339 if(res != RES_OK) { 340 fprintf(stderr, "Couldn't create the Star Schiff device.\n"); 341 goto error; 342 } 343 344 time_current(&t0); 345 346 /* Invoke the Schiff integration */ 347 nwlens = sa_size(args->wavelengths); 348 res = sschiff_integrate(sschiff, rng, distrib, args->wavelengths, nwlens, 349 sschiff_uniform_scattering_angles, args->nangles, args->nrealisations, 350 args->ninsamps, args->discard_large_angles, &estimator); 351 if(res != RES_OK) { 352 fprintf(stderr, "Schiff integration error.\n"); 353 goto error; 354 } 355 356 time_sub(&t0, time_current(&t1), &t0); 357 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 358 fprintf(stderr, "Computation performed in %s.\n", dump); 359 360 /* Write results */ 361 write_cross_sections(stream, estimator, args->wavelengths, nwlens); 362 write_phase_function_descriptors 363 (stream, estimator, args->wavelengths, nwlens, args->nangles_inv); 364 write_phase_functions(stream, estimator, args->wavelengths, nwlens); 365 write_phase_functions_cum(stream, estimator, args->wavelengths, nwlens); 366 write_phase_functions_invcum 367 (stream, estimator, args->wavelengths, nwlens, args->nangles_inv); 368 369 exit: 370 if(estimator) SSCHIFF(estimator_ref_put(estimator)); 371 if(sschiff) SSCHIFF(device_ref_put(sschiff)); 372 return res; 373 error: 374 goto exit; 375 } 376 377 static res_T 378 run(const struct schiff_args* args) 379 { 380 char dump[128]; 381 struct time t0, t1; 382 FILE* fp = stdout; 383 struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION; 384 struct ssp_rng* rng = NULL; 385 struct s3d_device* s3d = NULL; 386 res_T res = RES_OK; 387 ASSERT(args); 388 389 if(args->output_filename) { 390 fp = fopen(args->output_filename, "w"); 391 if(!fp) { 392 fprintf(stderr, 393 "Couldn't open the output file `%s'.\n", args->output_filename); 394 res = RES_IO_ERR; 395 goto error; 396 } 397 } 398 399 res = ssp_rng_create(&mem_default_allocator, SSP_RNG_THREEFRY, &rng); 400 if(res != RES_OK) { 401 fprintf(stderr, "Couldn't create the Random Number Generator.\n"); 402 goto error; 403 } 404 405 res = s3d_device_create(NULL, &mem_default_allocator, 0, &s3d); 406 if(res != RES_OK) { 407 fprintf(stderr, "Couldn't create the Star-3D device.\n"); 408 goto error; 409 } 410 411 time_current(&t0); 412 413 res = schiff_geometry_distribution_init(&distrib, s3d, args->geoms, 414 sa_size(args->geoms), args->characteristic_length, args->ran_geoms, 415 args->properties); 416 if(res != RES_OK) { 417 fprintf(stderr, 418 "Couldn't initialise the Star Schiff geometry distribution.\n"); 419 goto error; 420 } 421 422 time_sub(&t0, time_current(&t1), &t0); 423 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 424 fprintf(stderr, "Initialized the geometry distribution in %s.\n", dump); 425 426 if(!args->ngeoms_dump) { 427 res = run_integration(args, s3d, &distrib, rng, fp); 428 } else { 429 time_current(&t0); 430 res = dump_geometries(args, s3d, &distrib, rng, fp); 431 time_sub(&t0, time_current(&t1), &t0); 432 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 433 fprintf(stderr, "Dump geometries in %s.\n", dump); 434 } 435 436 /* Release before the check of the integration error code since if an error 437 * occurred the function will return without releasing the distribution. */ 438 schiff_geometry_distribution_release(&distrib); 439 440 if(res != RES_OK) 441 goto error; 442 443 exit: 444 if(fp && fp != stdout) fclose(fp); 445 if(rng) SSP(rng_ref_put(rng)); 446 if(s3d) S3D(device_ref_put(s3d)); 447 return res; 448 error: 449 goto exit; 450 } 451 452 /******************************************************************************* 453 * Entry point 454 ******************************************************************************/ 455 int 456 main(int argc, char** argv) 457 { 458 struct schiff_args args = SCHIFF_ARGS_NULL; 459 int err = 0; 460 res_T res; 461 462 res = schiff_args_init(&args, argc, argv); 463 if(res != RES_OK) goto error; 464 if(!args.ngeoms_dump && !args.wavelengths) goto exit; 465 res = run(&args); 466 if(res != RES_OK) goto error; 467 468 exit: 469 schiff_args_release(&args); 470 if(mem_allocated_size() != 0) { 471 fprintf(stderr, "Memory leaks: %lu Bytes\n", 472 (unsigned long) mem_allocated_size()); 473 } 474 return err; 475 error: 476 err = 1; 477 goto exit; 478 } 479