commit 351eb0beac9155956fccc84b704cc318b7dadbba
parent 94465b104794cf711ed3a42798eb6eaceea093c4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 27 Jan 2020 15:06:07 +0100
Refactor of the htsky_cloud.c file
Diffstat:
| M | src/htsky_cloud.c | | | 724 | ++++++++++++++++++++++++++++++++++++++++++++++--------------------------------- |
1 file changed, 419 insertions(+), 305 deletions(-)
diff --git a/src/htsky_cloud.c b/src/htsky_cloud.c
@@ -55,177 +55,431 @@ aabb_intersect
&& !(aabb0_upp[2] < aabb1_low[2]) && !(aabb0_low[2] > aabb1_upp[2]);
}
-static void
-cloud_vox_get_particle
- (const size_t xyz[3],
- float vox[],
- const struct build_tree_context* ctx)
+static res_T
+setup_svx2htcp_z_lut(struct htsky* sky, double* cloud_upp_z)
+{
+ double org, height;
+ double len_z;
+ size_t nsplits;
+ size_t iz, iz2;
+ res_T res = RES_OK;
+ ASSERT(sky && sky->htcp_desc.irregular_z && cloud_upp_z);
+
+ /* Find the min voxel size along Z and compute the length of a Z column */
+ len_z = 0;
+ sky->lut_cell_sz = DBL_MAX;
+ FOR_EACH(iz, 0, sky->htcp_desc.spatial_definition[2]) {
+ len_z += sky->htcp_desc.vxsz_z[iz];
+ sky->lut_cell_sz = MMIN(sky->lut_cell_sz, sky->htcp_desc.vxsz_z[iz]);
+ }
+
+ /* Allocate the svx2htcp LUT. This LUT is a regular table whose absolute
+ * size is greater or equal to a Z column in the htcp file. The size of its
+ * cells is the minimal voxel size in Z of the htcp file */
+ nsplits = (size_t)ceil(len_z / sky->lut_cell_sz);
+ 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");
+ goto error;
+ }
+
+ /* Setup the svx2htcp LUT. Each LUT entry stores the index of the current Z
+ * voxel in the htcp file that overlaps the entry lower bound as well as the
+ * lower bound in Z of the next htcp voxel. */
+ iz2 = 0;
+ org = sky->htcp_desc.lower[2];
+ height = org + sky->htcp_desc.vxsz_z[iz2];
+ FOR_EACH(iz, 0, nsplits) {
+ const double upp_z = (double)(iz + 1) * sky->lut_cell_sz + org;
+ darray_split_data_get(&sky->svx2htcp_z)[iz].index = iz2;
+ darray_split_data_get(&sky->svx2htcp_z)[iz].height = height;
+ if(upp_z >= height && iz + 1 < nsplits) {
+ ASSERT(iz2 + 1 < sky->htcp_desc.spatial_definition[2]);
+ height += sky->htcp_desc.vxsz_z[++iz2];
+ }
+ }
+ ASSERT(eq_eps(height - org, len_z, 1.e-6));
+
+exit:
+ *cloud_upp_z = height;
+ return res;
+error:
+ darray_split_clear(&sky->svx2htcp_z);
+ height = -1;
+ goto exit;
+}
+
+static INLINE void
+compute_k_bounds_regular_z
+ (const struct htsky* sky,
+ const size_t iband,
+ const double vox_svx_low[3], /* Lower bound of the SVX voxel */
+ const double vox_svx_upp[3], /* Upper bound of the SVX voxel */
+ double ka_bounds[2],
+ double ks_bounds[2],
+ double kext_bounds[2])
{
- const struct htcp_desc* htcp_desc;
size_t ivox[3];
size_t igrid_low[3], igrid_upp[3];
- double vox_low[3], vox_upp[3];
+ size_t i;
+ double ka, ks, kext;
double low[3], upp[3];
+ double ipart;
+ double rho_da; /* Dry air density */
double rct;
+ double Ca_avg;
+ double Cs_avg;
+ 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];
+
+ /* Fetch the optical properties of the spectral band */
+ Ca_avg = sky->sw_bands[i].Ca_avg;
+ Cs_avg = sky->sw_bands[i].Cs_avg;
+
+ /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
+ * the SVX voxel */
+ low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x;
+ low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x;
+ low[2] = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0];
+ upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y;
+ upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y;
+ upp[2] = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0];
+ igrid_low[0] = (size_t)low[0];
+ igrid_low[1] = (size_t)low[1];
+ igrid_low[2] = (size_t)low[2];
+ igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0);
+ igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0);
+ igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0);
+ ASSERT(igrid_low[0] <= igrid_upp[0]);
+ ASSERT(igrid_low[1] <= igrid_upp[1]);
+ ASSERT(igrid_low[2] <= igrid_upp[2]);
+
+ /* Prepare the iteration over the grid voxels overlapped by the SVX voxel */
+ ka_bounds[0] = ks_bounds[0] = kext_bounds[0] = DBL_MAX;
+ ka_bounds[1] = ks_bounds[1] = kext_bounds[1] =-DBL_MAX;
+
+ /* Iterate over the grid voxels overlapped by the SVX voxel */
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+ FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
+ /* Compute the radiative properties */
+ rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox);
+ rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
+ ka = Ca_avg * rho_da * rct;
+ ks = Cs_avg * rho_da * rct;
+ kext = ka + ks;
+ /* Update the boundaries of the radiative properties */
+ ka_bounds[0] = MMIN(ka_bounds[0], ka);
+ ka_bounds[1] = MMAX(ka_bounds[1], ka);
+ ks_bounds[0] = MMIN(ks_bounds[0], ks);
+ ks_bounds[1] = MMAX(ks_bounds[1], ks);
+ kext_bounds[0] = MMIN(kext_bounds[0], kext);
+ kext_bounds[1] = MMAX(kext_bounds[1], kext);
+ #ifndef NDEBUG
+ {
+ double tmp_low[3], tmp_upp[3];
+ htcp_desc_get_voxel_aabb
+ (&sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp);
+ ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_svx_low, vox_svx_upp));
+ }
+ #endif
+ }
+ }
+ }
+}
+
+
+static void
+compute_k_bounds_irregular_z
+ (const struct htsky* sky,
+ const size_t iband,
+ const double vox_svx_low[3], /* Lower bound of the SVX voxel */
+ const double vox_svx_upp[3], /* Upper bound of the SVX voxel */
+ double ka_bounds[2],
+ double ks_bounds[2],
+ double kext_bounds[2])
+{
+ size_t ivox[3];
+ size_t igrid_low[2], igrid_upp[2];
+ size_t i;
double ka, ks, kext;
- double ka_min, ka_max;
- double ks_min, ks_max;
- double kext_min, kext_max;
+ double low[2], upp[2];
+ double ipart;
double rho_da; /* Dry air density */
+ double rct;
double Ca_avg;
double Cs_avg;
- double ipart;
- size_t i;
- ASSERT(xyz && vox && ctx);
+ double pos_z;
+ size_t ilut_low, ilut_upp;
+ size_t ilut;
+ size_t ivox_z_prev = SIZE_MAX;
+ ASSERT(sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp);
+ ASSERT(ka_bounds && ks_bounds && kext_bounds);
- i = ctx->iband - ctx->sky->sw_bands_range[0];
- htcp_desc = &ctx->sky->htcp_desc;
+ i = iband - sky->sw_bands_range[0];
/* Fetch the optical properties of the spectral band */
- Ca_avg = ctx->sky->sw_bands[i].Ca_avg;
- Cs_avg = ctx->sky->sw_bands[i].Cs_avg;
+ Ca_avg = sky->sw_bands[i].Ca_avg;
+ Cs_avg = sky->sw_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 */
+ low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x;
+ low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x;
+ upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y;
+ upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y;
+ igrid_low[0] = (size_t)low[0];
+ igrid_low[1] = (size_t)low[1];
+ igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0);
+ igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0);
+ ASSERT(igrid_low[0] <= igrid_upp[0]);
+ ASSERT(igrid_low[1] <= igrid_upp[1]);
+
+ /* Compute the *inclusive* bounds of the indices of the LUT cells
+ * overlapped by the SVX voxel */
+ ilut_low = (size_t)((vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
+ ilut_upp = (size_t)((vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
+ ASSERT(ilut_low <= ilut_upp);
+
+ /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
+ ka_bounds[0] = ks_bounds[0] = kext_bounds[0] = DBL_MAX;
+ ka_bounds[1] = ks_bounds[1] = kext_bounds[1] =-DBL_MAX;
+ ivox_z_prev = SIZE_MAX;
+ pos_z = vox_svx_low[2];
+ ASSERT(pos_z >= (double)ilut_low * sky->lut_cell_sz);
+ ASSERT(pos_z <= (double)ilut_upp * sky->lut_cell_sz);
+
+ /* Iterate over the LUT cells overlapped by the voxel */
+ FOR_EACH(ilut, ilut_low, ilut_upp+1) {
+ const struct split* split = darray_split_cdata_get(&sky->svx2htcp_z)+ilut;
+ ASSERT(ilut < darray_split_size_get(&sky->svx2htcp_z));
+
+ ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
+ if(ivox[2] >= sky->htcp_desc.spatial_definition[2]
+ && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
+ ivox[2] = split->index;
+ }
+
+ /* Compute the upper bound of the *next* LUT cell clamped to the voxel
+ * upper bound. Note that the upper bound of the current LUT cell is
+ * the lower bound of the next cell, i.e. (ilut+1)*lut_cell_sz. The
+ * upper bound of the next cell is thus the lower bound of the cell
+ * following the next cell, i.e. (ilut+2)*lut_cell_sz */
+ pos_z = MMIN((double)(ilut+2)*sky->lut_cell_sz, vox_svx_upp[2]);
+
+ /* Does the LUT cell overlap an already handled LES voxel? */
+ if(ivox[2] == ivox_z_prev) continue;
+ ivox_z_prev = ivox[2];
+
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+
+ /* Compute the radiative properties */
+ rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox);
+ rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
+ ka = Ca_avg * rho_da * rct;
+ ks = Cs_avg * rho_da * rct;
+ kext = ka + ks;
+
+ /* Update the boundaries of the radiative properties */
+ ka_bounds[0] = MMIN(ka_bounds[0], ka);
+ ka_bounds[1] = MMAX(ka_bounds[1], ka);
+ ks_bounds[0] = MMIN(ks_bounds[0], ks);
+ ks_bounds[1] = MMAX(ks_bounds[1], ks);
+ kext_bounds[0] = MMIN(kext_bounds[0], kext);
+ kext_bounds[1] = MMAX(kext_bounds[1], kext);
+ }
+ }
+ }
+}
+
+static void
+cloud_vox_get_particle
+ (const size_t xyz[3],
+ float vox[],
+ const struct build_tree_context* ctx)
+{
+ double vox_low[3], vox_upp[3];
+ double ka[2], ks[2], kext[2];
+ ASSERT(xyz && vox && ctx);
/* Compute the AABB of the SVX voxel */
- vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + htcp_desc->lower[0];
- vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + htcp_desc->lower[1];
- vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + htcp_desc->lower[2];
+ vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + ctx->sky->htcp_desc.lower[0];
+ vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + ctx->sky->htcp_desc.lower[1];
+ vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + ctx->sky->htcp_desc.lower[2];
vox_upp[0] = vox_low[0] + ctx->vxsz[0];
vox_upp[1] = vox_low[1] + ctx->vxsz[1];
vox_upp[2] = vox_low[2] + ctx->vxsz[2];
+ if(!ctx->sky->htcp_desc.irregular_z) { /* 1 LES voxel <=> 1 SVX voxel */
+ compute_k_bounds_regular_z
+ (ctx->sky, ctx->iband, vox_low, vox_upp, ka, ks, kext);
+ } else {
+ ASSERT(xyz[2] < darray_split_size_get(&ctx->sky->svx2htcp_z));
+ compute_k_bounds_irregular_z
+ (ctx->sky, ctx->iband, vox_low, vox_upp, ka, ks, kext);
+ }
+
+ /* Ensure that the single precision bounds include their double precision
+ * version. */
+ if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX);
+ if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX);
+ if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX);
+ if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX);
+ if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX);
+ if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX);
+
+ vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ka, HTSKY_SVX_MIN, (float)ka[0]);
+ vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ka, HTSKY_SVX_MAX, (float)ka[1]);
+ vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ks, HTSKY_SVX_MIN, (float)ks[0]);
+ vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ks, HTSKY_SVX_MAX, (float)ks[1]);
+ vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Kext, HTSKY_SVX_MIN, (float)kext[0]);
+ vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Kext, HTSKY_SVX_MAX, (float)kext[1]);
+}
+
+static void
+compute_xh2o_range_regular_z
+ (const struct htsky* sky,
+ const double vox_svx_low[3], /* Lower bound of the SVX voxel */
+ const double vox_svx_upp[3], /* Upper bound of the SVX voxel */
+ double x_h2o_range[2])
+{
+ size_t ivox[3];
+ size_t igrid_low[3], igrid_upp[3];
+ double low[3], upp[3];
+ double ipart;
+ ASSERT(!sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp);
+ ASSERT(x_h2o_range);
+
/* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
* the SVX voxel in X and Y */
- low[0] = (vox_low[0]-htcp_desc->lower[0])/htcp_desc->vxsz_x;
- low[1] = (vox_low[1]-htcp_desc->lower[1])/htcp_desc->vxsz_x;
- upp[0] = (vox_upp[0]-htcp_desc->lower[0])/htcp_desc->vxsz_y;
- upp[1] = (vox_upp[1]-htcp_desc->lower[1])/htcp_desc->vxsz_y;
+ low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x;
+ low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x;
+ low[2] = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0];
+ upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y;
+ upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y;
+ upp[2] = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->htcp_desc.vxsz_z[0];
igrid_low[0] = (size_t)low[0];
igrid_low[1] = (size_t)low[1];
+ igrid_low[2] = (size_t)low[2];
igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0);
igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0);
+ igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0);
ASSERT(igrid_low[0] <= igrid_upp[0]);
ASSERT(igrid_low[1] <= igrid_upp[1]);
-
- if(!ctx->sky->htcp_desc.irregular_z) { /* 1 LES voxel <=> 1 SVX voxel */
- /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
- * the SVX voxel along the Z axis */
- low[2] = (vox_low[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
- upp[2] = (vox_upp[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
- igrid_low[2] = (size_t)low[2];
- igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0);
- ASSERT(igrid_low[2] <= igrid_upp[2]);
-
- /* Prepare the iteration over the grid voxels overlapped by the SVX voxel */
- ka_min = ks_min = kext_min = DBL_MAX;
- ka_max = ks_max = kext_max =-DBL_MAX;
-
- /* Iterate over the grid voxels overlapped by the SVX voxel */
- FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
- FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
- FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
- /* Compute the radiative properties */
- rho_da = cloud_dry_air_density(htcp_desc, ivox);
- rct = htcp_desc_RCT_at(htcp_desc, ivox[0], ivox[1], ivox[2], 0);
- ka = Ca_avg * rho_da * rct;
- ks = Cs_avg * rho_da * rct;
- kext = ka + ks;
- /* Update the boundaries of the radiative properties */
- ka_min = MMIN(ka_min, ka);
- ka_max = MMAX(ka_max, ka);
- ks_min = MMIN(ks_min, ks);
- ks_max = MMAX(ks_max, ks);
- kext_min = MMIN(kext_min, kext);
- kext_max = MMAX(kext_max, kext);
- #ifndef NDEBUG
- {
- double tmp_low[3], tmp_upp[3];
- htcp_desc_get_voxel_aabb
- (&ctx->sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp);
- ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_low, vox_upp));
- }
- #endif
+ ASSERT(igrid_low[2] <= igrid_upp[2]);
+
+ /* Prepare the iteration overt the grid voxels overlapped by the SVX voxel */
+ x_h2o_range[0] = DBL_MAX;
+ x_h2o_range[1] =-DBL_MAX;
+
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+ FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
+ double x_h2o;
+
+ /* Compute the xH2O for the current LES voxel */
+ x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
+
+ /* Update the xH2O range of the SVX voxel */
+ x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
+ x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
+ #ifndef NDEBUG
+ {
+ double tmp_low[3], tmp_upp[3];
+ htcp_desc_get_voxel_aabb
+ (&sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp);
+ ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_svx_low, vox_svx_upp));
}
+ #endif
}
}
- } else {
- double pos_z;
- size_t ilut_low, ilut_upp;
- size_t ilut;
- size_t ivox_z_prev = SIZE_MAX;
-
- /* Compute the *inclusive* bounds of the indices of the LUT cells
- * overlapped by the SVX voxel */
- ilut_low = (size_t)((vox_low[2]-htcp_desc->lower[2])/ctx->sky->lut_cell_sz);
- ilut_upp = (size_t)((vox_upp[2]-htcp_desc->lower[2])/ctx->sky->lut_cell_sz);
- ASSERT(ilut_low <= ilut_upp);
-
- /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
- ka_min = ks_min = kext_min = DBL_MAX;
- ka_max = ks_max = kext_max =-DBL_MAX;
- ivox_z_prev = SIZE_MAX;
- pos_z = vox_low[2];
- ASSERT(pos_z >= (double)ilut_low * ctx->sky->lut_cell_sz);
- ASSERT(pos_z <= (double)ilut_upp * ctx->sky->lut_cell_sz);
-
- /* Iterate over the LUT cells overlapped by the voxel */
- FOR_EACH(ilut, ilut_low, ilut_upp+1) {
- const struct split* split = darray_split_cdata_get(&ctx->sky->svx2htcp_z)+ilut;
- ASSERT(ilut < darray_split_size_get(&ctx->sky->svx2htcp_z));
-
- ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
- if(ivox[2] >= ctx->sky->htcp_desc.spatial_definition[2]
- && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
- ivox[2] = split->index;
- }
+ }
+}
- /* Compute the upper bound of the *next* LUT cell clamped to the voxel
- * upper bound. Note that the upper bound of the current LUT cell is
- * the lower bound of the next cell, i.e. (ilut+1)*lut_cell_sz. The
- * upper bound of the next cell is thus the lower bound of the cell
- * following the next cell, i.e. (ilut+2)*lut_cell_sz */
- pos_z = MMIN((double)(ilut+2)*ctx->sky->lut_cell_sz, vox_upp[2]);
-
- /* Does the LUT cell overlap an already handled LES voxel? */
- if(ivox[2] == ivox_z_prev) continue;
- ivox_z_prev = ivox[2];
-
- FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
- FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
-
- /* Compute the radiative properties */
- rho_da = cloud_dry_air_density(htcp_desc, ivox);
- rct = htcp_desc_RCT_at(htcp_desc, ivox[0], ivox[1], ivox[2], 0);
- ka = Ca_avg * rho_da * rct;
- ks = Cs_avg * rho_da * rct;
- kext = ka + ks;
-
- /* Update the boundaries of the radiative properties */
- ka_min = MMIN(ka_min, ka);
- ka_max = MMAX(ka_max, ka);
- ks_min = MMIN(ks_min, ks);
- ks_max = MMAX(ks_max, ks);
- kext_min = MMIN(kext_min, kext);
- kext_max = MMAX(kext_max, kext);
- }
+static void
+compute_xh2o_range_irregular_z
+ (const struct htsky* sky,
+ const double vox_svx_low[3], /* Lower bound of the SVX voxel */
+ const double vox_svx_upp[3], /* Upper bound of the SVX voxel */
+ double x_h2o_range[2])
+{
+ size_t ivox[3];
+ size_t igrid_low[2], igrid_upp[2];
+ size_t ilut_low, ilut_upp;
+ size_t ilut;
+ size_t ivox_z_prev;
+ double low[2], upp[2];
+ double pos_z;
+ double ipart;
+ ASSERT(sky->htcp_desc.irregular_z && vox_svx_low && vox_svx_upp);
+ ASSERT(x_h2o_range);
+
+ /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
+ * the SVX voxel in X and Y */
+ low[0] = (vox_svx_low[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_x;
+ low[1] = (vox_svx_low[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_x;
+ upp[0] = (vox_svx_upp[0]-sky->htcp_desc.lower[0])/sky->htcp_desc.vxsz_y;
+ upp[1] = (vox_svx_upp[1]-sky->htcp_desc.lower[1])/sky->htcp_desc.vxsz_y;
+ igrid_low[0] = (size_t)low[0];
+ igrid_low[1] = (size_t)low[1];
+ igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0);
+ igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0);
+ ASSERT(igrid_low[0] <= igrid_upp[0]);
+ ASSERT(igrid_low[1] <= igrid_upp[1]);
+
+ /* Compute the *inclusive* bounds of the indices of the LUT cells
+ * overlapped by the SVX voxel */
+ ilut_low = (size_t)((vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
+ ilut_upp = (size_t)((vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
+ ASSERT(ilut_low <= ilut_upp);
+
+ /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
+ x_h2o_range[0] = DBL_MAX;
+ x_h2o_range[1] =-DBL_MAX;
+ ivox_z_prev = SIZE_MAX;
+ pos_z = vox_svx_low[2];
+ ASSERT(pos_z >= (double)ilut_low * sky->lut_cell_sz);
+ ASSERT(pos_z <= (double)ilut_upp * sky->lut_cell_sz);
+
+ /* Iterate over the LUT cells overlapped by the voxel */
+ FOR_EACH(ilut, ilut_low, ilut_upp+1) {
+ const struct split* split = darray_split_cdata_get(&sky->svx2htcp_z) + ilut;
+ ASSERT(ilut < darray_split_size_get(&sky->svx2htcp_z));
+
+ ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
+ if(ivox[2] >= sky->htcp_desc.spatial_definition[2]
+ && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
+ ivox[2] = split->index;
+ }
+
+ /* Compute the upper bound of the *next* LUT cell clamped to the voxel
+ * upper bound. Note that the upper bound of the current LUT cell is
+ * the lower bound of the next cell, i.e. (ilut+1)*lut_cell_sz. The
+ * upper bound of the next cell is thus the lower bound of the cell
+ * following the next cell, i.e. (ilut+2)*lut_cell_sz */
+ pos_z = MMIN((double)(ilut+2)*sky->lut_cell_sz, vox_svx_upp[2]);
+
+ /* Does the LUT voxel overlap an already handled LES voxel? */
+ if(ivox[2] == ivox_z_prev) continue;
+ ivox_z_prev = ivox[2];
+
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+ double x_h2o;
+
+ /* Compute the xH2O for the current LES voxel */
+ x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
+
+ /* Update the xH2O range of the SVX voxel */
+ x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
+ x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
}
}
}
-
- /* Ensure that the single precision bounds include their double precision
- * version. */
- if(ka_min != (float)ka_min) ka_min = nextafterf((float)ka_min,-FLT_MAX);
- if(ka_max != (float)ka_max) ka_max = nextafterf((float)ka_max, FLT_MAX);
- if(ks_min != (float)ks_min) ks_min = nextafterf((float)ks_min,-FLT_MAX);
- if(ks_max != (float)ks_max) ks_max = nextafterf((float)ks_max, FLT_MAX);
- if(kext_min != (float)kext_min) kext_min = nextafterf((float)kext_min,-FLT_MAX);
- if(kext_max != (float)kext_max) kext_max = nextafterf((float)kext_max, FLT_MAX);
-
- vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ka, HTSKY_SVX_MIN, (float)ka_min);
- vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ka, HTSKY_SVX_MAX, (float)ka_max);
- vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ks, HTSKY_SVX_MIN, (float)ks_min);
- vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Ks, HTSKY_SVX_MAX, (float)ks_max);
- vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Kext, HTSKY_SVX_MIN, (float)kext_min);
- vox_set(vox, HTSKY_CPNT_PARTICLES, HTSKY_Kext, HTSKY_SVX_MAX, (float)kext_max);
}
static void
@@ -235,20 +489,13 @@ cloud_vox_get_gas
const struct build_tree_context* ctx)
{
const struct htcp_desc* htcp_desc;
- size_t ivox[3];
- size_t igrid_low[3], igrid_upp[3];
struct htgop_layer layer;
struct htgop_layer_sw_spectral_interval band;
size_t ilayer;
size_t layer_range[2];
double x_h2o_range[2];
- double low[3], upp[3];
double vox_low[3], vox_upp[3]; /* Voxel AABB */
- double ka_min, ka_max;
- double ks_min, ks_max;
- double kext_min, kext_max;
- double x_h2o;
- double ipart;
+ double ka[2], ks[2], kext[2];
ASSERT(xyz && vox && ctx);
htcp_desc = &ctx->sky->htcp_desc;
@@ -261,117 +508,20 @@ cloud_vox_get_gas
vox_upp[1] = vox_low[1] + ctx->vxsz[1];
vox_upp[2] = vox_low[2] + ctx->vxsz[2];
- /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
- * the SVX voxel in X and Y */
- low[0] = (vox_low[0]-htcp_desc->lower[0])/htcp_desc->vxsz_x;
- low[1] = (vox_low[1]-htcp_desc->lower[1])/htcp_desc->vxsz_x;
- upp[0] = (vox_upp[0]-htcp_desc->lower[0])/htcp_desc->vxsz_y;
- upp[1] = (vox_upp[1]-htcp_desc->lower[1])/htcp_desc->vxsz_y;
- igrid_low[0] = (size_t)low[0];
- igrid_low[1] = (size_t)low[1];
- igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0);
- igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0);
- ASSERT(igrid_low[0] <= igrid_upp[0]);
- ASSERT(igrid_low[1] <= igrid_upp[1]);
-
/* Define the xH2O range from the LES data */
if(!ctx->sky->htcp_desc.irregular_z) { /* 1 LES voxel <=> 1 SVX voxel */
- /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
- * the SVX voxel in Z */
- low[2] = (vox_low[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
- upp[2] = (vox_upp[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
- igrid_low[2] = (size_t)low[2];
- igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0);
- ASSERT(igrid_low[2] <= igrid_upp[2]);
-
- /* Prepare the iteration overt the grid voxels overlapped by the SVX voxel */
- x_h2o_range[0] = DBL_MAX;
- x_h2o_range[1] =-DBL_MAX;
-
- FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
- FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
- FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
-
- /* Compute the xH2O for the current LES voxel */
- x_h2o = cloud_water_vapor_molar_fraction(htcp_desc, ivox);
-
- /* Update the xH2O range of the SVX voxel */
- x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
- x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
- #ifndef NDEBUG
- {
- double tmp_low[3], tmp_upp[3];
- htcp_desc_get_voxel_aabb
- (&ctx->sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp);
- ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_low, vox_upp));
- }
- #endif
- }
- }
- }
+ compute_xh2o_range_regular_z(ctx->sky, vox_low, vox_upp, x_h2o_range);
} else { /* A SVX voxel can be overlapped by 2 LES voxels */
- double pos_z;
- size_t ilut_low, ilut_upp;
- size_t ilut;
- size_t ivox_z_prev;
ASSERT(xyz[2] < darray_split_size_get(&ctx->sky->svx2htcp_z));
-
- /* Compute the *inclusive* bounds of the indices of the LUT cells
- * overlapped by the SVX voxel */
- ilut_low = (size_t)((vox_low[2] - htcp_desc->lower[2]) / ctx->sky->lut_cell_sz);
- ilut_upp = (size_t)((vox_upp[2] - htcp_desc->lower[2]) / ctx->sky->lut_cell_sz);
- ASSERT(ilut_low <= ilut_upp);
-
- /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
- x_h2o_range[0] = DBL_MAX;
- x_h2o_range[1] =-DBL_MAX;
- ivox_z_prev = SIZE_MAX;
- pos_z = vox_low[2];
- ASSERT(pos_z >= (double)ilut_low * ctx->sky->lut_cell_sz);
- ASSERT(pos_z <= (double)ilut_upp * ctx->sky->lut_cell_sz);
-
- /* Iterate over the LUT cells overlapped by the voxel */
- FOR_EACH(ilut, ilut_low, ilut_upp+1) {
- const struct split* split = darray_split_cdata_get(&ctx->sky->svx2htcp_z)+ilut;
- ASSERT(ilut < darray_split_size_get(&ctx->sky->svx2htcp_z));
-
- ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
- if(ivox[2] >= ctx->sky->htcp_desc.spatial_definition[2]
- && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
- ivox[2] = split->index;
- }
-
- /* Compute the upper bound of the *next* LUT cell clamped to the voxel
- * upper bound. Note that the upper bound of the current LUT cell is
- * the lower bound of the next cell, i.e. (ilut+1)*lut_cell_sz. The
- * upper bound of the next cell is thus the lower bound of the cell
- * following the next cell, i.e. (ilut+2)*lut_cell_sz */
- pos_z = MMIN((double)(ilut+2)*ctx->sky->lut_cell_sz, vox_upp[2]);
-
- /* Does the LUT voxel overlap an already handled LES voxel? */
- if(ivox[2] == ivox_z_prev) continue;
- ivox_z_prev = ivox[2];
-
- FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
- FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
-
- /* Compute the xH2O for the current LES voxel */
- x_h2o = cloud_water_vapor_molar_fraction(&ctx->sky->htcp_desc, ivox);
-
- /* Update the xH2O range of the SVX voxel */
- x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
- x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
- }
- }
- }
+ compute_xh2o_range_irregular_z(ctx->sky, vox_low, vox_upp, x_h2o_range);
}
/* Define the atmospheric layers overlapped by the SVX voxel */
HTGOP(position_to_layer_id(ctx->sky->htgop, vox_low[2], &layer_range[0]));
HTGOP(position_to_layer_id(ctx->sky->htgop, vox_upp[2], &layer_range[1]));
- ka_min = ks_min = kext_min = DBL_MAX;
- ka_max = ks_max = kext_max =-DBL_MAX;
+ ka[0] = ks[0] = kext[0] = DBL_MAX;
+ 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) {
@@ -387,33 +537,33 @@ cloud_vox_get_gas
/* ... 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_min = MMIN(ka_min, k[0]);
- ka_max = MMAX(ka_max, k[1]);
+ 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_min = MMIN(ks_min, k[0]);
- ks_max = MMAX(ks_max, k[1]);
+ 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_min = MMIN(kext_min, k[0]);
- kext_max = MMAX(kext_max, k[1]);
+ kext[0] = MMIN(kext[0], k[0]);
+ kext[1] = MMAX(kext[1], k[1]);
}
/* Ensure that the single precision bounds include their double precision
* version. */
- if(ka_min != (float)ka_min) ka_min = nextafterf((float)ka_min,-FLT_MAX);
- if(ka_max != (float)ka_max) ka_max = nextafterf((float)ka_max, FLT_MAX);
- if(ks_min != (float)ks_min) ks_min = nextafterf((float)ks_min,-FLT_MAX);
- if(ks_max != (float)ks_max) ks_max = nextafterf((float)ks_max, FLT_MAX);
- if(kext_min != (float)kext_min) kext_min = nextafterf((float)kext_min,-FLT_MAX);
- if(kext_max != (float)kext_max) kext_max = nextafterf((float)kext_max, FLT_MAX);
-
- vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MIN, (float)ka_min);
- vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MAX, (float)ka_max);
- vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MIN, (float)ks_min);
- vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MAX, (float)ks_max);
- vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MIN, (float)kext_min);
- vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MAX, (float)kext_max);
+ if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX);
+ if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX);
+ if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX);
+ if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX);
+ if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX);
+ if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX);
+
+ vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MIN, (float)ka[0]);
+ vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MAX, (float)ka[1]);
+ vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MIN, (float)ks[0]);
+ vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MAX, (float)ks[1]);
+ vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MIN, (float)kext[0]);
+ vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MAX, (float)kext[1]);
}
static void
@@ -495,45 +645,9 @@ cloud_setup
/* Regular voxel size along the Z dimension: compute its upper boundary as
* the others dimensions */
upp[2] = low[2] + (double)raw_def[2] * sky->htcp_desc.vxsz_z[0];
-
- /* TODO move the following block in a separate function */
} else { /* Irregular voxel size along Z */
- double len_z;
- size_t nsplits;
- size_t iz, iz2;;
-
- /* Find the min voxel size along Z and compute the length of a Z column */
- len_z = 0;
- sky->lut_cell_sz = DBL_MAX;
- FOR_EACH(iz, 0, sky->htcp_desc.spatial_definition[2]) {
- len_z += sky->htcp_desc.vxsz_z[iz];
- sky->lut_cell_sz = MMIN(sky->lut_cell_sz, sky->htcp_desc.vxsz_z[iz]);
- }
- /* Allocate the svx2htcp LUT. This LUT is a regular table whose absolute
- * size is greater or equal to a Z column in the htcp file. The size of its
- * cells is the minimal voxel size in Z of the htcp file */
- nsplits = (size_t)ceil(len_z / sky->lut_cell_sz);
- 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");
- goto error;
- }
- /* Setup the svx2htcp LUT. Each LUT entry stores the index of the current Z
- * voxel in the htcp file that overlaps the entry lower bound as well as the
- * lower bound in Z of the next htcp voxel. */
- iz2 = 0;
- upp[2] = low[2] + sky->htcp_desc.vxsz_z[iz2];
- FOR_EACH(iz, 0, nsplits) {
- const double upp_z = (double)(iz + 1) * sky->lut_cell_sz + low[2];
- darray_split_data_get(&sky->svx2htcp_z)[iz].index = iz2;
- darray_split_data_get(&sky->svx2htcp_z)[iz].height = upp[2];
- if(upp_z >= upp[2] && iz + 1 < nsplits) {
- ASSERT(iz2 + 1 < sky->htcp_desc.spatial_definition[2]);
- upp[2] += sky->htcp_desc.vxsz_z[++iz2];
- }
- }
- ASSERT(eq_eps(upp[2] - low[2], len_z, 1.e-6));
+ res = setup_svx2htcp_z_lut(sky, &upp[2]);
+ if(res != RES_OK) goto error;
}
/* Setup the build context */
@@ -554,7 +668,7 @@ cloud_setup
goto error;
}
- /* Compute how many octree are going to be built */
+ /* 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];
