htsky

Load and structure a vertically stratified atmosphere
git clone git://git.meso-star.fr/htsky.git
Log | Files | Refs | README | LICENSE

commit 6465dbb2c83cce91b550779557d081660258c2f7
parent b754af0e37dbd1c9732e6d310e0a0f7402dabbff
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  2 Jun 2020 12:36:48 +0200

Update how <short|long>wave is handled

Add the htsky_spectral_type enumerate

Diffstat:
Msrc/htsky.c | 164++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------
Msrc/htsky.h | 16+++++++++++-----
Msrc/htsky_atmosphere.c | 28+++++++++++++++++++---------
Msrc/htsky_c.h | 4+++-
Msrc/htsky_cloud.c | 52+++++++++++++++++++++++++++++++++++-----------------
Msrc/htsky_dump_cloud_vtk.c | 12++++++++----
6 files changed, 175 insertions(+), 101 deletions(-)

diff --git a/src/htsky.c b/src/htsky.c @@ -56,7 +56,21 @@ check_args(const struct htsky_args* args) && args->grid_max_definition[2] && args->name && args->nthreads - && args->optical_thickness >= 0; + && args->optical_thickness >= 0 + && (unsigned)args->spectral_type < HTSKY_SPECTRAL_TYPES_COUNT__ + && args->wlen_range[0] <= args->wlen_range[1]; +} + +static INLINE const char* +spectral_type_string(const enum htsky_spectral_type type) +{ + const char* str = NULL; + switch(type) { + case HTSKY_SPECTRAL_LW: str = "longwave"; break; + case HTSKY_SPECTRAL_SW: str = "shortwave"; break; + default: FATAL("Unreachable code.\n"); break; + } + return str; } static res_T @@ -73,7 +87,7 @@ setup_bands_properties(struct htsky* sky) 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"); + spectral_type_string(sky->spectral_type)); res = RES_MEM_ERR; goto error; } @@ -82,10 +96,14 @@ setup_bands_properties(struct htsky* sky) double band_wlens[2]; const size_t i = iband - sky->bands_range[0]; - if(sky->is_long_wave) { - HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); - } else { - HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } band_wlens[0] = wavenumber_to_wavelength(band.wave_numbers[1]); band_wlens[1] = wavenumber_to_wavelength(band.wave_numbers[0]); @@ -167,7 +185,7 @@ gas_fetch_raw_property HTGOP(position_to_layer_id(sky->htgop, pos[2], &ilayer)); HTGOP(get_layer(sky->htgop, ilayer, &layer)); - if(sky->is_long_wave) { + if(sky->spectral_type == HTSKY_SPECTRAL_LW) { struct htgop_layer_lw_spectral_interval band; HTGOP(layer_get_lw_spectral_interval(&layer, iband, &band)); @@ -204,8 +222,9 @@ gas_fetch_raw_property } } else { struct htgop_layer_sw_spectral_interval band; - HTGOP(layer_get_sw_spectral_interval(&layer, iband, &band)); + ASSERT(sky->spectral_type == HTSKY_SPECTRAL_SW); + HTGOP(layer_get_sw_spectral_interval(&layer, iband, &band)); if(!in_clouds) { /* Pos is outside the clouds. Directly fetch the nominal optical * properties */ @@ -319,7 +338,7 @@ setup_cache_stream 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->spectral_type, 1); WRITE(sky->bands_range, 2); #undef WRITE CHK(fflush(fp) == 0); @@ -328,7 +347,7 @@ setup_cache_stream struct stat htgop_statbuf2; struct stat htmie_statbuf2; int cache_version; - int is_long_wave; + enum htsky_spectral_type spectral_type; size_t bands_range[2]; /* Read the cache header */ @@ -359,7 +378,7 @@ setup_cache_stream READ(&htgop_statbuf2.st_mtim, 1); READ(&htmie_statbuf2.st_ino, 1); READ(&htmie_statbuf2.st_mtim, 1); - READ(&is_long_wave, 1); + READ(&spectral_type, 1); READ(bands_range, 2); #undef READ @@ -382,7 +401,7 @@ setup_cache_stream /* 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 + if(spectral_type != sky->spectral_type || 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", @@ -422,12 +441,16 @@ print_spectral_info(const struct htsky* sky) 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband_low, &band_low)); + HTGOP(get_lw_spectral_interval(sky->htgop, iband_upp, &band_upp)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband_low, &band_low)); + HTGOP(get_sw_spectral_interval(sky->htgop, iband_upp, &band_upp)); + break; + default: FATAL("Unreachable code.\n"); break; } log_info(sky, "Sky data defined in [%g, %g] nanometers over %lu %s.\n", @@ -441,10 +464,14 @@ print_spectral_info(const struct htsky* sky) 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } naccels += band.quadrature_length; } @@ -486,6 +513,7 @@ htsky_create struct time t0, t1; struct mem_allocator* allocator = NULL; struct htsky* sky = NULL; + double wnums[2]; char buf[128]; int nthreads_max; int force_cache_upd = 0; @@ -516,7 +544,7 @@ htsky_create ref_init(&sky->ref); sky->allocator = allocator; sky->verbose = args->verbose; - sky->is_long_wave = args->is_long_wave ; + sky->spectral_type = args->spectral_type ; sky->repeat_clouds = args->repeat_clouds; sky->is_cloudy = args->htcp_filename != NULL; darray_split_init(sky->allocator, &sky->svx2htcp_z); @@ -565,29 +593,20 @@ htsky_create goto error; } - - /* If the spectral range is degenerated use the default wavelength range, - * i.e. the short wave wavelengths of the CIE XYZ color space */ - if(args->wlen_range[0] > args->wlen_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 { - /* the range is not degenerated ; it is either sw or lw */ - double wnums[2]; - sky->is_long_wave = args->is_long_wave; /* 0 by default */ - wnums[0] = wavelength_to_wavenumber(args->wlen_range[1]); - wnums[1] = wavelength_to_wavenumber(args->wlen_range[0]); - if (sky->is_long_wave) { + /* Retrieve the spectral bands */ + wnums[0] = wavelength_to_wavenumber(args->wlen_range[1]); + wnums[1] = wavelength_to_wavenumber(args->wlen_range[0]); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: res = htgop_get_lw_spectral_intervals(sky->htgop, wnums, sky->bands_range); - } else { + break; + case HTSKY_SPECTRAL_SW: res = htgop_get_sw_spectral_intervals(sky->htgop, wnums, sky->bands_range); - } - if(res != RES_OK) goto error; - } - + break; + default: FATAL("Unreachable code.\n"); break; + } + if(res != RES_OK) goto error; + print_spectral_info(sky); /* Setup the atmopshere */ @@ -949,10 +968,14 @@ htsky_get_spectral_band_quadrature_length ASSERT(sky); 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } return band.quadrature_length; } @@ -967,13 +990,16 @@ htsky_get_spectral_band_bounds res_T res = RES_OK; ASSERT(sky && wavelengths); - 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; + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + res = htgop_get_lw_spectral_interval(sky->htgop, iband, &specint); + break; + case HTSKY_SPECTRAL_SW: + res = htgop_get_sw_spectral_interval(sky->htgop, iband, &specint); + break; + default: FATAL("Unreachable code.\n"); break; } + 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]); @@ -1003,11 +1029,11 @@ htsky_get_raw_spectral_bounds(const struct htsky* sky, double wavelengths[2]) return RES_OK; } -int -htsky_is_long_wave(const struct htsky* htsky) +enum htsky_spectral_type +htsky_get_spectral_type(const struct htsky* htsky) { ASSERT(htsky); - return htsky->is_long_wave; + return htsky->spectral_type; } size_t @@ -1016,10 +1042,14 @@ htsky_find_spectral_band(const struct htsky* sky, const double wavelength) 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(find_lw_spectral_interval_id(sky->htgop, wnum, &iband)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(find_sw_spectral_interval_id(sky->htgop, wnum, &iband)); + break; + default: FATAL("Unreachable code.\n"); break; } return iband; } @@ -1035,10 +1065,14 @@ htsky_spectral_band_sample_quadrature 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } HTGOP(spectral_interval_sample_quadrature(&band, r, &iquad)); return iquad; diff --git a/src/htsky.h b/src/htsky.h @@ -63,13 +63,19 @@ enum htsky_svx_op { HTSKY_SVX_OPS_COUNT__ }; +enum htsky_spectral_type { + HTSKY_SPECTRAL_LW, /* Longwave */ + HTSKY_SPECTRAL_SW, /* Longwave */ + HTSKY_SPECTRAL_TYPES_COUNT__ +}; + struct htsky_args { const char* htcp_filename; const char* htgop_filename; 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 */ - int is_long_wave ; + enum htsky_spectral_type spectral_type; double wlen_range[2]; /* Spectral range to handle. In nm */ unsigned grid_max_definition[3]; /* Maximum definition of the grid */ double optical_thickness; /* Threshold used during octree building */ @@ -84,8 +90,8 @@ struct htsky_args { NULL, /* htmie filename */ \ NULL, /* cache filename */ \ "sky", /* Name */ \ - 0, /* is_long_wave */ \ - {DBL_MAX,-DBL_MAX}, /* Spectral integration. Degenerated <=> short wave */ \ + HTSKY_SPECTRAL_TYPES_COUNT__, /* spectral type */ \ + {DBL_MAX,-DBL_MAX}, /* Spectral integration range */ \ {UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \ 1, /* Optical thickness a*/ \ (unsigned)~0, /* #threads */ \ @@ -215,8 +221,8 @@ htsky_get_raw_spectral_bounds (const struct htsky* sky, double wavelengths[2]); -HTSKY_API int -htsky_is_long_wave +HTSKY_API enum htsky_spectral_type +htsky_get_spectral_type (const struct htsky* sky); /* Return the index of the band containing the submitted wavelength or SIZE_MAX diff --git a/src/htsky_atmosphere.c b/src/htsky_atmosphere.c @@ -58,7 +58,7 @@ 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 ... */ - if(ctx->sky->is_long_wave) { + if(ctx->sky->spectral_type == HTSKY_SPECTRAL_LW) { FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { struct htgop_layer layer; struct htgop_layer_lw_spectral_interval band; @@ -83,6 +83,8 @@ atmosphere_vox_get(const size_t xyz[3], void* dst, void* context) kext_min = ka_min; kext_max = ka_max; } else { + ASSERT(ctx->sky->spectral_type == HTSKY_SPECTRAL_SW); + FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { struct htgop_layer layer; struct htgop_layer_sw_spectral_interval band; @@ -211,10 +213,14 @@ atmosphere_setup(struct htsky* sky, const double optical_thickness_threshold) struct htgop_spectral_interval band; ctx.iband = i + sky->bands_range[0]; - 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, ctx.iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, ctx.iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } sky->atmosphere[i] = MEM_CALLOC(sky->allocator, @@ -273,10 +279,14 @@ atmosphere_clean(struct htsky* sky) 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } FOR_EACH(iquad, 0, band.quadrature_length) { diff --git a/src/htsky_c.h b/src/htsky_c.h @@ -17,6 +17,8 @@ #ifndef HTSKY_C_H #define HTSKY_C_H +#include "htsky.h" + #include <high_tune/htcp.h> #include <rsys/logger.h> @@ -81,7 +83,7 @@ struct htsky { /* Ids and optical properties of spectral bands loaded by HTGOP and currently * handled by the sky */ - int is_long_wave; + enum htsky_spectral_type spectral_type; size_t bands_range[2]; /* Inclusive band ids */ struct band_prop* bands; diff --git a/src/htsky_cloud.c b/src/htsky_cloud.c @@ -523,7 +523,7 @@ cloud_vox_get_gas ka[1] = ks[1] = kext[1] =-DBL_MAX; /* For each atmospheric layer that overlaps the SVX voxel ... */ - if(ctx->sky->is_long_wave) { + if(ctx->sky->spectral_type == HTSKY_SPECTRAL_LW) { FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { struct htgop_layer_lw_spectral_interval band; double k[2]; @@ -546,6 +546,8 @@ cloud_vox_get_gas kext[0] = ka[0]; kext[1] = ka[1]; } else { + ASSERT(ctx->sky->spectral_type == HTSKY_SPECTRAL_SW); + FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) { struct htgop_layer_sw_spectral_interval band; double k[2]; @@ -754,10 +756,14 @@ cloud_write_octrees(const struct htsky* sky, FILE* stream) size_t iquad; 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } ASSERT(sky->clouds[i]); @@ -793,10 +799,14 @@ cloud_read_octrees(struct htsky* sky, const char* stream_name, FILE* stream) size_t iquad; 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } ASSERT(sky->clouds[i]); @@ -896,10 +906,14 @@ cloud_setup const size_t iband = i + sky->bands_range[0]; size_t iquad; - if(sky->is_long_wave) { - HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); - } else { - HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } sky->clouds[i] = MEM_CALLOC(sky->allocator, @@ -990,10 +1004,14 @@ cloud_clean(struct htsky* sky) size_t iquad; 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)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &band)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band)); + break; + default: FATAL("Unreachable code.\n"); break; } if(!sky->clouds[i]) continue; diff --git a/src/htsky_dump_cloud_vtk.c b/src/htsky_dump_cloud_vtk.c @@ -204,10 +204,14 @@ htsky_dump_cloud_vtk ASSERT(ncells == cloud->octree_desc.nleaves); /* Fetch the spectral interval descriptor */ - if(sky->is_long_wave) { - HTGOP(get_lw_spectral_interval(sky->htgop, iband, &specint)); - } else { - HTGOP(get_sw_spectral_interval(sky->htgop, iband, &specint)); + switch(sky->spectral_type) { + case HTSKY_SPECTRAL_LW: + HTGOP(get_lw_spectral_interval(sky->htgop, iband, &specint)); + break; + case HTSKY_SPECTRAL_SW: + HTGOP(get_sw_spectral_interval(sky->htgop, iband, &specint)); + break; + default: FATAL("Unreachable code.\n"); break; } /* Write headers */