htrdr

Solving radiative transfer in heterogeneous media
git clone git://git.meso-star.fr/htrdr.git
Log | Files | Refs | README | LICENSE

commit 9d002482a618421a40eaf90931850c85b7b6903b
parent 7259298582ab56de48e9481ff3b4608284e88741
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 28 Oct 2022 15:55:35 +0200

htrdr-planeto: setup the spectral distribution

Diffstat:
Msrc/core/htrdr_cie_xyz.h | 2+-
Msrc/planeto/htrdr_planeto.c | 102++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---
Msrc/planeto/htrdr_planeto_c.h | 8++++++++
3 files changed, 108 insertions(+), 4 deletions(-)

diff --git a/src/core/htrdr_cie_xyz.h b/src/core/htrdr_cie_xyz.h @@ -34,7 +34,7 @@ HTRDR_CORE_API res_T htrdr_cie_xyz_create (struct htrdr* htrdr, const double range[2], /* Must be included in [380, 780] nanometers */ - const size_t nbands, /* # bands used to discretisze the CIE tristimulus s*/ + const size_t nbands, /* # bands used to discretize the CIE tristimulus */ struct htrdr_cie_xyz** cie); HTRDR_CORE_API void diff --git a/src/planeto/htrdr_planeto.c b/src/planeto/htrdr_planeto.c @@ -15,9 +15,11 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#define _POSIX_C_SOURCE 200112L /* fdopen */ +#define _POSIX_C_SOURCE 200112L /* fdopen, nextafter, rint */ #include "core/htrdr.h" +#include "core/htrdr_cie_xyz.h" +#include "core/htrdr_ran_wlen.h" #include "core/htrdr_log.h" #include "planeto/htrdr_planeto.h" @@ -35,12 +37,78 @@ #include <rsys/mem_allocator.h> #include <fcntl.h> /* open */ +#include <math.h> /* nextafter, rint */ #include <unistd.h> /* close */ #include <sys/stat.h> /******************************************************************************* * Helper function ******************************************************************************/ +/* Calculate the minimum length of the atmospheric spectral bands for the + * spectral domain considered */ +static double +compute_min_band_len(const struct htrdr_planeto* cmd) +{ + const double* range = NULL; /* In nm */ + double len = DBL_MAX; + size_t ibands[2]; + size_t i; + ASSERT(cmd); + + range = cmd->spectral_domain.wlen_range; + + /* The spectral range is degenerate to a wavelength */ + if(eq_eps(range[0], range[1], 1.e-6)) { + return 0; + } + + RNATM(find_bands(cmd->atmosphere, cmd->spectral_domain.wlen_range, ibands)); + + /* At least one band must be overlaped by the spectral domain */ + ASSERT(ibands[0]<=ibands[1]); + FOR_EACH(i, ibands[0], ibands[1]+1) { + double band_range[2]; + RNATM(band_get_range(cmd->atmosphere, i, band_range)); + + /* Make the upper bound inclusive */ + band_range[1] = nextafter(band_range[1], 0); + + /* Clamp the band range to the spectral domain */ + band_range[0] = MMAX(band_range[0], range[0]); + band_range[1] = MMIN(band_range[1], range[1]); + len = MMIN(band_range[1] - band_range[0], len); + } + return len; +} + +/* Calculate the number of fixed size spectral intervals to use for the + * cumulative */ +static size_t +compute_nintervals_for_spectral_cdf(const struct htrdr_planeto* cmd) +{ + double range_size; + double interval_len; + size_t nintervals; + ASSERT(cmd); + + range_size = + cmd->spectral_domain.wlen_range[1] + - cmd->spectral_domain.wlen_range[0]; + + /* Initially assume ~one interval per nanometer */ + nintervals = (size_t)rint(range_size); + + /* Calculate the minimum length of the atmospheric spectral bands fixed to + * the spectral integration domain. We ensure that an interval of the + * spectral cdf cannot be greater than this length */ + interval_len = compute_min_band_len(cmd); + if(interval_len < (range_size / (double)nintervals)) { + nintervals = (size_t)ceil(range_size / interval_len); + } + + return nintervals; +} + static res_T setup_octree_storage (struct htrdr_planeto* cmd, @@ -167,9 +235,36 @@ setup_spectral_domain (struct htrdr_planeto* cmd, const struct htrdr_planeto_args* args) { + size_t nintervals; + res_T res = RES_OK; ASSERT(cmd && args); + cmd->spectral_domain = args->spectral_domain; - return RES_OK; + + /* Configure the spectral distribution */ + nintervals = compute_nintervals_for_spectral_cdf(cmd); + switch(cmd->spectral_domain.spectral_type) { + /* Planck distribution */ + case HTRDR_SPECTRAL_LW: + case HTRDR_SPECTRAL_SW: + res = htrdr_ran_wlen_create(cmd->htrdr, cmd->spectral_domain.wlen_range, + nintervals, cmd->spectral_domain.ref_temperature, + &cmd->ran_wlen_planck); + break; + /* CIE XYZ distribution */ + case HTRDR_SPECTRAL_SW_CIE_XYZ: + res = htrdr_cie_xyz_create(cmd->htrdr, cmd->spectral_domain.wlen_range, + nintervals, &cmd->ran_wlen_cie); + break; + + default: FATAL("Unreachable code\n"); break; + } + if(res != RES_OK) goto error; + +exit: + return res; +error: + goto exit; } static res_T @@ -343,12 +438,13 @@ planeto_release(ref_T* ref) { struct htrdr_planeto* cmd = CONTAINER_OF(ref, struct htrdr_planeto, ref); struct htrdr* htrdr = NULL; - ASSERT(ref); if(cmd->atmosphere) RNATM(ref_put(cmd->atmosphere)); if(cmd->ground) RNGRD(ref_put(cmd->ground)); if(cmd->source) htrdr_planeto_source_ref_put(cmd->source); + if(cmd->ran_wlen_cie) htrdr_cie_xyz_ref_put(cmd->ran_wlen_cie); + if(cmd->ran_wlen_planck) htrdr_ran_wlen_ref_put(cmd->ran_wlen_planck); if(cmd->octrees_storage) CHK(fclose(cmd->octrees_storage) == 0); if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0); if(cmd->buf) htrdr_buffer_ref_put(cmd->buf); diff --git a/src/planeto/htrdr_planeto_c.h b/src/planeto/htrdr_planeto_c.h @@ -29,7 +29,9 @@ /* Forward declarations */ struct htrdr; +struct htrdr_cie; struct htrdr_pixel_format; +struct htrdr_ran_wlen; struct rnatm; struct rngrd; struct scam; @@ -68,6 +70,8 @@ struct htrdr_planeto { struct htrdr_planeto_source* source; struct htrdr_args_spectral spectral_domain; + struct htrdr_cie_xyz* ran_wlen_cie; /* HTRDR_SPECTRAL_SW_CIE_XYZ */ + struct htrdr_ran_wlen* ran_wlen_planck; /* HTRDR_SPECTRAL_<LW|SW> */ FILE* octrees_storage; @@ -85,6 +89,10 @@ struct htrdr_planeto { struct htrdr* htrdr; }; +extern LOCAL_SYM res_T +planeto_draw_map + (struct htrdr_planeto* cmd); + extern LOCAL_SYM void planeto_get_pixel_format (const struct htrdr_planeto* cmd,