star-phor

Radiative transfer solver for photoreactors.
git clone https://www.edstar.cnrs.fr/git/star-phor.git
Log | Files | Refs | README | LICENSE

commit b9bdd278372b200207ee9de1d61c1fceabcbea9f
parent 59b76e8ff5019cafa97baf2eeb4c65c604cb3712
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Mon,  8 Sep 2025 11:48:51 +0200

Implement spectral computations

The spectral properties implemented in star-phor-input provide the
foundation for wavelength-dependent computations in star-phor. This
commit introduces the structures and functions required for this type of
computation, as well as their integration with some of the functions
exposed by star-phor-input, in order to retrieve spectral properties
during path tracing.

For the data processing required by these computations, the normalized
discrete source emission spectrum provided by sphin is transformed into
a probability density function (PDF) using star-sampling. Each spectrum
is modeled as a piecewise linear function, and the resulting PDF is
included as part of the source_view structure.

Diffstat:
Msrc/sphor_compute_mvrea.c | 52++++++++++++++++++++++++++++++++++++++++++++++------
Msrc/sphor_sources.c | 56+++++++++++++++++++++++++++++++++++++++++++++++++++-----
Msrc/sphor_sources.h | 18++++++++++++++++--
3 files changed, 113 insertions(+), 13 deletions(-)

diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c @@ -44,13 +44,16 @@ compute_MVREA_realization double* weight/* TODO make this a list in sphor, since we can have multiple sensors. Each sensor should have its own weight */) { + double concentration = 0; double ka = 0; /* Absorption coefficient of current volume */ + double reflectivity = 0; /* Surface reflectivity */ + double response_function = 0; + double sigma_a = 0; + double wavelength = 0; double normal[3] = {0}; /* Geometry normal */ double dir[3] = {0}; /* Current direction in path sampling */ double pos[3] = {0}; /* Current position in path sampling */ double sample[3] = {0}; - double reflectivity = 0; /* Surface reflectivity */ - double response_function = 0; double s = 0; /* Random number */ enum sphin_brdf_type brdf_type = SPHIN_BRDF_NONE__; enum sphin_side side_id = SPHIN_SIDE_NONE__; @@ -61,8 +64,10 @@ compute_MVREA_realization float range[2] = {0.f, FLT_MAX}; float st[2] = {0}; /* Parametric coordinates of a point in a s3d_primitive */ int i = 0; /* Iterator */ + size_t j = 0; size_t triangle_id; size_t scn_prim_id; + size_t prop_rad_count; struct interface interface = INTERFACE_NULL; struct s3d_attrib attrib; struct s3d_hit hit = S3D_HIT_NULL; @@ -70,10 +75,15 @@ compute_MVREA_realization struct s3d_primitive src_prim = S3D_PRIMITIVE_NULL; struct source_view* source_view = NULL; struct sphin_brdf* brdf = NULL; + struct sphin_prop_rad* prop_rad = NULL; + struct sphin_scatterer* scatterer = NULL; struct sphin_sensor_volume* sensor_volume = NULL; struct sphin_source_surface* source = NULL; + struct sphin_spectral_property* abs_cross_sec = NULL; struct sphin_surface* surface = NULL; struct sphin_volume* volume = NULL; + struct sphin_source_surface_flux_density flux_density + = SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL; res_T res = RES_OK; ASSERT(NULL != sphor); @@ -94,8 +104,6 @@ compute_MVREA_realization d3_set_f3(pos, attrib.value); - /* TODO Sample a wavelength when spectral information available in sphin */ - res = sphin_config_get_surface (sphor->config, source_view->sphin_id, &surface); if (RES_OK != res){ goto error; } @@ -108,6 +116,15 @@ compute_MVREA_realization if (RES_OK != res){ goto error; } d3_normalize(normal, d3_minus(normal, d3_set_f3(normal, attrib.value))); + /* Sample a wavelength according to the source emission spectrum */ + res = sphin_source_surface_get_flux_density(source, &flux_density); + if (RES_OK != res){ goto error; } + + if (NULL != flux_density.emission_spectrum) { + wavelength = ssp_ranst_piecewise_linear_get + (source_view->emission_spectrum_pdf, rng); + } + /* Sample a direction according to the source direction distribution */ res = source_surface_sample_direction(source, rng, &src_prim, dir); if (RES_OK != res){ goto error; } @@ -157,10 +174,33 @@ compute_MVREA_realization (sphor->config, interface.volumes[side_id], &volume); if (RES_OK != res){ goto error; } - res = sphin_volume_get_ka(volume, &ka); + /* Get ka corresponding to the current wavelength */ + res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count); if (RES_OK != res){ goto error; } - /* Sample free path length according to the volume absorption coefficent */ + ka = 0.; + + FOR_EACH(j, 0, prop_rad_count) { + res = sphin_volume_get_prop_rad(volume, j, &prop_rad); + if (RES_OK != res){ goto error; } + + res = sphin_prop_rad_get_scatterer(prop_rad, &scatterer); + if (RES_OK != res){ goto error; } + + res = sphin_scatterer_get_concentration(scatterer, &concentration); + if (RES_OK != res){ goto error; } + + res = sphin_scatterer_get_abs_cross_sec(scatterer, &abs_cross_sec); + if (RES_OK != res){ goto error; } + + res = sphin_spectral_property_interpolate_at_wavelength + (abs_cross_sec, wavelength, SPHIN_INTERPOLATION_LINEAR, &sigma_a); + if (RES_OK != res){ goto error; } + + ka += concentration * sigma_a; + } + + /* Sample free path length according to volume absorption coefficent */ s = ssp_ran_exp(rng, ka);} else { s = FLT_MAX; diff --git a/src/sphor_sources.c b/src/sphor_sources.c @@ -100,6 +100,53 @@ error: } static res_T +setup_emission_spectrum_pdf + (struct sphor* sphor, + struct sphin_surface* surface, + struct source_view* source) +{ + struct sphin_source_surface* sphin_source = NULL; + struct sphin_spectral_property_descriptor spectrum = + SPHIN_SPECTRAL_PROPERTY_DESCRIPTOR_NULL; + struct sphin_source_surface_flux_density density = + SPHIN_SOURCE_SURFACE_FLUX_DENSITY_NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != surface); + ASSERT(NULL != source); + + SPHIN(surface_get_source(surface, &sphin_source)); + if(NULL == sphin_source) { res = RES_BAD_ARG; goto error; } + + SPHIN(source_surface_get_flux_density(sphin_source, &density)); + + /* TODO Discuss what to do if the spectrum is not given in the config file */ + if (NULL == density.emission_spectrum) { goto exit; } /* Nothing to do */ + + SPHIN(spectral_property_get_desc(density.emission_spectrum, &spectrum)); + + res = ssp_ranst_piecewise_linear_create + (sphor->allocator, &source->emission_spectrum_pdf); + if(RES_OK != res) { goto error; } + + res = ssp_ranst_piecewise_linear_setup + (source->emission_spectrum_pdf, + spectrum.wavelengths, + spectrum.values, + spectrum.data_count); + if(RES_OK != res) { goto error; } + +exit: + return res; +error: + if(source->emission_spectrum_pdf != NULL) { + SSP(ranst_piecewise_linear_ref_put(source->emission_spectrum_pdf)); + } + goto exit; +} + +static res_T setup_source_view (struct sphor* sphor, const struct sphor_create_args* args, @@ -132,6 +179,9 @@ setup_source_view res = setup_geometry_accel_struct(sphor, suniq, S3D_SAMPLE, &source->view); if (RES_OK != res) { goto error; } + res = setup_emission_spectrum_pdf(sphor, surface, source); + if (RES_OK != res) { goto error; } + exit: if (NULL != suniq) { SUNIQ(ref_put(suniq)); } return res; @@ -149,7 +199,6 @@ register_source_view struct sphin_surface* surface = NULL; struct sphin_source_surface* source_surface = NULL; struct source_view source = SOURCE_VIEW_NULL; - double power = 0; res_T res = RES_OK; ASSERT(NULL != sphor); @@ -161,9 +210,6 @@ register_source_view SPHIN(surface_get_source(surface, &source_surface)); ASSERT(NULL != source_surface); - res = sphin_surface_source_get_power(surface, &power); - if (RES_OK != res) { goto error; } - source.sphin_id = isurface; res = setup_source_view(sphor, args, surface, suniq_scene, &source); @@ -225,6 +271,7 @@ setup_source_distrib_power (struct sphor* sphor, const struct sphor_create_args* args) { + double power = 0; size_t i_surface = 0; size_t surface_count = 0; struct darray_double powers; @@ -250,7 +297,6 @@ setup_source_distrib_power if (RES_OK != res) { goto error; } if (NULL != source_surface) { - double power = 0; res = sphin_surface_source_get_power(surface, &power); if (RES_OK != res) { goto error; } res = darray_double_push_back(&powers, &power); diff --git a/src/sphor_sources.h b/src/sphor_sources.h @@ -25,6 +25,7 @@ #define SPHOR_SOURCES_H #include <star/s3d.h> +#include <star/ssp.h> #include <rsys/dynamic_array.h> #include <rsys/dynamic_array_size_t.h> @@ -37,8 +38,11 @@ struct ssp_rng; struct source_view { struct s3d_scene_view* view; /* View of the source */ - struct darray_size_t src2scn; /* Map the prim id of a source to its prim in - the scene */ + /* Map the prim id of a source to its prim in the scene */ + struct darray_size_t src2scn; + /* Normalized probability density function derived from emission spectrum */ + struct ssp_ranst_piecewise_linear* emission_spectrum_pdf; + size_t sphin_id; /* Surface identifier in star-phor-input */ }; @@ -67,6 +71,10 @@ source_release S3D(scene_view_ref_put(source->view)); source->view = NULL; } + if(NULL != source->emission_spectrum_pdf) { + SSP(ranst_piecewise_linear_ref_put(source->emission_spectrum_pdf)); + source->emission_spectrum_pdf = NULL; + } darray_size_t_release(&source->src2scn); } @@ -81,6 +89,10 @@ source_copy dst->sphin_id = src->sphin_id; S3D(scene_view_ref_get(src->view)); dst->view = src->view; + if(src->emission_spectrum_pdf != NULL) { + SSP(ranst_piecewise_linear_ref_get(src->emission_spectrum_pdf)); + dst->emission_spectrum_pdf = src->emission_spectrum_pdf; + } return darray_size_t_copy(&dst->src2scn, &src->src2scn); } @@ -94,8 +106,10 @@ source_copy_and_release dst->sphin_id = src->sphin_id; dst->view = src->view; + dst->emission_spectrum_pdf = src->emission_spectrum_pdf; src->view = NULL; + src->emission_spectrum_pdf = NULL; return darray_size_t_copy_and_release(&dst->src2scn, &src->src2scn); }