commit 885505948d3147cad6d73dfaf6114f66b02adef5
parent baa6ed5557c2e3641827012560512d0d258ec02c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 4 May 2020 12:09:48 +0200
Merge branch 'release_0.1'
Diffstat:
11 files changed, 686 insertions(+), 753 deletions(-)
diff --git a/README.md b/README.md
@@ -25,24 +25,18 @@ subsequent initialisation steps of the same sky. We point out that this file is
valid as long as the provided HTGOP, HTCP and HTMie files are the ones used to
build the cached structures. If not, an error is returned on sky creation.
-Note that this library provides functions compatible with the
-[Star-MTL](https://gitlab.Com/meso-star/star-mtl/) specification. Consequently,
-HTSky can be used as it, to describe a semi-transparent inhomogeneous medium in
-a Star-MTL file.
-
## How to build
This library is compatible GNU/Linux 64-bits. It relies on the
[CMake](http://www.cmake.org) and the
-[RCMake](https://gitlab.com/vaplv/rcmake/) packages to build. It also depends
-on the [HTCP](https://gitlab.com/meso-star/htcp/),
+[RCMake](https://gitlab.com/vaplv/rcmake/) packages to build.
+It also depends on the
+[HTCP](https://gitlab.com/meso-star/htcp/),
[HTGOP](https://gitlab.com/meso-star/htgop/),
[HTMie](https://gitlab.com/meso-star/htmie/),
[RSys](https://gitlab.com/vaplv/rsys/) and
[Star-VX](https://gitlab.com/meso-star/star-vx/) libraries, and on
-[OpenMP](http://www.openmp.org) 1.2 to parallelize its computations. If the
-[Star-MTL](https://gitlab.com/meso-star/star-mtl) library is found, HTSky
-will provide functions compatible with the Star-MTL specification.
+[OpenMP](http://www.openmp.org) 1.2 to parallelize its computations.
First ensure that CMake is installed on your system. Then install the RCMake
package as well as the aforementioned prerequisites. Finally generate the
@@ -52,6 +46,24 @@ resulting project can be edited, built, tested and installed as any CMake
project. Refer to the [CMake documentation](https://cmake.org/documentation)
for further informations on CMake.
+## Release notes
+
+### Version 0.1
+
+- Add the support of long waves. Add the `double wlen_lw_range[2]` member
+ variable to the `struct htsky_args` data structure that, once correctly
+ defined, is used to setup the sky data for the provided long wave range. By
+ default this range is degenerated meaning that the sky is setup for the short
+ wave range [380, 780] nm.
+- Add the `htsky_find_spectral_band` function: it returns the spectral band
+ that includes the submitted wavelength.
+- Remove the `htsky_sample_sw_spectral_data_CIE_1931_<X|Y|Z>` functions that
+ explicitly rely on the CIE XYZ color space.
+- Add the
+ `htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter` that
+ returns the Henyey-Greenstein phase function parameter for a given
+ wavelength.
+
## License
Copyright (C) 2018, 2019, 2020 [|Meso|Star](http://www.meso-star.com)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -21,13 +21,12 @@ set(HTSKY_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src)
################################################################################
# Check dependencies
################################################################################
-find_package(HTCP REQUIRED)
-find_package(HTGOP REQUIRED)
-find_package(HTMIE REQUIRED)
+find_package(HTCP 0.0.2 REQUIRED)
+find_package(HTGOP 0.1 REQUIRED)
+find_package(HTMIE 0.0.2 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
find_package(RCMake 0.3 REQUIRED)
find_package(RSys 0.7 REQUIRED)
-find_package(StarMTL)
find_package(StarVX 0.1 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
@@ -45,7 +44,7 @@ include_directories(
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 0)
+set(VERSION_MINOR 1)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -66,14 +65,6 @@ set(HTSKY_FILES_INC_API htsky.h)
set(HTSKY_FILES_DOC COPYING README.md)
-if(NOT StarMTL_FOUND)
- message(WARNING "Cannot find StarMTL. "
- "Build the HTSky library without the StarMTL binding.")
-else()
- set(HTSKY_FILES_SRC ${HTSKY_FILES_SRC} htsky_smtl.c)
- set(HTSKY_FILES_INC_API ${HTSKY_FILES_INC_API} htsky_smtl.h)
-endif()
-
# Prepend each file in the `HTSKY_FILES_<SRC|INC>' list by `HTSKY_SOURCE_DIR'
rcmake_prepend_path(HTSKY_FILES_SRC ${HTSKY_SOURCE_DIR})
rcmake_prepend_path(HTSKY_FILES_INC ${HTSKY_SOURCE_DIR})
diff --git a/src/htsky.c b/src/htsky.c
@@ -38,6 +38,11 @@
#include <omp.h>
+/* Current version the cache structure. One should increment it and perform a
+ * version management onto serialized data when the cache data data structure
+ * is updated. */
+static const int CACHE_VERSION = 0;
+
/*******************************************************************************
* Helper function
******************************************************************************/
@@ -55,72 +60,190 @@ check_args(const struct htsky_args* args)
}
static res_T
-setup_sw_bands_properties(struct htsky* sky)
+setup_bands_properties(struct htsky* sky)
{
- res_T res = RES_OK;
size_t nbands;
- size_t i;
+ size_t iband;
+ res_T res = RES_OK;
ASSERT(sky);
- nbands = htsky_get_sw_spectral_bands_count(sky);
+ nbands = htsky_get_spectral_bands_count(sky);
ASSERT(nbands);
- sky->sw_bands = MEM_CALLOC(sky->allocator, nbands, sizeof(*sky->sw_bands));
- if(!sky->sw_bands) {
- log_err(sky, "could not allocate the list of SW band properties.\n");
+
+ sky->bands = MEM_CALLOC(sky->allocator, nbands, sizeof(*sky->bands));
+ if(!sky->bands) {
+ log_err(sky, "Could not allocate the list of %s band properties.\n",
+ sky->is_long_wave ? "long wave" : "short wave");
res = RES_MEM_ERR;
goto error;
}
- FOR_EACH(i, 0, nbands) {
+ FOR_EACH(iband, sky->bands_range[0], sky->bands_range[1]+1) {
struct htgop_spectral_interval band;
double band_wlens[2];
+ const size_t i = iband - sky->bands_range[0];
- HTGOP(get_sw_spectral_interval(sky->htgop, i+sky->sw_bands_range[0], &band));
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
band_wlens[0] = wavenumber_to_wavelength(band.wave_numbers[1]);
band_wlens[1] = wavenumber_to_wavelength(band.wave_numbers[0]);
ASSERT(band_wlens[0] < band_wlens[1]);
- sky->sw_bands[i].Ca_avg = htmie_compute_xsection_absorption_average
+ sky->bands[i].Ca_avg = htmie_compute_xsection_absorption_average
(sky->htmie, band_wlens, HTMIE_FILTER_LINEAR);
- sky->sw_bands[i].Cs_avg = htmie_compute_xsection_scattering_average
+ sky->bands[i].Cs_avg = htmie_compute_xsection_scattering_average
(sky->htmie, band_wlens, HTMIE_FILTER_LINEAR);
- sky->sw_bands[i].g_avg = htmie_compute_asymmetry_parameter_average
+ sky->bands[i].g_avg = htmie_compute_asymmetry_parameter_average
(sky->htmie, band_wlens, HTMIE_FILTER_LINEAR);
- ASSERT(sky->sw_bands[i].Ca_avg > 0);
- ASSERT(sky->sw_bands[i].Cs_avg > 0);
- ASSERT(sky->sw_bands[i].g_avg > 0);
+ ASSERT(sky->bands[i].Ca_avg > 0);
+ ASSERT(sky->bands[i].Cs_avg > 0);
+ ASSERT(sky->bands[i].g_avg > 0);
}
exit:
return res;
error:
- if(sky->sw_bands) {
- MEM_RM(sky->allocator, sky->sw_bands);
- sky->sw_bands = NULL;
+ if(sky->bands) {
+ MEM_RM(sky->allocator, sky->bands);
+ sky->bands = NULL;
}
goto exit;
}
-static void
-sample_sw_spectral_data
- (struct htgop* htgop,
- const double r0,
- const double r1,
- res_T (*sample_sw_band)(const struct htgop*, const double, size_t*),
- size_t* ispectral_band,
- size_t* iquadrature_pt)
+static INLINE double
+particle_fetch_raw_property
+ (const struct htsky* sky,
+ const enum htsky_property prop,
+ const size_t iband,
+ const size_t iquad,
+ const size_t ivox[3])
{
- struct htgop_spectral_interval specint;
- res_T res = RES_OK;
- ASSERT(htgop && sample_sw_band && ispectral_band && iquadrature_pt);
- ASSERT(ispectral_band && iquadrature_pt);
- ASSERT(r0 >= 0 && r0 < 1);
- ASSERT(r1 >= 0 && r1 < 1);
- (void)res;
- res = sample_sw_band(htgop, r0, ispectral_band);
- ASSERT(res == RES_OK);
- HTGOP(get_sw_spectral_interval(htgop, *ispectral_band, &specint));
- HTGOP(spectral_interval_sample_quadrature(&specint, r1, iquadrature_pt));
+ double rho_da = 0; /* Dry air density */
+ double rct = 0; /* #droplets in kg of water per kg of dry air */
+ double ql = 0; /* Droplet density In kg.m^-3 */
+ double Ca = 0; /* Massic absorption cross section in m^2.kg^-1 */
+ double Cs = 0; /* Massic scattering cross section in m^2.kg^-1 */
+ double k_particle = 0;
+ size_t i = 0;
+ (void)iquad;
+ ASSERT(sky && ivox);
+
+ /* Compute he dry air density */
+ rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox);
+
+ /* Compute the droplet density */
+ rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
+ ql = rho_da * rct;
+
+ i = iband - sky->bands_range[0];
+
+ /* Use the average cross section of the current spectral band */
+ if(prop == HTSKY_Ka || prop == HTSKY_Kext) Ca = sky->bands[i].Ca_avg;
+ if(prop == HTSKY_Ks || prop == HTSKY_Kext) Cs = sky->bands[i].Cs_avg;
+
+ k_particle = ql*(Ca + Cs);
+ return k_particle;
+}
+
+static INLINE double
+gas_fetch_raw_property
+ (const struct htsky* sky,
+ const enum htsky_property prop,
+ const size_t iband,
+ const size_t iquad,
+ const int in_clouds,
+ const double pos[3],
+ const size_t ivox[3])
+{
+ struct htgop_layer layer;
+ size_t ilayer = 0;
+ double k_gas = 0;
+ ASSERT(sky && pos && ivox);
+
+ /* Retrieve the quadrature point into the spectral band of the layer into
+ * which `pos' lies */
+ HTGOP(position_to_layer_id(sky->htgop, pos[2], &ilayer));
+ HTGOP(get_layer(sky->htgop, ilayer, &layer));
+
+ if(sky->is_long_wave) {
+ struct htgop_layer_lw_spectral_interval band;
+ HTGOP(layer_get_lw_spectral_interval(&layer, iband, &band));
+
+ if(!in_clouds) {
+ /* Pos is outside the clouds. Directly fetch the nominal optical
+ * properties */
+ ASSERT(iquad < band.quadrature_length);
+ switch(prop) {
+ case HTSKY_Ka:
+ case HTSKY_Kext:
+ k_gas = band.ka_nominal[iquad];
+ break;
+ case HTSKY_Ks: k_gas = 0; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ } else {
+ /* Pos is inside the clouds. Compute the water vapor molar fraction at
+ * the current voxel */
+ const double x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
+ struct htgop_layer_lw_spectral_interval_tab tab;
+
+ /* Retrieve the tabulated data for the quadrature point */
+ HTGOP(layer_lw_spectral_interval_get_tab(&band, iquad, &tab));
+
+ /* Fetch the optical properties wrt x_h2o */
+ switch(prop) {
+ case HTSKY_Ka:
+ case HTSKY_Kext:
+ HTGOP(layer_lw_spectral_interval_tab_fetch_ka(&tab, x_h2o, &k_gas));
+ break;
+ case HTSKY_Ks: k_gas = 0; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+ } else {
+ struct htgop_layer_sw_spectral_interval band;
+ HTGOP(layer_get_sw_spectral_interval(&layer, iband, &band));
+
+ if(!in_clouds) {
+ /* Pos is outside the clouds. Directly fetch the nominal optical
+ * properties */
+ ASSERT(iquad < band.quadrature_length);
+ switch(prop) {
+ case HTSKY_Ka: k_gas = band.ka_nominal[iquad]; break;
+ case HTSKY_Ks: k_gas = band.ks_nominal[iquad]; break;
+ case HTSKY_Kext:
+ k_gas = band.ka_nominal[iquad] + band.ks_nominal[iquad];
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ } else {
+ /* Pos is inside the clouds. Compute the water vapor molar fraction at
+ * the current voxel */
+ const double x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
+ struct htgop_layer_sw_spectral_interval_tab tab;
+
+ /* Retrieve the tabulated data for the quadrature point */
+ HTGOP(layer_sw_spectral_interval_get_tab(&band, iquad, &tab));
+
+ /* Fetch the optical properties wrt x_h2o */
+ switch(prop) {
+ case HTSKY_Ka:
+ HTGOP(layer_sw_spectral_interval_tab_fetch_ka(&tab, x_h2o, &k_gas));
+ break;
+ case HTSKY_Ks:
+ HTGOP(layer_sw_spectral_interval_tab_fetch_ks(&tab, x_h2o, &k_gas));
+ break;
+ case HTSKY_Kext:
+ HTGOP(layer_sw_spectral_interval_tab_fetch_kext(&tab, x_h2o, &k_gas));
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+ }
+ return k_gas;
}
static res_T
@@ -152,7 +275,7 @@ setup_cache_stream
}
if(fd < 0) {
- log_err(sky, "unexpected error while opening the cache file `%s'.\n",
+ log_err(sky, "Unexpected error while opening the cache file `%s'.\n",
cache_filename);
res = RES_IO_ERR;
goto error;
@@ -160,7 +283,7 @@ setup_cache_stream
fp = fdopen(fd, "w+");
if(!fp) {
- log_err(sky, "could not open the cache file `%s'.\n", cache_filename);
+ log_err(sky, "Could not open the cache file `%s'.\n", cache_filename);
res = RES_IO_ERR;
goto error;
}
@@ -190,18 +313,24 @@ setup_cache_stream
goto error; \
} \
} (void)0
+ WRITE(&CACHE_VERSION, 1);
WRITE(&htcp_statbuf.st_ino, 1);
WRITE(&htcp_statbuf.st_mtim, 1);
WRITE(&htgop_statbuf.st_ino, 1);
WRITE(&htgop_statbuf.st_mtim, 1);
WRITE(&htmie_statbuf.st_ino, 1);
WRITE(&htmie_statbuf.st_mtim, 1);
+ WRITE(&sky->is_long_wave, 1);
+ WRITE(sky->bands_range, 2);
#undef WRITE
CHK(fflush(fp) == 0);
} else {
struct stat htcp_statbuf2;
struct stat htgop_statbuf2;
struct stat htmie_statbuf2;
+ int cache_version;
+ int is_long_wave;
+ size_t bands_range[2];
/* Read the cache header */
#define READ(Var, N) { \
@@ -217,12 +346,22 @@ setup_cache_stream
goto error; \
} \
} (void)0
+ READ(&cache_version, 1);
+ if(cache_version != CACHE_VERSION) {
+ log_err(sky,
+ "%s: invalid cache in version %d. Expecting a cache in version %d.\n",
+ cache_filename, cache_version, CACHE_VERSION);
+ res = RES_BAD_ARG;
+ goto error;
+ }
READ(&htcp_statbuf2.st_ino, 1);
READ(&htcp_statbuf2.st_mtim, 1);
READ(&htgop_statbuf2.st_ino, 1);
READ(&htgop_statbuf2.st_mtim, 1);
READ(&htmie_statbuf2.st_ino, 1);
READ(&htmie_statbuf2.st_mtim, 1);
+ READ(&is_long_wave, 1);
+ READ(bands_range, 2);
#undef READ
/* Compare the cache header with the input file status to check that the
@@ -241,6 +380,17 @@ setup_cache_stream
CHK_STAT(&htgop_statbuf, &htgop_statbuf2);
CHK_STAT(&htmie_statbuf, &htmie_statbuf2);
#undef CHK_STAT
+
+ /* Compare the handled spectral bands with the bands to handled to check
+ * that the cached octress are the expected ones */
+ if(is_long_wave != sky->is_long_wave
+ || bands_range[0] != sky->bands_range[0]
+ || bands_range[1] != sky->bands_range[1]) {
+ log_err(sky, "%s: invalid cache regarding the wavelengths to handle.\n",
+ cache_filename);
+ res = RES_BAD_ARG;
+ goto error;
+ }
}
exit:
@@ -258,6 +408,53 @@ error:
}
static void
+print_spectral_info(const struct htsky* sky)
+{
+ struct htgop_spectral_interval band_low, band_upp;
+ size_t iband_low, iband_upp;
+ size_t nbands;
+ size_t i;
+ size_t naccels = 0;
+ ASSERT(sky);
+
+ nbands = htsky_get_spectral_bands_count(sky);
+
+ iband_low = htsky_get_spectral_band_id(sky, 0);
+ iband_upp = htsky_get_spectral_band_id(sky, nbands-1);
+
+ /* Retrieve the spectral interval boundaries */
+ if(htsky_is_long_wave(sky)) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband_low, &band_low));
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband_upp, &band_upp));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband_low, &band_low));
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband_upp, &band_upp));
+ }
+
+ log_info(sky, "Sky data defined in [%g, %g] nanometers over %lu %s.\n",
+ wavenumber_to_wavelength(band_upp.wave_numbers[1]),
+ wavenumber_to_wavelength(band_low.wave_numbers[0]),
+ (unsigned long)nbands,
+ nbands > 1 ? "bands" : "band");
+
+ /* Compute the overall number of sky acceleration structures to build */
+ FOR_EACH(i, 0, nbands) {
+ struct htgop_spectral_interval band;
+ const size_t iband = htsky_get_spectral_band_id(sky, i);
+
+ if(htsky_is_long_wave(sky)) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
+ naccels += band.quadrature_length;
+ }
+
+ log_info(sky, "Number of clouds partitionning structures: %lu\n",
+ (unsigned long)naccels);
+}
+
+static void
release_sky(ref_T* ref)
{
struct htsky* sky;
@@ -269,7 +466,7 @@ release_sky(ref_T* ref)
if(sky->htcp) HTCP(ref_put(sky->htcp));
if(sky->htgop) HTGOP(ref_put(sky->htgop));
if(sky->htmie) HTMIE(ref_put(sky->htmie));
- if(sky->sw_bands) MEM_RM(sky->allocator, sky->sw_bands);
+ if(sky->bands) MEM_RM(sky->allocator, sky->bands);
darray_split_release(&sky->svx2htcp_z);
str_release(&sky->name);
ASSERT(MEM_ALLOCATED_SIZE(&sky->svx_allocator) == 0);
@@ -324,8 +521,8 @@ htsky_create
sky->is_cloudy = args->htcp_filename != NULL;
darray_split_init(sky->allocator, &sky->svx2htcp_z);
str_init(sky->allocator, &sky->name);
- sky->sw_bands_range[0] = 1;
- sky->sw_bands_range[1] = 0;
+ sky->bands_range[0] = 1;
+ sky->bands_range[1] = 0;
sky->nthreads = MMIN(args->nthreads, (unsigned)nthreads_max);
if(logger) {
@@ -336,14 +533,14 @@ htsky_create
res = str_set(&sky->name, args->name);
if(res != RES_OK) {
- log_err(sky, "cannot setup the material name to `%s'.\n", args->name);
+ log_err(sky, "Cannot setup the material name to `%s'.\n", args->name);
goto error;
}
/* Setup an allocator specific to the SVX library */
res = mem_init_proxy_allocator(&sky->svx_allocator, sky->allocator);
if(res != RES_OK) {
- log_err(sky, "cannot init the allocator used to manage the Star-VX data.\n");
+ log_err(sky, "Cannot init the allocator used to manage the Star-VX data.\n");
goto error;
}
@@ -351,26 +548,40 @@ htsky_create
res = svx_device_create
(sky->logger, &sky->svx_allocator, sky->verbose, &sky->svx);
if(res != RES_OK) {
- log_err(sky, "error creating the Star-VX library device.\n");
+ log_err(sky, "Error creating the Star-VX library device.\n");
goto error;
}
/* Load the gas optical properties */
res = htgop_create(sky->logger, sky->allocator, sky->verbose, &sky->htgop);
if(res != RES_OK) {
- log_err(sky, "could not create the gas optical properties loader.\n");
+ log_err(sky, "Could not create the gas optical properties loader.\n");
goto error;
}
res = htgop_load(sky->htgop, args->htgop_filename);
if(res != RES_OK) {
- log_err(sky, "error loading the gas optical properties -- `%s'.\n",
+ log_err(sky, "Error loading the gas optical properties -- `%s'.\n",
args->htgop_filename);
goto error;
}
+ /* If the long wave range is degenerated use the default wavelength range,
+ * i.e. the short wave wavelengths of the CIE XYZ color space */
+ if(args->wlen_lw_range[0] > args->wlen_lw_range[1]) {
+ double wnums[2];
+ wnums[0] = wavelength_to_wavenumber(780.0);
+ wnums[1] = wavelength_to_wavenumber(380.0);
+ res = htgop_get_sw_spectral_intervals(sky->htgop, wnums, sky->bands_range);
+ if(res != RES_OK) goto error;
+ } else {
+ double wnums[2];
+ wnums[0] = wavelength_to_wavenumber(args->wlen_lw_range[1]);
+ wnums[1] = wavelength_to_wavenumber(args->wlen_lw_range[0]);
+ res = htgop_get_lw_spectral_intervals(sky->htgop, wnums, sky->bands_range);
+ if(res != RES_OK) goto error;
+ sky->is_long_wave = 1;
+ }
- /* Fetch short wave bands range */
- res = htgop_get_sw_spectral_intervals_CIE_XYZ(sky->htgop, sky->sw_bands_range);
- if(res != RES_OK) goto error;
+ print_spectral_info(sky);
/* Setup the atmopshere */
time_current(&t0);
@@ -378,19 +589,19 @@ htsky_create
if(res != RES_OK) goto error;
time_sub(&t0, time_current(&t1), &t0);
time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
- log_info(sky, "setup atmosphere in %s\n", buf);
+ log_info(sky, "Setup atmosphere in %s\n", buf);
/* Nothing more to do */
if(!sky->is_cloudy) goto exit;
if(!args->htmie_filename) {
- log_err(sky, "missing the HTMie filename.\n");
+ log_err(sky, "Missing the HTMie filename.\n");
res = RES_BAD_ARG;
goto error;
}
if(!args->htcp_filename) {
- log_err(sky, "missing the HTCP filename.\n");
+ log_err(sky, "Missing the HTCP filename.\n");
res = RES_BAD_ARG;
goto error;
}
@@ -398,27 +609,28 @@ htsky_create
/* Load MIE data */
res = htmie_create(sky->logger, sky->allocator, sky->verbose, &sky->htmie);
if(res != RES_OK) {
- log_err(sky, "could not create the Mie's data loader.\n");
+ log_err(sky, "Could not create the Mie's data loader.\n");
goto error;
}
res = htmie_load(sky->htmie, args->htmie_filename);
if(res != RES_OK) {
- log_err(sky, "error loading the Mie's data -- `%s'.\n", args->htmie_filename);
+ log_err(sky, "Error loading the Mie's data -- `%s'.\n", args->htmie_filename);
goto error;
}
- res = setup_sw_bands_properties(sky);
+ /* Setup the properties of the Short/Long wave bands */
+ res = setup_bands_properties(sky);
if(res != RES_OK) goto error;
/* Load clouds properties */
res = htcp_create(sky->logger, sky->allocator, sky->verbose, &sky->htcp);
if(res != RES_OK) {
- log_err(sky, "could not create the loader of cloud properties.\n");
+ log_err(sky, "Could not create the loader of cloud properties.\n");
goto error;
}
res = htcp_load(sky->htcp, args->htcp_filename);
if(res != RES_OK) {
- log_err(sky, "error loading the cloud properties -- `%s'.\n",
+ log_err(sky, "Error loading the cloud properties -- `%s'.\n",
args->htcp_filename);
goto error;
}
@@ -435,7 +647,7 @@ htsky_create
if(res != RES_OK) goto error;
time_sub(&t0, time_current(&t1), &t0);
time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
- log_info(sky, "setup clouds in %s\n", buf);
+ log_info(sky, "Setup clouds in %s\n", buf);
if(sky->verbose) {
log_svx_memory_usage(sky);
@@ -484,14 +696,28 @@ htsky_fetch_particle_phase_function_asymmetry_parameter
{
size_t i;
ASSERT(sky);
- ASSERT(ispectral_band >= sky->sw_bands_range[0]);
- ASSERT(ispectral_band <= sky->sw_bands_range[1]);
+ ASSERT(ispectral_band >= sky->bands_range[0]);
+ ASSERT(ispectral_band <= sky->bands_range[1]);
(void)iquad;
if(!sky->is_cloudy) {
return 0;
} else {
- i = ispectral_band - sky->sw_bands_range[0];
- return sky->sw_bands[i].g_avg;
+ i = ispectral_band - sky->bands_range[0];
+ return sky->bands[i].g_avg;
+ }
+}
+
+double
+htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter
+ (const struct htsky* sky,
+ const double wavelength) /* In nanometer */
+{
+ ASSERT(sky);
+ if(!sky->is_cloudy) {
+ return 0;
+ } else {
+ return htmie_fetch_asymmetry_parameter
+ (sky->htmie, wavelength, HTMIE_FILTER_LINEAR);
}
}
@@ -518,11 +744,11 @@ htsky_fetch_raw_property
double k_gas = 0;
double k = 0;
ASSERT(sky && pos);
- ASSERT(iband >= sky->sw_bands_range[0]);
- ASSERT(iband <= sky->sw_bands_range[1]);
+ ASSERT(iband >= sky->bands_range[0]);
+ ASSERT(iband <= sky->bands_range[1]);
ASSERT(comp_mask & HTSKY_CPNT_MASK_ALL);
- i = iband - sky->sw_bands_range[0];
+ i = iband - sky->bands_range[0];
cloud_desc = sky->is_cloudy ? &sky->clouds[i][iquad].octree_desc : NULL;
atmosphere_desc = &sky->atmosphere[i][iquad].bitree_desc;
ASSERT(atmosphere_desc->frame[0] == SVX_AXIS_Z);
@@ -560,6 +786,8 @@ htsky_fetch_raw_property
/* Not in atmopshere => No gas */
if(!in_atmosphere) comp_mask &= ~HTSKY_CPNT_FLAG_GAS;
} else {
+ const double* upp;
+ const size_t* def;
world_to_cloud(sky, pos, pos_cs);
/* Compute the index of the voxel to fetch */
@@ -577,114 +805,153 @@ htsky_fetch_raw_property
((pos_cs[2] - cloud_desc->lower[2]) / sky->lut_cell_sz);
ivox[2] = splits[ilut].index + (pos_cs[2] > splits[ilut].height);
}
+
+ /* Handle numerical issues that may lead to a position lying onto the cloud
+ * upper boundaries */
+ def = sky->htcp_desc.spatial_definition;
+ upp = cloud_desc->upper;
+ if(ivox[0] == def[0] && eq_eps(pos_cs[0], upp[0], 1.e-6)) ivox[0] = def[0]-1;
+ if(ivox[1] == def[1] && eq_eps(pos_cs[1], upp[1], 1.e-6)) ivox[1] = def[1]-1;
+ if(ivox[2] == def[2] && eq_eps(pos_cs[2], upp[2], 1.e-6)) ivox[2] = def[2]-1;
+ if(ivox[0] >= def[0]) FATAL("Out of bound X voxel coordinate\n");
+ if(ivox[1] >= def[1]) FATAL("Out of bound Y voxel coordinate\n");
+ if(ivox[2] >= def[2]) FATAL("Out of bound Z voxel coordinate\n");
}
if(comp_mask & HTSKY_CPNT_FLAG_PARTICLES) {
- double rho_da = 0; /* Dry air density */
- double rct = 0; /* #droplets in kg of water per kg of dry air */
- double ql = 0; /* Droplet density In kg.m^-3 */
- double Ca = 0; /* Massic absorption cross section in m^2.kg^-1 */
- double Cs = 0; /* Massic scattering cross section in m^2.kg^-1 */
- ASSERT(in_clouds);
+ k_particle = particle_fetch_raw_property(sky, prop, iband, iquad, ivox);
+ }
+ if(comp_mask & HTSKY_CPNT_FLAG_GAS) {
+ k_gas = gas_fetch_raw_property
+ (sky, prop, iband, iquad, in_clouds, pos, ivox);
+ }
+ k = k_particle + k_gas;
+ ASSERT(k >= k_min && k <= k_max);
+ (void)k_min, (void)k_max;
+ return k;
+}
- /* Compute he dry air density */
- rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox);
+double
+htsky_fetch_temperature(const struct htsky* sky, const double pos[3])
+{
+ double temperature = 0;
+ struct htgop_level lvl0, lvl1;
+ size_t nlvls = 0;
+ int in_clouds = 0;
+ int in_atmosphere = 0;
+ ASSERT(sky && pos);
- /* Compute the droplet density */
- rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
- ql = rho_da * rct;
+ /* Is the position inside the atmosphere */
+ HTGOP(get_levels_count(sky->htgop, &nlvls));
+ ASSERT(nlvls > 1);
+ HTGOP(get_level(sky->htgop, 0, &lvl0));
+ HTGOP(get_level(sky->htgop, nlvls-1, &lvl1));
+ in_atmosphere = pos[2] >= lvl0.height && pos[2] <= lvl1.height;
- /* Use the average cross section of the current spectral band */
- if(prop == HTSKY_Ka || prop == HTSKY_Kext) Ca = sky->sw_bands[i].Ca_avg;
- if(prop == HTSKY_Ks || prop == HTSKY_Kext) Cs = sky->sw_bands[i].Cs_avg;
+ /* Is the position inside the clouds? */
+ if(!sky->is_cloudy) {
+ in_clouds = 0;
+ } else if(sky->repeat_clouds) {
+ in_clouds =
+ pos[2] >= sky->htcp_desc.lower[2]
+ && pos[2] <= sky->htcp_desc.upper[2];
+ } else {
+ in_clouds =
+ pos[0] >= sky->htcp_desc.lower[0]
+ && pos[1] >= sky->htcp_desc.lower[1]
+ && pos[2] >= sky->htcp_desc.lower[2]
+ && pos[0] <= sky->htcp_desc.upper[0]
+ && pos[1] <= sky->htcp_desc.upper[1]
+ && pos[2] <= sky->htcp_desc.upper[2];
+ }
- k_particle = ql*(Ca + Cs);
+ if(in_clouds) {
+ /* Clouds have priority on atmosphere */
+ in_atmosphere = 0;
}
- if(comp_mask & HTSKY_CPNT_FLAG_GAS) {
- struct htgop_layer layer;
- struct htgop_layer_sw_spectral_interval band;
- size_t ilayer = 0;
- ASSERT(in_atmosphere);
+ if(in_atmosphere) {
+ double u;
+ size_t ilayer;
- /* Retrieve the quadrature point into the spectral band of the layer into
- * which `pos' lies */
+ /* Find the layer into which pos is included */
HTGOP(position_to_layer_id(sky->htgop, pos[2], &ilayer));
- HTGOP(get_layer(sky->htgop, ilayer, &layer));
- HTGOP(layer_get_sw_spectral_interval(&layer, iband, &band));
+ ASSERT(ilayer < nlvls-1);
- if(!in_clouds) {
- /* Pos is outside the clouds. Directly fetch the nominal optical
- * properties */
- ASSERT(iquad < band.quadrature_length);
- switch(prop) {
- case HTSKY_Ka: k_gas = band.ka_nominal[iquad]; break;
- case HTSKY_Ks: k_gas = band.ks_nominal[iquad]; break;
- case HTSKY_Kext:
- k_gas = band.ka_nominal[iquad] + band.ks_nominal[iquad];
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- } else {
- /* Pos is inside the clouds. Compute the water vapor molar fraction at
- * the current voxel */
- const double x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
- struct htgop_layer_sw_spectral_interval_tab tab;
+ /* Fetch the levels enclosing the current layer */
+ HTGOP(get_level(sky->htgop, ilayer+0, &lvl0));
+ HTGOP(get_level(sky->htgop, ilayer+1, &lvl1));
+ ASSERT(lvl0.height < lvl1.height);
+ ASSERT(lvl0.height <= pos[2] && pos[2] <= lvl1.height);
- /* Retrieve the tabulated data for the quadrature point */
- HTGOP(layer_sw_spectral_interval_get_tab(&band, iquad, &tab));
+ /* Linearly interpolate the temperature of the levels into which pos lies */
+ u = (pos[2] - lvl0.height) / (lvl1.height - lvl0.height);
+ temperature = u * (lvl1.temperature - lvl0.temperature) + lvl0.temperature;
- /* Fetch the optical properties wrt x_h2o */
- switch(prop) {
- case HTSKY_Ka:
- HTGOP(layer_sw_spectral_interval_tab_fetch_ka(&tab, x_h2o, &k_gas));
- break;
- case HTSKY_Ks:
- HTGOP(layer_sw_spectral_interval_tab_fetch_ks(&tab, x_h2o, &k_gas));
- break;
- case HTSKY_Kext:
- HTGOP(layer_sw_spectral_interval_tab_fetch_kext(&tab, x_h2o, &k_gas));
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
+ } else if(in_clouds) {
+ const struct htcp_desc* desc = &sky->htcp_desc;
+ size_t ivox[3];
+ double pos_cs[3];
+
+ /* Transform the submitted position in local cloud space */
+ world_to_cloud(sky, pos, pos_cs);
+
+ /* Compute the index of the voxel to fetch */
+ ivox[0] = (size_t)((pos_cs[0] - desc->lower[0])/desc->vxsz_x);
+ ivox[1] = (size_t)((pos_cs[1] - desc->lower[1])/desc->vxsz_y);
+ if(!desc->irregular_z) {
+ /* The voxels along the Z dimension have the same size */
+ ivox[2] = (size_t)((pos_cs[2] - desc->lower[2])/desc->vxsz_z[0]);
+ } else {
+ /* Irregular voxel size along the Z dimension. Compute the index of the Z
+ * position in the svx2htcp_z Look Up Table and use the LUT to define the
+ * voxel index into the HTCP descripptor */
+ const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z);
+ const size_t ilut = (size_t)((pos_cs[2] - desc->lower[2]) / sky->lut_cell_sz);
+ ivox[2] = splits[ilut].index + (pos_cs[2] > splits[ilut].height);
}
+
+ /* Fetch the cloud temperature */
+ temperature = htcp_desc_T_at(desc, ivox[0], ivox[1], ivox[2], 0);
}
- k = k_particle + k_gas;
- ASSERT(k >= k_min && k <= k_max);
- (void)k_min, (void)k_max;
- return k;
+ return temperature;
}
size_t
-htsky_get_sw_spectral_bands_count(const struct htsky* sky)
+htsky_get_spectral_bands_count(const struct htsky* sky)
{
- ASSERT(sky && sky->sw_bands_range[0] <= sky->sw_bands_range[1]);
- return sky->sw_bands_range[1] - sky->sw_bands_range[0] + 1;
+ ASSERT(sky && sky->bands_range[0] <= sky->bands_range[1]);
+ return sky->bands_range[1] - sky->bands_range[0] + 1;
}
size_t
-htsky_get_sw_spectral_band_id
+htsky_get_spectral_band_id
(const struct htsky* sky, const size_t i)
{
- ASSERT(sky && i < htsky_get_sw_spectral_bands_count(sky));
- return sky->sw_bands_range[0] + i;
+ ASSERT(sky);
+ ASSERT(i < htsky_get_spectral_bands_count(sky));
+ return sky->bands_range[0] + i;
}
size_t
-htsky_get_sw_spectral_band_quadrature_length
+htsky_get_spectral_band_quadrature_length
(const struct htsky* sky, const size_t iband)
{
struct htgop_spectral_interval band;
ASSERT(sky);
- ASSERT(iband >= sky->sw_bands_range[0]);
- ASSERT(iband <= sky->sw_bands_range[1]);
- HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ ASSERT(iband >= sky->bands_range[0]);
+ ASSERT(iband <= sky->bands_range[1]);
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
return band.quadrature_length;
}
res_T
-htsky_get_sw_spectral_band_bounds
+htsky_get_spectral_band_bounds
(const struct htsky* sky,
const size_t iband,
double wavelengths[2])
@@ -693,52 +960,58 @@ htsky_get_sw_spectral_band_bounds
res_T res = RES_OK;
ASSERT(sky && wavelengths);
- res = htgop_get_sw_spectral_interval(sky->htgop, iband, &specint);
- if(res != RES_OK) return res;
-
+ if(sky->is_long_wave) {
+ res = htgop_get_lw_spectral_interval(sky->htgop, iband, &specint);
+ if(res != RES_OK) return res;
+ } else {
+ res = htgop_get_sw_spectral_interval(sky->htgop, iband, &specint);
+ if(res != RES_OK) return res;
+ }
wavelengths[0] = wavenumber_to_wavelength(specint.wave_numbers[1]);
wavelengths[1] = wavenumber_to_wavelength(specint.wave_numbers[0]);
ASSERT(wavelengths[0] < wavelengths[1]);
return RES_OK;
}
-void
-htsky_sample_sw_spectral_data_CIE_1931_X
- (const struct htsky* sky,
- const double r0, /* Random number in [0, 1[ */
- const double r1, /* Random number in [0, 1[ */
- size_t* ispectral_band,
- size_t* iquadrature_pt)
+int
+htsky_is_long_wave(const struct htsky* htsky)
{
- sample_sw_spectral_data
- (sky->htgop, r0, r1, htgop_sample_sw_spectral_interval_CIE_1931_X,
- ispectral_band, iquadrature_pt);
+ ASSERT(htsky);
+ return htsky->is_long_wave;
}
-void
-htsky_sample_sw_spectral_data_CIE_1931_Y
- (const struct htsky* sky,
- const double r0, /* Random number in [0, 1[ */
- const double r1, /* Random number in [0, 1[ */
- size_t* ispectral_band,
- size_t* iquadrature_pt)
+size_t
+htsky_find_spectral_band(const struct htsky* sky, const double wavelength)
{
- sample_sw_spectral_data
- (sky->htgop, r0, r1, htgop_sample_sw_spectral_interval_CIE_1931_Y,
- ispectral_band, iquadrature_pt);
+ const double wnum = wavelength_to_wavenumber(wavelength);
+ size_t iband;
+ ASSERT(sky);
+ if(sky->is_long_wave) {
+ HTGOP(find_lw_spectral_interval_id(sky->htgop, wnum, &iband));
+ } else {
+ HTGOP(find_sw_spectral_interval_id(sky->htgop, wnum, &iband));
+ }
+ return iband;
}
-void
-htsky_sample_sw_spectral_data_CIE_1931_Z
+size_t
+htsky_spectral_band_sample_quadrature
(const struct htsky* sky,
- const double r0, /* Random number in [0, 1[ */
- const double r1, /* Random number in [0, 1[ */
- size_t* ispectral_band,
- size_t* iquadrature_pt)
+ const double r,
+ const size_t iband)
{
- sample_sw_spectral_data
- (sky->htgop, r0, r1, htgop_sample_sw_spectral_interval_CIE_1931_Z,
- ispectral_band, iquadrature_pt);
+ struct htgop_spectral_interval band;
+ size_t iquad;
+ ASSERT(sky);
+ ASSERT(sky->bands_range[0] <= iband || iband <= sky->bands_range[1]);
+
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
+ HTGOP(spectral_interval_sample_quadrature(&band, r, &iquad));
+ return iquad;
}
/*******************************************************************************
@@ -751,6 +1024,7 @@ world_to_cloud
double out_pos_cs[3])
{
double cloud_sz[2];
+ double upper[2];
double pos_cs[3];
double pos_cs_n[2];
ASSERT(sky && pos_ws && out_pos_cs);
@@ -769,14 +1043,18 @@ world_to_cloud
return d3_set(out_pos_cs, pos_ws);
}
- cloud_sz[0] = sky->htcp_desc.upper[0] - sky->htcp_desc.lower[0];
- cloud_sz[1] = sky->htcp_desc.upper[1] - sky->htcp_desc.lower[1];
+ /* The cloud upper bound is not inclusive. Define the inclusive upper bound
+ * of the cloud */
+ upper[0] = nextafter(sky->htcp_desc.upper[0], sky->htcp_desc.lower[0]);
+ upper[1] = nextafter(sky->htcp_desc.upper[1], sky->htcp_desc.lower[1]);
+ cloud_sz[0] = upper[0] - sky->htcp_desc.lower[0];
+ cloud_sz[1] = upper[1] - sky->htcp_desc.lower[1];
/* Transform pos in normalize local cloud space */
pos_cs_n[0] = (pos_ws[0] - sky->htcp_desc.lower[0]) / cloud_sz[0];
pos_cs_n[1] = (pos_ws[1] - sky->htcp_desc.lower[1]) / cloud_sz[1];
- pos_cs_n[0] -= (int)pos_cs_n[0]; /* Keep fractional part */
- pos_cs_n[1] -= (int)pos_cs_n[1]; /* Keep fractional part */
+ pos_cs_n[0] -= (int)pos_cs_n[0]; /* Get fractional part */
+ pos_cs_n[1] -= (int)pos_cs_n[1]; /* Get fractional part */
if(pos_cs_n[0] < 0) pos_cs_n[0] += 1;
if(pos_cs_n[1] < 0) pos_cs_n[1] += 1;
diff --git a/src/htsky.h b/src/htsky.h
@@ -69,6 +69,7 @@ struct htsky_args {
const char* htmie_filename;
const char* cache_filename; /* May be NULL <=> no cached data structure */
const char* name; /* Name of the sky. Used by the Star-MTL binding */
+ double wlen_lw_range[2]; /* Long wave wavelength range to handle. In nm */
unsigned grid_max_definition[3]; /* Maximum definition of the grid */
double optical_thickness; /* Threshold used during octree building */
unsigned nthreads; /* Hint on the number of threads to use */
@@ -82,6 +83,7 @@ struct htsky_args {
NULL, /* htmie filename */ \
NULL, /* cache filename */ \
"sky", /* Name */ \
+ {DBL_MAX,-DBL_MAX}, /* Long wave wavelengths. Degenerated <=> short wave */ \
{UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \
1, /* Optical thickness a*/ \
(unsigned)~0, /* #threads */ \
@@ -130,6 +132,11 @@ htsky_fetch_particle_phase_function_asymmetry_parameter
const size_t iquad);
HTSKY_API double
+htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter
+ (const struct htsky* sky,
+ const double wavelength); /* In nanometer */
+
+HTSKY_API double
htsky_fetch_raw_property
(const struct htsky* sky,
const enum htsky_property prop,
@@ -142,6 +149,11 @@ htsky_fetch_raw_property
const double k_max);
HTSKY_API double
+htsky_fetch_temperature
+ (const struct htsky* sky,
+ const double pos[3]);
+
+HTSKY_API double
htsky_fetch_svx_property
(const struct htsky* sky,
const enum htsky_property prop,
@@ -175,48 +187,42 @@ htsky_trace_ray
struct svx_hit* hit);
HTSKY_API size_t
-htsky_get_sw_spectral_bands_count
+htsky_get_spectral_bands_count
(const struct htsky* sky);
HTSKY_API size_t
-htsky_get_sw_spectral_band_id
+htsky_get_spectral_band_id
(const struct htsky* sky,
const size_t i);
HTSKY_API size_t
-htsky_get_sw_spectral_band_quadrature_length
+htsky_get_spectral_band_quadrature_length
(const struct htsky* sky,
const size_t iband);
HTSKY_API res_T
-htsky_get_sw_spectral_band_bounds
+htsky_get_spectral_band_bounds
(const struct htsky* sky,
const size_t iband,
double wavelengths[2]);
-HTSKY_API void
-htsky_sample_sw_spectral_data_CIE_1931_X
- (const struct htsky* sky,
- const double r0, /* Random number in [0, 1[ */
- const double r1, /* Random number in [0, 1[ */
- size_t* ispectral_band,
- size_t* iquadrature_pt);
+HTSKY_API int
+htsky_is_long_wave
+ (const struct htsky* sky);
-HTSKY_API void
-htsky_sample_sw_spectral_data_CIE_1931_Y
+/* Return the index of the band containing the submitted wavelength or SIZE_MAX
+ * if `wavelength' is not included in any band */
+HTSKY_API size_t
+htsky_find_spectral_band
(const struct htsky* sky,
- const double r0, /* Random number in [0, 1[ */
- const double r1, /* Random number in [0, 1[ */
- size_t* ispectral_band,
- size_t* iquadrature_pt);
+ const double wavelength); /* In nanometer */
-HTSKY_API void
-htsky_sample_sw_spectral_data_CIE_1931_Z
+/* Return the sampled quadrature point */
+HTSKY_API size_t
+htsky_spectral_band_sample_quadrature
(const struct htsky* sky,
- const double r0, /* Random number in [0, 1[ */
- const double r1, /* Random number in [0, 1[ */
- size_t* ispectral_band,
- size_t* iquadrature_pt);
+ const double r, /* Random number in [0, 1[ */
+ const size_t ispectral_band);
HTSKY_API res_T
htsky_dump_cloud_vtk
diff --git a/src/htsky_atmosphere.c b/src/htsky_atmosphere.c
@@ -58,27 +58,53 @@ atmosphere_vox_get(const size_t xyz[3], void* dst, void* context)
ka_max = ks_max = kext_max =-DBL_MAX;
/* For each atmospheric layer that overlaps the SVX voxel ... */
- FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) {
- struct htgop_layer layer;
- struct htgop_layer_sw_spectral_interval band;
- size_t iquad;
-
- HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer));
-
- /* ... retrieve the considered spectral interval */
- HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band));
-
- /* ... and update the radiative properties bound with the per quadrature
- * point nominal values */
- ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
- ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
- FOR_EACH(iquad, ctx->quadrature_range[0], ctx->quadrature_range[1]+1) {
- ka_min = MMIN(ka_min, band.ka_nominal[iquad]);
- ka_max = MMAX(ka_max, band.ka_nominal[iquad]);
- ks_min = MMIN(ks_min, band.ks_nominal[iquad]);
- ks_max = MMAX(ks_max, band.ks_nominal[iquad]);
- kext_min = MMIN(kext_min, band.ka_nominal[iquad]+band.ks_nominal[iquad]);
- kext_max = MMAX(kext_max, band.ka_nominal[iquad]+band.ks_nominal[iquad]);
+ if(ctx->sky->is_long_wave) {
+ FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) {
+ struct htgop_layer layer;
+ struct htgop_layer_lw_spectral_interval band;
+ size_t iquad;
+
+ HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer));
+
+ /* ... retrieve the considered spectral interval */
+ HTGOP(layer_get_lw_spectral_interval(&layer, ctx->iband, &band));
+
+ /* ... and update the radiative properties bound with the per quadrature
+ * point nominal values */
+ ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
+ ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
+ FOR_EACH(iquad, ctx->quadrature_range[0], ctx->quadrature_range[1]+1) {
+ ka_min = MMIN(ka_min, band.ka_nominal[iquad]);
+ ka_max = MMAX(ka_max, band.ka_nominal[iquad]);
+ }
+ }
+ ks_min = 0;
+ ks_max = 0;
+ kext_min = ka_min;
+ kext_max = ka_max;
+ } else {
+ FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) {
+ struct htgop_layer layer;
+ struct htgop_layer_sw_spectral_interval band;
+ size_t iquad;
+
+ HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer));
+
+ /* ... retrieve the considered spectral interval */
+ HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band));
+
+ /* ... and update the radiative properties bound with the per quadrature
+ * point nominal values */
+ ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
+ ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
+ FOR_EACH(iquad, ctx->quadrature_range[0], ctx->quadrature_range[1]+1) {
+ ka_min = MMIN(ka_min, band.ka_nominal[iquad]);
+ ka_max = MMAX(ka_max, band.ka_nominal[iquad]);
+ ks_min = MMIN(ks_min, band.ks_nominal[iquad]);
+ ks_max = MMAX(ks_max, band.ks_nominal[iquad]);
+ kext_min = MMIN(kext_min, band.ka_nominal[iquad]+band.ks_nominal[iquad]);
+ kext_max = MMAX(kext_max, band.ka_nominal[iquad]+band.ks_nominal[iquad]);
+ }
}
}
@@ -170,12 +196,12 @@ atmosphere_setup(struct htsky* sky, const double optical_thickness_threshold)
/* Create as many atmospheric data structure than considered SW spectral
* bands */
- nbands = htsky_get_sw_spectral_bands_count(sky);
+ nbands = htsky_get_spectral_bands_count(sky);
sky->atmosphere = MEM_CALLOC
(sky->allocator, nbands, sizeof(*sky->atmosphere));
if(!sky->atmosphere) {
log_err(sky,
- "could not create the list of per SW band atmospheric data structure.\n");
+ "Could not create the list of per SW band atmospheric data structure.\n");
res = RES_MEM_ERR;
goto error;
}
@@ -183,15 +209,19 @@ atmosphere_setup(struct htsky* sky, const double optical_thickness_threshold)
FOR_EACH(i, 0, nbands) {
size_t iquad;
struct htgop_spectral_interval band;
- ctx.iband = i + sky->sw_bands_range[0];
+ ctx.iband = i + sky->bands_range[0];
- HTGOP(get_sw_spectral_interval(sky->htgop, ctx.iband, &band));
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, ctx.iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, ctx.iband, &band));
+ }
sky->atmosphere[i] = MEM_CALLOC(sky->allocator,
band.quadrature_length, sizeof(*sky->atmosphere[i]));
if(!sky->atmosphere[i]) {
log_err(sky,
- "could not create the list of per quadrature point atmospheric data "
+ "Could not create the list of per quadrature point atmospheric data "
"for the band %lu.\n", (unsigned long)ctx.iband);
res = RES_MEM_ERR;
goto error;
@@ -207,7 +237,7 @@ atmosphere_setup(struct htsky* sky, const double optical_thickness_threshold)
res = svx_bintree_create(sky->svx, low, upp, definition, SVX_AXIS_Z,
&vox_desc, &sky->atmosphere[i][iquad].bitree);
if(res != RES_OK) {
- log_err(sky, "could not create the binary tree of the "
+ log_err(sky, "Could not create the binary tree of the "
"atmospheric properties for the band %lu.\n", (unsigned long)ctx.iband);
goto error;
}
@@ -234,17 +264,21 @@ atmosphere_clean(struct htsky* sky)
if(!sky->atmosphere) return;
- nbands = htsky_get_sw_spectral_bands_count(sky);
+ nbands = htsky_get_spectral_bands_count(sky);
FOR_EACH(i, 0, nbands) {
struct htgop_spectral_interval band;
size_t iband;
size_t iquad;
- iband = sky->sw_bands_range[0] + i;
- HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
-
if(!sky->atmosphere[i]) continue;
+ iband = sky->bands_range[0] + i;
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
+
FOR_EACH(iquad, 0, band.quadrature_length) {
if(sky->atmosphere[i][iquad].bitree) {
SVX(tree_ref_put(sky->atmosphere[i][iquad].bitree));
diff --git a/src/htsky_c.h b/src/htsky_c.h
@@ -49,7 +49,7 @@ struct split {
#include <rsys/dynamic_array.h>
/* Properties of a short wave spectral band */
-struct sw_band_prop {
+struct band_prop {
/* Average cross section in the band */
double Ca_avg; /* Absorption cross section */
double Cs_avg; /* Scattering cross section */
@@ -79,10 +79,11 @@ struct htsky {
struct darray_split svx2htcp_z;
double lut_cell_sz; /* Size of a svx2htcp_z cell */
- /* Ids and optical properties of the short wave spectral bands loaded by
- * HTGOP and that overlap the CIE XYZ color space */
- size_t sw_bands_range[2];
- struct sw_band_prop* sw_bands;
+ /* Ids and optical properties of spectral bands loaded by HTGOP and currently
+ * handled by the sky */
+ int is_long_wave;
+ size_t bands_range[2]; /* Inclusive band ids */
+ struct band_prop* bands;
int repeat_clouds; /* Make clouds infinite in X and Y */
int is_cloudy; /* The sky has clouds */
diff --git a/src/htsky_cloud.c b/src/htsky_cloud.c
@@ -81,7 +81,7 @@ setup_svx2htcp_z_lut(struct htsky* sky, double* cloud_upp_z)
res = darray_split_resize(&sky->svx2htcp_z, nsplits);
if(res != RES_OK) {
log_err(sky,
- "could not allocate the table mapping regular to irregular Z.\n");
+ "Could not allocate the table mapping regular to irregular Z.\n");
goto error;
}
@@ -134,11 +134,11 @@ compute_k_bounds_regular_z
ASSERT(!sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp);
ASSERT(ka_bounds && ks_bounds && kext_bounds);
- i = iband - sky->sw_bands_range[0];
+ i = iband - sky->bands_range[0];
/* Fetch the optical properties of the spectral band */
- Ca_avg = sky->sw_bands[i].Ca_avg;
- Cs_avg = sky->sw_bands[i].Cs_avg;
+ Ca_avg = sky->bands[i].Ca_avg;
+ Cs_avg = sky->bands[i].Cs_avg;
/* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
* the SVX voxel */
@@ -219,11 +219,11 @@ compute_k_bounds_irregular_z
ASSERT(sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp);
ASSERT(ka_bounds && ks_bounds && kext_bounds);
- i = iband - sky->sw_bands_range[0];
+ i = iband - sky->bands_range[0];
/* Fetch the optical properties of the spectral band */
- Ca_avg = sky->sw_bands[i].Ca_avg;
- Cs_avg = sky->sw_bands[i].Cs_avg;
+ Ca_avg = sky->bands[i].Ca_avg;
+ Cs_avg = sky->bands[i].Cs_avg;
/* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
* the SVX voxel along the X and Y axes */
@@ -490,7 +490,6 @@ cloud_vox_get_gas
{
const struct htcp_desc* htcp_desc;
struct htgop_layer layer;
- struct htgop_layer_sw_spectral_interval band;
size_t ilayer;
size_t layer_range[2];
double x_h2o_range[2];
@@ -524,29 +523,54 @@ cloud_vox_get_gas
ka[1] = ks[1] = kext[1] =-DBL_MAX;
/* For each atmospheric layer that overlaps the SVX voxel ... */
- FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) {
- double k[2];
-
- HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer));
-
- /* ... retrieve the considered spectral interval */
- HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band));
- ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
- ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
-
- /* ... and compute the radiative properties and upd their bounds */
- HTGOP(layer_sw_spectral_interval_quadpoints_get_ka_bounds
- (&band, ctx->quadrature_range, x_h2o_range, k));
- ka[0] = MMIN(ka[0], k[0]);
- ka[1] = MMAX(ka[1], k[1]);
- HTGOP(layer_sw_spectral_interval_quadpoints_get_ks_bounds
- (&band, ctx->quadrature_range, x_h2o_range, k));
- ks[0] = MMIN(ks[0], k[0]);
- ks[1] = MMAX(ks[1], k[1]);
- HTGOP(layer_sw_spectral_interval_quadpoints_get_kext_bounds
- (&band, ctx->quadrature_range, x_h2o_range, k));
- kext[0] = MMIN(kext[0], k[0]);
- kext[1] = MMAX(kext[1], k[1]);
+ if(ctx->sky->is_long_wave) {
+ FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) {
+ struct htgop_layer_lw_spectral_interval band;
+ double k[2];
+
+ HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer));
+
+ /* ... retrieve the considered spectral interval */
+ HTGOP(layer_get_lw_spectral_interval(&layer, ctx->iband, &band));
+ ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
+ ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
+
+ /* ... and compute the radiative properties and upd their bounds */
+ HTGOP(layer_lw_spectral_interval_quadpoints_get_ka_bounds
+ (&band, ctx->quadrature_range, x_h2o_range, k));
+ ka[0] = MMIN(ka[0], k[0]);
+ ka[1] = MMAX(ka[1], k[1]);
+ }
+ ks[0] = 0;
+ ks[1] = 0;
+ kext[0] = ka[0];
+ kext[1] = ka[1];
+ } else {
+ FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) {
+ struct htgop_layer_sw_spectral_interval band;
+ double k[2];
+
+ HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer));
+
+ /* ... retrieve the considered spectral interval */
+ HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band));
+ ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
+ ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
+
+ /* ... and compute the radiative properties and upd their bounds */
+ HTGOP(layer_sw_spectral_interval_quadpoints_get_ka_bounds
+ (&band, ctx->quadrature_range, x_h2o_range, k));
+ ka[0] = MMIN(ka[0], k[0]);
+ ka[1] = MMAX(ka[1], k[1]);
+ HTGOP(layer_sw_spectral_interval_quadpoints_get_ks_bounds
+ (&band, ctx->quadrature_range, x_h2o_range, k));
+ ks[0] = MMIN(ks[0], k[0]);
+ ks[1] = MMAX(ks[1], k[1]);
+ HTGOP(layer_sw_spectral_interval_quadpoints_get_kext_bounds
+ (&band, ctx->quadrature_range, x_h2o_range, k));
+ kext[0] = MMIN(kext[0], k[0]);
+ kext[1] = MMAX(kext[1], k[1]);
+ }
}
/* Ensure that the single precision bounds include their double precision
@@ -646,7 +670,7 @@ cloud_build_octrees
struct build_tree_context ctx = BUILD_TREE_CONTEXT_NULL;
const size_t iband = darray_specdata_data_get(specdata)[ispecdata].iband;
const size_t iquad = darray_specdata_data_get(specdata)[ispecdata].iquad;
- const size_t id = iband - sky->sw_bands_range[0];
+ const size_t id = iband - sky->bands_range[0];
int32_t pcent;
size_t n;
res_T res_local = RES_OK;
@@ -723,14 +747,18 @@ cloud_write_octrees(const struct htsky* sky, FILE* stream)
res_T res = RES_OK;
ASSERT(sky && sky->clouds);
- nbands = htsky_get_sw_spectral_bands_count(sky);
+ nbands = htsky_get_spectral_bands_count(sky);
FOR_EACH(i, 0, nbands) {
struct htgop_spectral_interval band;
size_t iband;
size_t iquad;
- iband = sky->sw_bands_range[0] + i;
- HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ iband = sky->bands_range[0] + i;
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
ASSERT(sky->clouds[i]);
FOR_EACH(iquad, 0, band.quadrature_length) {
@@ -758,21 +786,25 @@ cloud_read_octrees(struct htsky* sky, const char* stream_name, FILE* stream)
ASSERT(sky && sky->clouds && stream_name);
log_info(sky, "Loading octrees from `%s'.\n", stream_name);
- nbands = htsky_get_sw_spectral_bands_count(sky);
+ nbands = htsky_get_spectral_bands_count(sky);
FOR_EACH(i, 0, nbands) {
struct htgop_spectral_interval band;
size_t iband;
size_t iquad;
- iband = sky->sw_bands_range[0] + i;
- HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ iband = sky->bands_range[0] + i;
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
ASSERT(sky->clouds[i]);
FOR_EACH(iquad, 0, band.quadrature_length) {
/* The octree was already setup */
if(sky->clouds[i][iquad].octree != NULL) continue;
- res = svx_tree_create_from_stream
+ res = svx_tree_create_from_stream
(sky->svx, stream, &sky->clouds[i][iquad].octree);
if(res != RES_OK) {
log_err(sky, "Could not read the octree %lu:%lu from `%s' -- %s.\n",
@@ -812,7 +844,7 @@ cloud_setup
size_t nbands;
size_t i;
ATOMIC res = RES_OK;
- ASSERT(sky && sky->sw_bands && optical_thickness_threshold >= 0);
+ ASSERT(sky && sky->bands && optical_thickness_threshold >= 0);
ASSERT(grid_max_definition);
ASSERT(grid_max_definition[0]);
ASSERT(grid_max_definition[1]);
@@ -822,7 +854,7 @@ cloud_setup
res = htcp_get_desc(sky->htcp, &sky->htcp_desc);
if(res != RES_OK) {
- log_err(sky, "could not retrieve the HTCP descriptor.\n");
+ log_err(sky, "Could not retrieve the HTCP descriptor.\n");
goto error;
}
@@ -849,11 +881,11 @@ cloud_setup
}
/* Create as many cloud data structure than considered SW spectral bands */
- nbands = htsky_get_sw_spectral_bands_count(sky);
+ nbands = htsky_get_spectral_bands_count(sky);
sky->clouds = MEM_CALLOC(sky->allocator, nbands, sizeof(*sky->clouds));
if(!sky->clouds) {
log_err(sky,
- "could not create the list of per SW band cloud data structure.\n");
+ "Could not create the list of per SW band cloud data structure.\n");
res = RES_MEM_ERR;
goto error;
}
@@ -861,10 +893,14 @@ cloud_setup
/* Compute how many octrees are going to be built */
FOR_EACH(i, 0, nbands) {
struct htgop_spectral_interval band;
- const size_t iband = i + sky->sw_bands_range[0];
+ const size_t iband = i + sky->bands_range[0];
size_t iquad;
- HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
sky->clouds[i] = MEM_CALLOC(sky->allocator,
band.quadrature_length, sizeof(*sky->clouds[i]));
@@ -947,14 +983,18 @@ cloud_clean(struct htsky* sky)
if(!sky->clouds) return;
- nbands = htsky_get_sw_spectral_bands_count(sky);
+ nbands = htsky_get_spectral_bands_count(sky);
FOR_EACH(i, 0, nbands) {
struct htgop_spectral_interval band;
size_t iband;
size_t iquad;
- iband = sky->sw_bands_range[0] + i;
- HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ iband = sky->bands_range[0] + i;
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ }
if(!sky->clouds[i]) continue;
diff --git a/src/htsky_dump_cloud_vtk.c b/src/htsky_dump_cloud_vtk.c
@@ -70,8 +70,8 @@ octree_data_init
struct octree_data* data)
{
ASSERT(data);
- ASSERT(iband >= sky->sw_bands_range[0]);
- ASSERT(iband <= sky->sw_bands_range[1]);
+ ASSERT(iband >= sky->bands_range[0]);
+ ASSERT(iband <= sky->bands_range[1]);
(void)iquad;
htable_vertex_init(sky->allocator, &data->vertex2id);
darray_double_init(sky->allocator, &data->vertices);
@@ -165,19 +165,19 @@ htsky_dump_cloud_vtk
goto error;
}
- if(iband >= sky->sw_bands_range[0] && iband <= sky->sw_bands_range[1]) {
+ if(iband < sky->bands_range[0] || iband > sky->bands_range[1]) {
log_err(sky,
"%s: invalid spectral band index `%lu'. "
"Valid short wave spectral bands lie in [%lu, %lu]\n",
FUNC_NAME,
(unsigned long)iband,
- (unsigned long)sky->sw_bands_range[0],
- (unsigned long)sky->sw_bands_range[1]);
+ (unsigned long)sky->bands_range[0],
+ (unsigned long)sky->bands_range[1]);
res = RES_BAD_ARG;
goto error;
}
- if(iquad >= htsky_get_sw_spectral_band_quadrature_length(sky, iband)) {
+ if(iquad >= htsky_get_spectral_band_quadrature_length(sky, iband)) {
log_err(sky,
"invalid quadrature point `%lu' for the spectral band `%lu'.\n",
(unsigned long)iquad, (unsigned long)iband);
@@ -190,7 +190,7 @@ htsky_dump_cloud_vtk
return RES_OK;
}
- i = iband - sky->sw_bands_range[0];
+ i = iband - sky->bands_range[0];
octree_data_init(sky, iband, iquad, &data);
cloud = &sky->clouds[i][iquad];
@@ -204,7 +204,11 @@ htsky_dump_cloud_vtk
ASSERT(ncells == cloud->octree_desc.nleaves);
/* Fetch the spectral interval descriptor */
- HTGOP(get_sw_spectral_interval(sky->htgop, iband, &specint));
+ if(sky->is_long_wave) {
+ HTGOP(get_lw_spectral_interval(sky->htgop, iband, &specint));
+ } else {
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &specint));
+ }
/* Write headers */
fprintf(stream, "# vtk DataFile Version 2.0\n");
diff --git a/src/htsky_smtl.c b/src/htsky_smtl.c
@@ -1,304 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * 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 /* getopt support */
-
-#include "htsky.h"
-#include "htsky_smtl.h"
-
-#include <rsys/cstr.h>
-
-#include <getopt.h>
-#include <string.h>
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-parse_grid_definition
- (const char* cmd,
- struct htsky_args* args,
- const char* str)
-{
- unsigned def[3];
- size_t len;
- res_T res = RES_OK;
- ASSERT(cmd && args && str);
-
- res = cstr_to_list_uint(str, ',', def, &len, 3);
- if(res == RES_OK && len != 3) res = RES_BAD_ARG;
- if(res != RES_OK) {
- fprintf(stderr, "%s: invalid grid definition `%s'.\n", cmd, str);
- goto error;
- }
-
- if(!def[0] || !def[1] || !def[2]) {
- fprintf(stderr,
- "%s: invalid null grid definition {%u, %u, %u}.\n", cmd, SPLIT3(def));
- res = RES_BAD_ARG;
- goto error;
- }
-
- args->grid_max_definition[0] = def[0];
- args->grid_max_definition[1] = def[1];
- args->grid_max_definition[2] = def[2];
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static void
-args_release(struct htsky_args* args)
-{
- ASSERT(args);
- *args = HTSKY_ARGS_DEFAULT;
-}
-
-static res_T
-args_init(struct htsky_args* args, int argc, char** argv)
-{
- int opt;
- res_T res = RES_OK;
- ASSERT(args);
-
- *args = HTSKY_ARGS_DEFAULT;
-
- /* Reset the getopt global variables */
- optind = 0;
- opterr = 0;
-
- while((opt = getopt(argc, argv, "a:c:m:n:O:rT:t:V:v")) != -1) {
- switch(opt) {
- case 'a': args->htgop_filename = optarg; break;
- case 'c': args->htcp_filename = optarg; break;
- case 'm': args->htmie_filename = optarg; break;
- case 'n': args->name = optarg; break;
- case 'O': args->cache_filename = optarg; break;
- case 'r': args->repeat_clouds = 1; break;
- case 'T':
- res = cstr_to_double(optarg, &args->optical_thickness);
- if(res == RES_OK && args->optical_thickness < 0) res = RES_BAD_ARG;
- break;
- case 't':
- res = cstr_to_uint(optarg, &args->nthreads);
- if(res == RES_OK && !args->nthreads) res = RES_BAD_ARG;
- break;
- case 'V': res = parse_grid_definition(argv[0], args, optarg); break;
- case 'v': args->verbose = 1; break;
- default: res = RES_BAD_ARG; break;
- }
- if(res != RES_OK) {
- if(optarg) {
- fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
- argv[0], optarg, opt);
- }
- goto error;
- }
- }
- if(!args->htgop_filename) {
- fprintf(stderr,
- "%s: missing the path of the gas optical properties file -- option '-a'\n",
- argv[0]);
- res = RES_BAD_ARG;
- goto error;
- }
- if(args->htcp_filename && !args->htmie_filename) {
- fprintf(stderr,
- "%s: missing the path toward the file of the Mie's data -- option '-m'\n",
- argv[0]);
- res = RES_BAD_ARG;
- goto error;
- }
-
-exit:
- return res;
-error:
- args_release(args);
- goto exit;
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-smtl_program_init
- (struct logger* logger, /* NULL <=> use default logger */
- struct mem_allocator* allocator, /* NULL <=> use default allocator */
- int argc,
- char* argv[],
- void** out_prog)
-{
- struct htsky* sky = NULL;
- struct htsky_args args = HTSKY_ARGS_DEFAULT;
- res_T res = RES_OK;
-
- if(!argc || !argv || !out_prog) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- res = args_init(&args, argc, argv);
- if(res != RES_OK) goto error;
-
- res = htsky_create(logger, allocator, &args, &sky);
- if(res != RES_OK) goto error;
-
-exit:
- if(out_prog) *out_prog = sky;
- return res;
-error:
- if(sky) {
- HTSKY(ref_put(sky));
- sky = NULL;
- }
- goto exit;
-}
-
-void
-smtl_program_release(void* program)
-{
- HTSKY(ref_put(program));
-}
-
-const char*
-smtl_program_get_mtl_name(void* program)
-{
- return htsky_get_name(program);
-}
-
-enum smtl_mtl_type
-smtl_prgram_get_mtl_type(void* program)
-{
- (void)program;
- return SMTL_MTL_SEMI_TRANSPARENT;
-}
-
-size_t
-smtl_program_get_cpnts_count(void* program)
-{
- (void)program;
- return HTSKY_CPNTS_COUNT__;
-}
-
-int
-smtl_program_find_cpnt(void* program, const char* cpnt_name)
-{
- int i;
- FOR_EACH(i, 0, HTSKY_CPNTS_COUNT__) {
- if(!strcmp(cpnt_name, smtl_program_cpnt_get_name(program, i)))
- return i;
- }
- return -1;
-}
-
-const char*
-smtl_program_cpnt_get_name(void* program, const int cpnt_id)
-{
- const char* name = NULL;
- (void)program;
- switch(cpnt_id) {
- case HTSKY_CPNT_GAS: name = "gas"; break;
- case HTSKY_CPNT_PARTICLES: name = "particles"; break;
- default: FATAL("Unreachable code.\n"); break;
- }
- return name;
-}
-
-double
-smtl_program_cpnt_get_ka
- (void* program,
- const int cpnt_id,
- const struct smtl_spectral_data* spec,
- const struct smtl_vertex* vtx)
-{
- struct htsky* sky = program;
- ASSERT(program && (unsigned)cpnt_id < HTSKY_CPNTS_COUNT__);
- ASSERT(spec && spec->type == SMTL_SPECTRAL_KDISTRIB);
- return htsky_fetch_raw_property(sky, HTSKY_Ka, BIT(cpnt_id),
- spec->value.kdistrib.iband, spec->value.kdistrib.iquad, vtx->pos,
- -DBL_MAX, DBL_MAX);
-}
-
-enum smtl_spectral_type
-smtl_program_cpnt_get_ka_spectral_type(void* program, const int cpnt_id)
-{
- (void)program, (void)cpnt_id;
- return SMTL_SPECTRAL_KDISTRIB;
-}
-
-double
-smtl_program_cpnt_get_ks
- (void* program,
- const int cpnt_id,
- const struct smtl_spectral_data* spec,
- const struct smtl_vertex* vtx)
-{
- struct htsky* sky = program;
- ASSERT(program && (unsigned)cpnt_id < HTSKY_CPNTS_COUNT__);
- ASSERT(spec && spec->type == SMTL_SPECTRAL_KDISTRIB);
- return htsky_fetch_raw_property(sky, HTSKY_Ks, BIT(cpnt_id),
- spec->value.kdistrib.iband, spec->value.kdistrib.iquad, vtx->pos,
- -DBL_MAX, DBL_MAX);
-}
-
-enum smtl_spectral_type
-smtl_program_cpnt_get_ks_spectral_type(void* program, const int cpnt_id)
-{
- (void)program, (void)cpnt_id;
- return SMTL_SPECTRAL_KDISTRIB;
-}
-
-enum smtl_phasefn_type
-smtl_program_cpnt_get_phasefn_type(void* program, const int cpnt_id)
-{
- enum smtl_phasefn_type phasefn = SMTL_PHASEFN_NONE__;
- (void)program;
- switch(cpnt_id) {
- case HTSKY_CPNT_GAS:
- phasefn = SMTL_PHASEFN_RAYLEIGH;
- break;
- case HTSKY_CPNT_PARTICLES:
- phasefn = SMTL_PHASEFN_HG;
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- return phasefn;
-}
-
-double
-smtl_program_cpnt_phasefn_hg_get_g
- (void* program,
- const int cpnt_id,
- const struct smtl_spectral_data* spec,
- const struct smtl_vertex* vtx)
-{
- struct htsky* sky = program;
- (void)vtx, (void)cpnt_id;
- ASSERT(program && cpnt_id == HTSKY_CPNT_PARTICLES);
- ASSERT(spec && spec->type == SMTL_SPECTRAL_KDISTRIB);
- return htsky_fetch_particle_phase_function_asymmetry_parameter
- (sky, spec->value.kdistrib.iband, spec->value.kdistrib.iquad);
-}
-
-enum smtl_spectral_type
-smtl_program_cpnt_phasefn_hg_get_g_spectral_type(void* program, const int cpnt_id)
-{
- (void)program, (void)cpnt_id;
- return SMTL_SPECTRAL_KDISTRIB;
-}
-
diff --git a/src/htsky_smtl.h b/src/htsky_smtl.h
@@ -1,129 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTSKY_SMTL_H
-#define HTSKY_SMTL_H
-
-#include <star/smtl.h>
-#include <rsys/rsys.h>
-
-/*
- * Usage: htsky [OPTION]... -a ATMOSPHERE
- * Manage the data representing a clear/cloudy sky.
- *
- * -a ATMOSPHERE HTGOP file of the gas optical properties of the atmosphere
- * -c CLOUDS HTCP file storing the properties of the clouds.
- * -m MIE HTMie file storing the optical properties of the clouds.
- * -n NAME name of the sky. BY default it is "sky"
- * -O CACHE name the cache file used to store/restore the sky data
- * -r infinitely repeat the clouds along the X and Y axis.
- * -T THRESHOLD optical thickness used as threshold during the octree\n"
- * building. By default its value is 1.
- * -t THREADS hint on the number of threads to use. By default use as
- * -V X,Y,Z maximum definition of the cloud acceleration grids along
- * the 3 axis. By default use the definition of the raw data.
- * -v make htsky verbose.
- */
-
-BEGIN_DECLS
-
-/*******************************************************************************
- * Common functions
- ******************************************************************************/
-HTSKY_API res_T
-smtl_program_init
- (struct logger* logger, /* NULL <=> use default logger */
- struct mem_allocator* allocator, /* NULL <=> use default allocator */
- int argc,
- char* argv[],
- void** out_prog);
-
-HTSKY_API void
-smtl_program_release
- (void* program);
-
-/*******************************************************************************
- * General material attribs
- ******************************************************************************/
-HTSKY_API const char*
-smtl_program_get_mtl_name
- (void* program);
-
-HTSKY_API enum smtl_mtl_type
-smtl_prgram_get_mtl_type
- (void* program);
-
-/*******************************************************************************
- * Component attributes
- ******************************************************************************/
-HTSKY_API size_t
-smtl_program_get_cpnts_count
- (void* program);
-
-HTSKY_API int
-smtl_program_find_cpnt
- (void* program,
- const char* cpnt_name);
-
-HTSKY_API const char*
-smtl_program_cpnt_get_name
- (void* program,
- const int cpnt_id);
-
-HTSKY_API double
-smtl_program_cpnt_get_ka
- (void* program,
- const int cpnt_id,
- const struct smtl_spectral_data* spec,
- const struct smtl_vertex* vtx);
-
-HTSKY_API enum smtl_spectral_type
-smtl_program_cpnt_get_ka_spectral_type
- (void* program,
- const int cpnt_id);
-
-HTSKY_API double
-smtl_program_cpnt_get_ks
- (void* program,
- const int cpnt_id,
- const struct smtl_spectral_data* spec,
- const struct smtl_vertex* vtx);
-
-HTSKY_API enum smtl_spectral_type
-smtl_program_cpnt_get_ks_spectral_type
- (void* program,
- const int cpnt_id);
-
-HTSKY_API double
-smtl_program_cpnt_phasefn_hg_get_g
- (void* program,
- const int cpnt_id,
- const struct smtl_spectral_data* spec,
- const struct smtl_vertex* vtx);
-
-HTSKY_API enum smtl_spectral_type
-smtl_program_cpnt_phasefn_hg_get_g_spectral_type
- (void* program,
- const int cpnt_id);
-
-HTSKY_API enum smtl_phasefn_type
-smtl_program_cpnt_get_phasefn_type
- (void* program,
- const int cpnt_id);
-
-END_DECLS
-
-#endif /* HTSKY_SMTL_H */
diff --git a/src/htsky_svx.c b/src/htsky_svx.c
@@ -205,11 +205,11 @@ htsky_fetch_svx_property
int comp_mask = components_mask;
ASSERT(sky && pos);
ASSERT(comp_mask & HTSKY_CPNT_MASK_ALL);
- ASSERT(iband >= sky->sw_bands_range[0]);
- ASSERT(iband <= sky->sw_bands_range[1]);
- ASSERT(iquad < htsky_get_sw_spectral_band_quadrature_length(sky, iband));
+ ASSERT(iband >= sky->bands_range[0]);
+ ASSERT(iband <= sky->bands_range[1]);
+ ASSERT(iquad < htsky_get_spectral_band_quadrature_length(sky, iband));
- i = iband - sky->sw_bands_range[0];
+ i = iband - sky->bands_range[0];
cloud = sky->is_cloudy ? &sky->clouds[i][iquad] : NULL;
atmosphere = &sky->atmosphere[i][iquad];
@@ -317,19 +317,19 @@ htsky_trace_ray
goto error;
}
- if(iband < sky->sw_bands_range[0] || iband > sky->sw_bands_range[1]) {
+ if(iband < sky->bands_range[0] || iband > sky->bands_range[1]) {
log_err(sky,
"%s: invalid spectral band index `%lu'. "
"Valid short wave spectral bands lie in [%lu, %lu]\n",
FUNC_NAME,
(unsigned long)iband,
- (unsigned long)sky->sw_bands_range[0],
- (unsigned long)sky->sw_bands_range[1]);
+ (unsigned long)sky->bands_range[0],
+ (unsigned long)sky->bands_range[1]);
res = RES_BAD_ARG;
goto error;
}
- if(iquad >= htsky_get_sw_spectral_band_quadrature_length(sky, iband)) {
+ if(iquad >= htsky_get_spectral_band_quadrature_length(sky, iband)) {
log_err(sky,
"invalid quadrature point `%lu' for the spectral band `%lu'.\n",
(unsigned long)iquad, (unsigned long)iband);
@@ -338,7 +338,7 @@ htsky_trace_ray
}
/* Fetch the clouds/atmosphere corresponding to the submitted spectral data */
- i = iband - sky->sw_bands_range[0];
+ i = iband - sky->bands_range[0];
clouds = sky->is_cloudy ? sky->clouds[i][iquad].octree : NULL;
atmosphere = sky->atmosphere[i][iquad].bitree;
