commit b8624878c2def4826afd1afd90c2a0d49908bdd3
parent c8c0b159392054cb6fb15a28a36662ff6da7a49e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 12 Mar 2020 15:36:12 +0100
Add long waves support
Diffstat:
7 files changed, 524 insertions(+), 205 deletions(-)
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,48 +60,55 @@ 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;
}
@@ -123,6 +135,139 @@ sample_sw_spectral_data
HTGOP(spectral_interval_sample_quadrature(&specint, r1, 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])
+{
+ 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
setup_cache_stream
(struct htsky* sky,
@@ -190,18 +335,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;
+ double bands_range[2];
/* Read the cache header */
#define READ(Var, N) { \
@@ -217,12 +368,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 +402,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:
@@ -269,7 +441,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 +496,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) {
@@ -368,9 +540,19 @@ htsky_create
goto error;
}
- /* 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;
+ /* 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]) {
+ res = htgop_get_sw_spectral_intervals_CIE_XYZ(sky->htgop, 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;
+ }
/* Setup the atmopshere */
time_current(&t0);
@@ -407,7 +589,8 @@ htsky_create
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 */
@@ -484,14 +667,14 @@ 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;
}
}
@@ -518,11 +701,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);
@@ -580,111 +763,139 @@ htsky_fetch_raw_property
}
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,15 +904,26 @@ 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;
}
+int
+htsky_is_long_wave(const struct htsky* htsky)
+{
+ ASSERT(htsky);
+ return htsky->is_long_wave;
+}
+
void
htsky_sample_sw_spectral_data_CIE_1931_X
(const struct htsky* sky,
@@ -710,6 +932,7 @@ htsky_sample_sw_spectral_data_CIE_1931_X
size_t* ispectral_band,
size_t* iquadrature_pt)
{
+ ASSERT(!sky->is_long_wave);
sample_sw_spectral_data
(sky->htgop, r0, r1, htgop_sample_sw_spectral_interval_CIE_1931_X,
ispectral_band, iquadrature_pt);
@@ -723,6 +946,7 @@ htsky_sample_sw_spectral_data_CIE_1931_Y
size_t* ispectral_band,
size_t* iquadrature_pt)
{
+ ASSERT(!sky->is_long_wave);
sample_sw_spectral_data
(sky->htgop, r0, r1, htgop_sample_sw_spectral_interval_CIE_1931_Y,
ispectral_band, iquadrature_pt);
@@ -736,6 +960,7 @@ htsky_sample_sw_spectral_data_CIE_1931_Z
size_t* ispectral_band,
size_t* iquadrature_pt)
{
+ ASSERT(!sky->is_long_wave);
sample_sw_spectral_data
(sky->htgop, r0, r1, htgop_sample_sw_spectral_interval_CIE_1931_Z,
ispectral_band, iquadrature_pt);
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 */ \
@@ -142,6 +144,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,25 +182,33 @@ 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_lw_spectral_bands_count
+ (const struct htsky* sky);
+
+HTSKY_API size_t
+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 int
+htsky_is_long_wave
+ (const struct htsky* sky);
+
HTSKY_API void
htsky_sample_sw_spectral_data_CIE_1931_X
(const struct htsky* sky,
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,7 +196,7 @@ 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) {
@@ -183,9 +209,13 @@ 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]));
@@ -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
@@ -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]);
@@ -849,7 +881,7 @@ 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,
@@ -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_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;
