htsky_atmosphere.c (10853B)
1 /* Copyright (C) 2018, 2019, 2020, 2021 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2018, 2019 Centre National de la Recherche Scientifique 3 * Copyright (C) 2018, 2019 Université Paul Sabatier 4 * 5 * This program is free software: you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation, either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 17 18 #define _POSIX_C_SOURCE 200809L /* nextafterf */ 19 20 #include "htsky_c.h" 21 #include "htsky_atmosphere.h" 22 #include "htsky_log.h" 23 #include "htsky_svx.h" 24 25 #include <high_tune/htgop.h> 26 #include <star/svx.h> 27 28 #include <math.h> 29 30 /******************************************************************************* 31 * Helper functions 32 ******************************************************************************/ 33 static void 34 atmosphere_vox_get 35 (const size_t xyz[3], 36 const uint64_t mcode, 37 void* dst, 38 void* context) 39 { 40 float* vox = dst; 41 struct build_tree_context* ctx = context; 42 struct htgop_level level; 43 size_t layer_range[2]; 44 size_t nlevels; 45 double vox_low, vox_upp; 46 double ka_min, ks_min, kext_min; 47 double ka_max, ks_max, kext_max; 48 size_t ilayer; 49 ASSERT(xyz && dst && context); 50 (void)mcode; 51 52 /* Compute the boundaries of the SVX voxel */ 53 HTGOP(get_level(ctx->sky->htgop, 0, &level)); 54 vox_low = (double)xyz[2] * ctx->vxsz[2] + level.height; 55 HTGOP(get_levels_count(ctx->sky->htgop, &nlevels)); 56 HTGOP(get_level(ctx->sky->htgop, nlevels-1, &level)); 57 vox_upp = MMIN(vox_low + ctx->vxsz[2], level.height); /* Handle numerical issues */ 58 59 /* Define the atmospheric layers overlapped by the SVX voxel */ 60 HTGOP(position_to_layer_id(ctx->sky->htgop, vox_low, &layer_range[0])); 61 HTGOP(position_to_layer_id(ctx->sky->htgop, vox_upp, &layer_range[1])); 62 63 ka_min = ks_min = kext_min = DBL_MAX; 64 ka_max = ks_max = kext_max =-DBL_MAX; 65 66 /* For each atmospheric layer that overlaps the SVX voxel ... */ 67 if(ctx->sky->spectral_type == HTSKY_SPECTRAL_LW) { 68 FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { 69 struct htgop_layer layer; 70 struct htgop_layer_lw_spectral_interval band; 71 size_t iquad; 72 73 HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer)); 74 75 /* ... retrieve the considered spectral interval */ 76 HTGOP(layer_get_lw_spectral_interval(&layer, ctx->iband, &band)); 77 78 /* ... and update the radiative properties bound with the per quadrature 79 * point nominal values */ 80 ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]); 81 ASSERT(ctx->quadrature_range[1] < band.quadrature_length); 82 FOR_EACH(iquad, ctx->quadrature_range[0], ctx->quadrature_range[1]+1) { 83 ka_min = MMIN(ka_min, band.ka_nominal[iquad]); 84 ka_max = MMAX(ka_max, band.ka_nominal[iquad]); 85 } 86 } 87 ks_min = 0; 88 ks_max = 0; 89 kext_min = ka_min; 90 kext_max = ka_max; 91 } else { 92 ASSERT(ctx->sky->spectral_type == HTSKY_SPECTRAL_SW); 93 94 FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { 95 struct htgop_layer layer; 96 struct htgop_layer_sw_spectral_interval band; 97 size_t iquad; 98 99 HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer)); 100 101 /* ... retrieve the considered spectral interval */ 102 HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band)); 103 104 /* ... and update the radiative properties bound with the per quadrature 105 * point nominal values */ 106 ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]); 107 ASSERT(ctx->quadrature_range[1] < band.quadrature_length); 108 FOR_EACH(iquad, ctx->quadrature_range[0], ctx->quadrature_range[1]+1) { 109 ka_min = MMIN(ka_min, band.ka_nominal[iquad]); 110 ka_max = MMAX(ka_max, band.ka_nominal[iquad]); 111 ks_min = MMIN(ks_min, band.ks_nominal[iquad]); 112 ks_max = MMAX(ks_max, band.ks_nominal[iquad]); 113 kext_min = MMIN(kext_min, band.ka_nominal[iquad]+band.ks_nominal[iquad]); 114 kext_max = MMAX(kext_max, band.ka_nominal[iquad]+band.ks_nominal[iquad]); 115 } 116 } 117 } 118 119 /* Ensure that the single precision bounds include their double precision 120 * version. */ 121 if(ka_min != (float)ka_min) ka_min = nextafterf((float)ka_min,-FLT_MAX); 122 if(ka_max != (float)ka_max) ka_max = nextafterf((float)ka_max, FLT_MAX); 123 if(ks_min != (float)ks_min) ks_min = nextafterf((float)ks_min,-FLT_MAX); 124 if(ks_max != (float)ks_max) ks_max = nextafterf((float)ks_max, FLT_MAX); 125 if(kext_min != (float)kext_min) kext_min = nextafterf((float)kext_min,-FLT_MAX); 126 if(kext_max != (float)kext_max) kext_max = nextafterf((float)kext_max, FLT_MAX); 127 128 /* Setup gas optical properties of the voxel */ 129 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MIN, (float)ka_min); 130 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ka, HTSKY_SVX_MAX, (float)ka_max); 131 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MIN, (float)ks_min); 132 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Ks, HTSKY_SVX_MAX, (float)ks_max); 133 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MIN, (float)kext_min); 134 vox_set(vox, HTSKY_CPNT_GAS, HTSKY_Kext, HTSKY_SVX_MAX, (float)kext_max); 135 } 136 137 static void 138 atmosphere_vox_merge 139 (void* dst, 140 const void* voxs[], 141 const size_t nvoxs, 142 void* ctx) 143 { 144 ASSERT(dst && voxs && nvoxs); 145 (void)ctx; 146 vox_merge_component(dst, HTSKY_CPNT_GAS, (const float**)voxs, nvoxs); 147 } 148 149 static int 150 atmosphere_vox_challenge_merge 151 (const struct svx_voxel voxs[], 152 const size_t nvoxs, 153 void* ctx) 154 { 155 ASSERT(voxs); 156 return vox_challenge_merge_component(HTSKY_CPNT_GAS, voxs, nvoxs, ctx); 157 } 158 159 /******************************************************************************* 160 * Local functions 161 ******************************************************************************/ 162 res_T 163 atmosphere_setup(struct htsky* sky, const double optical_thickness_threshold) 164 { 165 struct build_tree_context ctx = BUILD_TREE_CONTEXT_NULL; 166 struct htgop_level lvl_low, lvl_upp; 167 struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL; 168 double low, upp; 169 size_t nlayers, nlevels; 170 size_t definition; 171 size_t nbands; 172 size_t i; 173 res_T res = RES_OK; 174 ASSERT(sky && optical_thickness_threshold >= 0); 175 176 HTGOP(get_layers_count(sky->htgop, &nlayers)); 177 HTGOP(get_levels_count(sky->htgop, &nlevels)); 178 HTGOP(get_level(sky->htgop, 0, &lvl_low)); 179 HTGOP(get_level(sky->htgop, nlevels-1, &lvl_upp)); 180 low = lvl_low.height; 181 upp = lvl_upp.height; 182 definition = nlayers; 183 184 /* Setup the build context */ 185 ctx.sky = sky; 186 ctx.tau_threshold = optical_thickness_threshold; 187 ctx.vxsz[0] = INF; 188 ctx.vxsz[1] = INF; 189 ctx.vxsz[2] = (upp-low)/(double)definition; 190 191 /* Setup the voxel descriptor for the atmosphere. Note that in contrats with 192 * the clouds, the voxel contains only NFLOATS_PER_CPNT floats and not 193 * NFLOATS_PER_VOXEL. Indeed, the atmosphere has only 1 component: the gas. 194 * Anyway, we still rely on the memory layout of the cloud voxels excepted 195 * that we assume that the optical properties of the particles are never 196 * fetched. We thus have to ensure that the gas properties are arranged 197 * before the particles, i.e. HTSKY_CPNT_GAS == 0 */ 198 vox_desc.get = atmosphere_vox_get; 199 vox_desc.merge = atmosphere_vox_merge; 200 vox_desc.challenge_merge = atmosphere_vox_challenge_merge; 201 vox_desc.context = &ctx; 202 vox_desc.size = sizeof(float) * NFLOATS_PER_CPNT; 203 { STATIC_ASSERT(HTSKY_CPNT_GAS == 0, Unexpected_enum_value); } 204 205 /* Create as many atmospheric data structure than considered SW spectral 206 * bands */ 207 nbands = htsky_get_spectral_bands_count(sky); 208 sky->atmosphere = MEM_CALLOC 209 (sky->allocator, nbands, sizeof(*sky->atmosphere)); 210 if(!sky->atmosphere) { 211 log_err(sky, 212 "Could not create the list of per SW band atmospheric data structure.\n"); 213 res = RES_MEM_ERR; 214 goto error; 215 } 216 217 FOR_EACH(i, 0, nbands) { 218 size_t iquad; 219 struct htgop_spectral_interval band; 220 ctx.iband = i + sky->bands_range[0]; 221 222 switch(sky->spectral_type) { 223 case HTSKY_SPECTRAL_LW: 224 HTGOP(get_lw_spectral_interval(sky->htgop, ctx.iband, &band)); 225 break; 226 case HTSKY_SPECTRAL_SW: 227 HTGOP(get_sw_spectral_interval(sky->htgop, ctx.iband, &band)); 228 break; 229 default: FATAL("Unreachable code.\n"); break; 230 } 231 232 sky->atmosphere[i] = MEM_CALLOC(sky->allocator, 233 band.quadrature_length, sizeof(*sky->atmosphere[i])); 234 if(!sky->atmosphere[i]) { 235 log_err(sky, 236 "Could not create the list of per quadrature point atmospheric data " 237 "for the band %lu.\n", (unsigned long)ctx.iband); 238 res = RES_MEM_ERR; 239 goto error; 240 } 241 242 /* Build an atmospheric binary tree for each quadrature point of the 243 * considered spectral band */ 244 FOR_EACH(iquad, 0, band.quadrature_length) { 245 ctx.quadrature_range[0] = iquad; 246 ctx.quadrature_range[1] = iquad; 247 248 /* Create the atmospheric binary tree */ 249 res = svx_bintree_create(sky->svx, low, upp, definition, SVX_AXIS_Z, 250 &vox_desc, &sky->atmosphere[i][iquad].bitree); 251 if(res != RES_OK) { 252 log_err(sky, "Could not create the binary tree of the " 253 "atmospheric properties for the band %lu.\n", (unsigned long)ctx.iband); 254 goto error; 255 } 256 257 /* Fetch the binary tree descriptor for future use */ 258 SVX(tree_get_desc(sky->atmosphere[i][iquad].bitree, 259 &sky->atmosphere[i][iquad].bitree_desc)); 260 } 261 } 262 263 exit: 264 return res; 265 error: 266 atmosphere_clean(sky); 267 goto exit; 268 } 269 270 void 271 atmosphere_clean(struct htsky* sky) 272 { 273 size_t nbands; 274 size_t i; 275 ASSERT(sky); 276 277 if(!sky->atmosphere) return; 278 279 nbands = htsky_get_spectral_bands_count(sky); 280 FOR_EACH(i, 0, nbands) { 281 struct htgop_spectral_interval band; 282 size_t iband; 283 size_t iquad; 284 285 if(!sky->atmosphere[i]) continue; 286 287 iband = sky->bands_range[0] + i; 288 switch(sky->spectral_type) { 289 case HTSKY_SPECTRAL_LW: 290 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); 291 break; 292 case HTSKY_SPECTRAL_SW: 293 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); 294 break; 295 default: FATAL("Unreachable code.\n"); break; 296 } 297 298 FOR_EACH(iquad, 0, band.quadrature_length) { 299 if(sky->atmosphere[i][iquad].bitree) { 300 SVX(tree_ref_put(sky->atmosphere[i][iquad].bitree)); 301 sky->atmosphere[i][iquad].bitree = NULL; 302 } 303 } 304 MEM_RM(sky->allocator, sky->atmosphere[i]); 305 } 306 MEM_RM(sky->allocator, sky->atmosphere); 307 sky->atmosphere = NULL; 308 } 309