solstice

Compute collected power and efficiencies of a solar plant
git clone git://git.meso-star.com/solstice.git
Log | Files | Refs | README | LICENSE

commit d7a8abb0333ea295235187d45538d0b18013b0dc
parent 45c734e0660122a604558c82e5ae26a73c336b8e
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Mon, 22 Jun 2026 11:20:08 +0200

Allow to define sun direction through location + time

Diffstat:
MMakefile.core | 31++++++++++++++++++++++++++-----
Mdoc/solstice-output.5 | 26+++++++++++++++++---------
Mdoc/solstice.1.in | 99+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------
Msrc/solstice.c | 102+++++++++++++++++++++++++++++++++++++++++++++++++------------------------------
Msrc/solstice.h | 32+++++++++++++++++++++++++++++---
Msrc/solstice_args.c | 172++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------
Msrc/solstice_args.h.in | 30+++++++++++++++++++++++++++---
Msrc/solstice_sun.c | 8+++++++-
Msrc/test_solstice_args.c | 195+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------
Asrc/test_solstice_simulation_time.c | 613+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
10 files changed, 1145 insertions(+), 163 deletions(-)

diff --git a/Makefile.core b/Makefile.core @@ -117,7 +117,8 @@ lint: ################################################################################ TEST_SRC =\ src/test_solstice_args.c\ - src/test_solstice_simulation.c + src/test_solstice_simulation.c\ + src/test_solstice_simulation_time.c TEST_OBJ =\ $(TEST_SRC:.c=.o) TEST_DEP =\ @@ -141,7 +142,7 @@ src/score-local.pc: score.pc.in src/.config_score_test # Regular cflags PKG_CONFIG_LOCAL = PKG_CONFIG_PATH="./src:$${PKG_CONFIG_PATH}" $(PKG_CONFIG) INCS_TEST = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys score-local sprs-local srcv-local) -LIBS_TEST = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys score-local sprs-local srcv-local)\ +LIBS_TEST = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys score-local sprs-local srcv-local scem)\ -lm tests: library src/score-local.pc $(TEST_DEP) $(TEST_TGT) @@ -186,9 +187,7 @@ test: tests \ for i in $(TEST_SRC); do \ test="$$(basename "$${i}" ".c")"; \ - if [ "$${test}" != "test_solstice_simulation" ]; then \ - check "$${test}" "$${test}"; \ - else \ + if [ "$${test}" = "test_solstice_simulation" ]; then \ check test_solstice_simulation_beam_down \ test_solstice_simulation ./solstice yaml/ beam_down; \ check test_solstice_simulation_test01 \ @@ -207,6 +206,27 @@ test: tests test_solstice_simulation ./solstice yaml/ test07; \ check test_solstice_simulation_test08 \ test_solstice_simulation ./solstice yaml/ test08; \ + elif [ "$${test}" = "test_solstice_simulation_time" ]; then \ + check test_solstice_simulation_time_beam_down \ + test_solstice_simulation_time ./solstice yaml/ beam_down; \ + check test_solstice_simulation_time_test01 \ + test_solstice_simulation_time ./solstice yaml/ test01; \ + check test_solstice_simulation_time_test02 \ + test_solstice_simulation_time ./solstice yaml/ test02; \ + check test_solstice_simulation_time_test03 \ + test_solstice_simulation_time ./solstice yaml/ test03; \ + check test_solstice_simulation_time_test04 \ + test_solstice_simulation_time ./solstice yaml/ test04; \ + check test_solstice_simulation_time_test05 \ + test_solstice_simulation_time ./solstice yaml/ test05; \ + check test_solstice_simulation_time_test06 \ + test_solstice_simulation_time ./solstice yaml/ test06; \ + check test_solstice_simulation_time_test07 \ + test_solstice_simulation_time ./solstice yaml/ test07; \ + check test_solstice_simulation_time_test08 \ + test_solstice_simulation_time ./solstice yaml/ test08; \ + else \ + check "$${test}" "$${test}"; \ fi \ done; \ \ @@ -226,5 +246,6 @@ $(TEST_OBJ) : src/.config_score_test src/score-local.pc test_solstice_args\ test_solstice_simulation\ +test_solstice_simulation_time\ : src/.config_score_test src/score-local.pc $(LIBNAME) $(CC) $(CFLAGS_TEST) -o $@ src/$@.o $(LDFLAGS_TEST) diff --git a/doc/solstice-output.5 b/doc/solstice-output.5 @@ -69,7 +69,7 @@ are on a single line until a closing quote mark. | <rendering-output> # -r option <simulation-output> - ::= <sun-direction> + ::= <sun-specification> <counts> <global> [ <receivers-list> ] @@ -79,22 +79,27 @@ are on a single line until a closing quote mark. [ <simulation-output> ... ] <dump-geometry-output> - ::= <sun-direction> + ::= <sun-specification> <geometry-data> [ <dump-geometry-output> ... ] <dump-radiative-paths-output> - ::= <sun-direction> + ::= <sun-specification> VTK-RADIATIVE-PATHS [ <dump-radiative-paths-output> ... ] <rendering-output> - ::= <sun-direction> + ::= <sun-specification> PPM-FILE # ASCII PPM with 8-bits per component [1] [ <rendering-output> ... ] -.Ed -.Bd -literal -<sun-direction> ::= "#--- Sun direction: <alpha> <beta> (<sun-vector>)" + +<sun-specification> + ::= <sun-direction> + | <sun-time> + +<sun-direction> ::= "#--- Sun direction: <azimuth> <elevation> (<sun-vector>)" + +<sun-time> ::= "#--- Sun location and time: <lat> <long> <time> (<sun-vector>)" <counts> ::= "<#globals> <#receivers> <#primaries> <#samples> <#failed>" @@ -169,9 +174,12 @@ are on a single line until a closing quote mark. <area> ::= REAL # in ]0, INF) <real3> ::= REAL REAL REAL -<alpha> ::= REAL # Degrees in [0, 360[ -<beta> ::= REAL # Degrees in [0, 90] +<azimuth> ::= REAL # Degrees in [0, 360[ +<elevation> ::= REAL # Degrees in [0, 90] <sun-vector> ::= <real3> +<lat> ::= REAL # Degrees in [-90, +90] +<long> ::= REAL # Degrees in [-180, +180] +<time> ::= STRING # Same YYYY-MM-DDThh:mm:ss format as on the CLI <incoming-flux> ::= <estimate> <in-if-no-mat-loss> ::= <estimate> diff --git a/doc/solstice.1.in b/doc/solstice.1.in @@ -21,7 +21,6 @@ .Nd compute the power collected by a concentrated solar plant .Sh SYNOPSIS .Nm -.Nm .Op Ar option ... .Op Ar file .Nm @@ -94,7 +93,7 @@ geometries into CAD files. The .Fl p option saves the sampled radiative paths used by the estimates, allowing -to visualise them externally, which may be a great help to identify a +the user to visualise them externally, which may be a great help to identify a design issue. Finally, the .Fl r @@ -104,26 +103,21 @@ they replace the default .Nm behaviour. .Pp -Any coordinate-related question in +All coordinates in .Nm -must be considered with the right-handed convention in mind. +follow the right-handed convention. .Sh OPTIONS .Bl -tag -width Ds -.It Fl D Ar alpha,beta Ns Op : Ns Ar ... -List of sun directions. -A direction is defined by two angles in degrees. -The first one, -.Ar alpha , -is an azimuthal angle in [0,\ 360[ and the second one, -.Ar beta , -is an elevation in [0,\ 90]. +.It Fl D Ar azimuth,elevation Ns +A sun direction, defined by two angles in degrees. +The first one is the azimuthal angle in [0,\ 360[, and the second one +is the elevation in [0,\ 90]. Each provided sun direction triggers a new computation whose results are concatenated to the -.Ar output -file. +.Ar output . .Pp Following the right-handed convention, azimuthal rotation is -counter-clockwise, with 0\(de on the X axis. +counter-clockwise, with 0\(de on the +X axis. Elevation starts from 0\(de for directions in the XY plane, up to 90\(de at zenith. Thus @@ -134,9 +128,9 @@ and .Fl D Ns Ar 270,0 produce solar vectors {-1,0,0}, {0,-1,0}, {+1,0,0} and {0,+1,0} respectively, while -.Fl D Ns Ar alpha , Ns 90 +.Fl D Ns Ar azimuth , Ns 90 produces {0,0,-1} regardless of the value of -.Ar alpha . +.Ar azimuth . .It Fl f Force overwrite of the output files, i.e.\& the .Ar output @@ -196,10 +190,22 @@ By default .Ar samples-count is set to .Sy @SOLSTICE_ARGS_DEFAULT_NREALISATIONS@ . +.It Fl L Ar latitude,longitude +Define the location of the solar plant. +The +.Ar latitude +must be in [-90, 90] degrees relative to the equator, counting +positive towards the north. +The +.Ar longitude +must be in [-180, 180] degrees relative to Greenwich, counting +positive towards the east. +A location must be defined before any time is defined. +It is then applied to any following time until a new location is defined. .It Fl o Ar output Write results to .Ar output -with respect to the +in the .Xr solstice-output 5 format. If not defined, write results to standard output. @@ -263,7 +269,7 @@ the sun is ignored and the only light source is positioned at the camera position. In .Cm pt -mode, the scene is rendered with the un-biased path-tracing Monte-Carlo +mode, the scene is rendered with the unbiased path-tracing Monte-Carlo algorithm; the materials described in the committed .Ar file as well as the submitted sun directions are correctly handled and a @@ -303,8 +309,27 @@ By default, is set to .Sy { Ns @SOLSTICE_ARGS_DEFAULT_CAMERA_UP@ Ns } . .El +.It Fl T Ar utc_time Ns +A time, defined by a string in Coordinated Universal Time +.Pq UTC+00:00 . +.Pp +The date must be in the format +.Dq YYYY-MM-DDThh:mm:ss , +i.e., as printed for the current date by the following date command: +.Bd -literal -offset indent +date -u +"%Y-%m-%dT%H:%M:%S" +.Ed +.Pp +Each provided +.Ar time +is used in conjunction with the last provided location to compute a +sun direction using the PSA algorithm +.Pq see Sy Blanco et al. in the SEE ALSO section +and triggers a new computation whose results are +concatenated to the +.Ar output . .It Fl t Ar threads-count -Hint on the number of threads to use. +Hint at the number of threads to use. By default, as many threads as CPU cores are used. .It Fl v Make @@ -326,16 +351,28 @@ are written to .Pa output even though this file already exists: .Bd -literal -offset indent -solstice -D45,70:50,75 -R rcvs.yaml -n 10000 -f -o output input.yaml +solstice -D45,70 -D50,75 -R rcvs.yaml -n 10000 -f -o output input.yaml +.Ed +.Pp +Launch one simulation at Toulouse, France, at 2 PM on May 1 2022 +.Pq UTC+2 . +The solar facility is described in +.Pa input.yaml +and the receivers on which the integrations must be performed are declared +in +.Pa rcvs.yaml . +The results are written to standard output: +.Bd -literal -offset indent +solstice -L43.605,1.445 -T"2022-05-01T12:00:00" -R rcvs.yaml input.yaml .Ed .Pp Generate a mesh for each geometry described in .Pa input.yaml and save them in .Pa output -with respect to the Alias Wavefront OBJ format. +in the Alias Wavefront OBJ format. The meshes are positioned according to their orientation constraints, -with respect to the sun direction whose azimuthal and elevation angles +for the sun direction whose azimuthal and elevation angles are {30,60}. Use .Xr csplit 1 @@ -359,7 +396,7 @@ csplit -f geom -b %02d.obj -z --suppress-matched output /^---$/ {*} .Pp Trace 100 radiative paths into the solar plant described in .Pa input.yaml , -with respect to the sun direction whose azimuthal and elevation angles +for the sun direction whose azimuthal and elevation angles are 0 and 90 degrees, respectively. Write the .Xr solstice-output 5 @@ -375,7 +412,7 @@ solstice -n 100 -D0,90 -R rcvs.yaml -p default input.yaml | sed '1d' > paths.vtk .Pp Use the path-tracing rendering algorithm to draw the solar plant .Pa solplant.yaml -with respect to the sun direction whose azimuthal and elevation angles +for the sun direction whose azimuthal and elevation angles are 180 and 45 degrees, respectively. Use 64 samples per pixel to estimate the per-pixel radiance and fix the camera up vector to {0,0,1}. @@ -397,6 +434,18 @@ solstice -D180,45 -r up=0,0,1:rmode=pt:spp=64 solplant.yaml | sed '1d' | feh - .Xr solstice-input 5 , .Xr solstice-output 5 , .Xr solstice-receiver 5 +.Rs +.%A Manuel Blanco-Muriel +.%A Diego C. Alarcón-Padilla +.%A Teodoro López-Moratalla +.%A Martín Lara-Coira +.%T Computing the solar vector +.%J Solar Energy +.%V 70 +.%N 5 +.%D 2001 +.%P 431-441 +.Re .Sh HISTORY .Nm was initially developed with funding from the diff --git a/src/solstice.c b/src/solstice.c @@ -151,7 +151,7 @@ auto_look_at double radius; double depth; res_T res; - ASSERT(scn && fov_x && proj_ratio && up); + ASSERT(scn && fov_x > 0 && proj_ratio > 0 && up); res = ssol_scene_compute_aabb(scn, flower, fupper); if(res != RES_OK) { @@ -316,9 +316,11 @@ error: goto exit; } +/* This function differs from scem_sun_position_to_sun_vector because 1) angles + * are in degrees VS radians, 2) the conventions on azimuth are different. */ static INLINE void spherical_to_cartesian_sun_dir - (struct solstice_args_spherical* spherical, double sun_dir[3]) + (const struct solstice_spherical* spherical, double sun_dir[3]) { double cos_azimuth; double sin_azimuth; @@ -344,38 +346,41 @@ spherical_to_cartesian_sun_dir static res_T setup_sun_dirs(struct solstice* solstice, const struct solstice_args* args) { - double* sun_dirs = NULL; - double* sun_angles = NULL; + struct solstice_sun_dir* sun_dirs = NULL; size_t i; res_T res = RES_OK; ASSERT(solstice && args); - res = darray_double_resize(&solstice->sun_dirs, args->nsun_dirs*3/*#dims*/); + res = darray_sun_dir_resize(&solstice->sun_dirs, args->nsun_dirs); if(res != RES_OK) { fprintf(stderr, "Could not reserve the list of %lu sun directions.\n", (unsigned long)args->nsun_dirs); goto error; } - res = darray_double_resize(&solstice->sun_angles, args->nsun_dirs*2/*#dims*/); - if(res != RES_OK) { - fprintf(stderr, - "Could not reserve the list of %lu sun angles.\n", - (unsigned long)args->nsun_dirs); - goto error; - } - sun_dirs = darray_double_data_get(&solstice->sun_dirs); - sun_angles = darray_double_data_get(&solstice->sun_angles); + sun_dirs = darray_sun_dir_data_get(&solstice->sun_dirs); FOR_EACH(i, 0, args->nsun_dirs) { - spherical_to_cartesian_sun_dir(args->sun_dirs + i, sun_dirs + i*3/*#dims*/); - d2(sun_angles + i*2, args->sun_dirs[i].azimuth, args->sun_dirs[i].elevation); + struct solstice_args_sun_dir* r = args->sun_dirs + i; + switch(r->type) { + case SOLSTICE_ARGS_SPHERICAL: + sun_dirs[i].type = SOLSTICE_SUN_DIR_SPHERICAL; + sun_dirs[i].u.spherical.azimuth = r->spherical.azimuth; + sun_dirs[i].u.spherical.elevation = r->spherical.elevation; + break; + case SOLSTICE_SUN_DIR_LOCATION_TIME: + sun_dirs[i].type = SOLSTICE_SUN_DIR_LOCATION_TIME; + sun_dirs[i].u.location_time.time = r->location_time.time; + sun_dirs[i].u.location_time.latitude = r->location_time.latitude; + sun_dirs[i].u.location_time.longitude = r->location_time.longitude; + break; + default: FATAL("Invalid enum value.\n"); + } } exit: return res; error: - darray_double_clear(&solstice->sun_dirs); - darray_double_clear(&solstice->sun_angles); + darray_sun_dir_clear(&solstice->sun_dirs); goto exit; } @@ -601,8 +606,7 @@ solstice_init darray_nodes_init(allocator, &solstice->roots); darray_nodes_init(allocator, &solstice->pivots); darray_receiver_init(allocator, &solstice->rcvs_list); - darray_double_init(allocator, &solstice->sun_dirs); - darray_double_init(allocator, &solstice->sun_angles); + darray_sun_dir_init(allocator, &solstice->sun_dirs); solstice->allocator = allocator ? allocator : &mem_default_allocator; @@ -741,16 +745,14 @@ solstice_release(struct solstice* solstice) darray_nodes_release(&solstice->roots); darray_nodes_release(&solstice->pivots); darray_receiver_release(&solstice->rcvs_list); - darray_double_release(&solstice->sun_dirs); - darray_double_release(&solstice->sun_angles); + darray_sun_dir_release(&solstice->sun_dirs); logger_release(&solstice->logger); } res_T solstice_run(struct solstice* solstice) { - const double* sun_dirs = NULL; - const double* sun_angles = NULL; + const struct solstice_sun_dir* sun_dirs = NULL; size_t nsun_dirs = 0; size_t i; int dump; @@ -758,11 +760,8 @@ solstice_run(struct solstice* solstice) res_T res = RES_OK; ASSERT(solstice); - sun_dirs = darray_double_cdata_get(&solstice->sun_dirs); - sun_angles = darray_double_cdata_get(&solstice->sun_angles); - nsun_dirs = darray_double_size_get(&solstice->sun_dirs); - ASSERT(nsun_dirs%3 == 0); - nsun_dirs /= 3/*#dims*/; + sun_dirs = darray_sun_dir_cdata_get(&solstice->sun_dirs); + nsun_dirs = darray_sun_dir_size_get(&solstice->sun_dirs); dump = solstice->dump_format != SOLSTICE_ARGS_DUMP_NONE; draw = solstice->framebuffer != NULL; @@ -788,20 +787,45 @@ solstice_run(struct solstice* solstice) } } else { FOR_EACH(i, 0, nsun_dirs) { - const double* sun_dir = sun_dirs + i*3/*#dims*/; - const double* sun_angle = sun_angles + i*2/*#angles*/; - fprintf(solstice->output, "#--- Sun direction: %g %g (%g %g %g)\n", - SPLIT2(sun_angle), SPLIT3(sun_dir)); + double sun_dir[3]; + if(sun_dirs[i].type == SOLSTICE_SUN_DIR_SPHERICAL) { + spherical_to_cartesian_sun_dir(&sun_dirs[i].u.spherical, sun_dir); + fprintf(solstice->output, "#--- Sun direction: %g %g (%g %g %g)\n", + sun_dirs[i].u.spherical.azimuth, sun_dirs[i].u.spherical.elevation, + SPLIT3(sun_dir)); + res = ssol_sun_set_direction(solstice->sun, sun_dir); + if(res != RES_OK) { + fprintf(stderr, "Could not update the sun direction.\n"); + goto error; + } + } else { + char buf[128]; + struct ssol_location loc; + const struct solstice_location_time* lt = &sun_dirs[i].u.location_time; + ASSERT(sun_dirs[i].type == SOLSTICE_SUN_DIR_LOCATION_TIME); + loc.latitude = lt->latitude; + loc.longitude = lt->longitude; + res = ssol_sun_set_location_and_date(solstice->sun, &loc, &lt->time); + if(res != RES_OK) { + fprintf(stderr, "Could not update the sun location and time.\n"); + goto error; + } + res = ssol_sun_get_direction(solstice->sun, sun_dir); + if(res != RES_OK) goto error; + if(0 == strftime(buf, sizeof(buf), + "%Y-%m-%dT%H:%M:%S", &sun_dirs[i].u.location_time.time)) { + res = RES_BAD_ARG; + goto error; + } + fprintf(solstice->output, + "#--- Sun location and time: %g %g %s (%g %g %g)\n", + lt->latitude, lt->longitude, buf, + SPLIT3(sun_dir)); + } res = solstice_update_entities(solstice, sun_dir); if(res != RES_OK) goto error; - res = ssol_sun_set_direction(solstice->sun, sun_dir); - if(res != RES_OK) { - fprintf(stderr, "Could not update the sun direction.\n"); - goto error; - } - if(draw) { res = solstice_draw(solstice); if(res != RES_OK) goto error; diff --git a/src/solstice.h b/src/solstice.h @@ -17,7 +17,6 @@ #ifndef SOLSTICE_H #define SOLSTICE_H -#include "parser/solparser_material.h" #include "receivers/srcvl.h" #include "solstice_args.h" @@ -161,6 +160,33 @@ struct solstice_primary { #define HTABLE_DATA struct solstice_primary #include <rsys/hash_table.h> +enum solstice_sun_dir_type { + SOLSTICE_SUN_DIR_SPHERICAL, + SOLSTICE_SUN_DIR_LOCATION_TIME +}; + +struct solstice_location_time { + struct tm time; + double latitude, longitude; +}; + +struct solstice_spherical { + double azimuth; /* In radians */ + double elevation; /* In radians */ +}; + +struct solstice_sun_dir { + enum solstice_sun_dir_type type; + union { + struct solstice_spherical spherical; + struct solstice_location_time location_time; + } u; +}; + +#define DARRAY_NAME sun_dir +#define DARRAY_DATA struct solstice_sun_dir +#include <rsys/dynamic_array.h> + struct solstice { struct ssol_device* ssol; struct ssol_scene* scene; @@ -195,8 +221,8 @@ struct solstice { /* Dump radiative paths mode */ struct ssol_path_tracker path_tracker; - struct darray_double sun_dirs; /* List of double3 */ - struct darray_double sun_angles; + /* Sun directions, either a spherical dir or by location+time */ + struct darray_sun_dir sun_dirs; size_t nexperiments; /* # MC experiments */ FILE* output; /* Output stream */ diff --git a/src/solstice_args.c b/src/solstice_args.c @@ -14,12 +14,14 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +#define _XOPEN_SOURCE /* strptime support */ #define _POSIX_C_SOURCE 2 #include "solstice_args.h" #include "solstice_version.h" #include <rsys/cstr.h> +#include <rsys/double2.h> #include <rsys/double3.h> #include <rsys/stretchy_array.h> @@ -31,6 +33,7 @@ #endif #include <string.h> +#include <time.h> /******************************************************************************* * Helper functions @@ -45,7 +48,7 @@ print_help(const char* program) "man page for more informations.\n\n", program); printf( -" -D <dirs> list of sun directions.\n"); +" -D <dir> request a computation for a sun direction.\n"); printf( " -f overwrite the output files if they already exist, i.e. the\n" " OUTPUT file and the output RNG state.\n"); @@ -59,6 +62,8 @@ print_help(const char* program) " -n SAMPLES number of Monte Carlo samples. Default is %lu.\n", SOLSTICE_ARGS_DEFAULT.nexperiments); printf( +" -L <location> set the location of the plant on earth surface.\n"); + printf( " -o OUTPUT write results to OUTPUT. If not defined, write results to\n" " standard output.\n"); printf( @@ -70,6 +75,9 @@ print_help(const char* program) printf( " -r <rendering> switch in rendering mode and configure it.\n"); printf( +" -T <time> request a computation at a given time at the last specified\n" +" location.\n"); + printf( " -t THREADS hint on the number of threads to use. By default use as\n" " many threads as CPU cores.\n"); printf( @@ -139,58 +147,116 @@ error: } static res_T -parse_sun_dir_list(const char* str, struct solstice_args* args) +parse_sun_spherical(const char* str, struct solstice_args* args) { - char buf[512]; - char* tk; - char* ctx; size_t len; + struct solstice_args_sun_dir sun_dir; + double tmp[2]; res_T res = RES_OK; ASSERT(str && args); - if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) { + res = cstr_to_list_double(str, ',', tmp, &len, 2); + if(res != RES_OK || len != 2) { + if(res == RES_OK) res = RES_BAD_ARG; + fprintf(stderr, "Invalid sun direction `%s'.\n", str); + goto error; + } + + if(tmp[0] < 0 || tmp[0] >= 360) { fprintf(stderr, - "Could not duplicate the list of sun directions `%s'.\n", str); - res = RES_MEM_ERR; + "Invalid azimuth angle `%g'. Azimuth must be in [0, 360[ degrees.\n", + tmp[0]); + res = RES_BAD_ARG; + goto error; + } + if(tmp[1] < 0 || tmp[1] > 90) { + fprintf(stderr, + "Invalid elevation angle `%g'. Elevation must be in [0, 90] degrees.\n", + tmp[1]); + res = RES_BAD_ARG; goto error; } - strncpy(buf, str, sizeof(buf)); - tk = strtok_r(buf, ":", &ctx); - while(tk) { - struct solstice_args_spherical spherical; - double tmp[2]; + sun_dir.type = SOLSTICE_ARGS_SPHERICAL; + sun_dir.spherical.azimuth = tmp[0]; + sun_dir.spherical.elevation = tmp[1]; + sa_push(args->sun_dirs, sun_dir); + args->nsun_dirs++; - res = cstr_to_list_double(tk, ',', tmp, &len, 2); - if(res == RES_OK && len != 2) res = RES_BAD_ARG; - if(res != RES_OK) { - fprintf(stderr, "Invalid sun direction `%s'.\n", tk); - goto error; - } +exit: + return res; +error: + if(args->sun_dirs) { + sa_release(args->sun_dirs); + args->sun_dirs = NULL; + args->nsun_dirs = 0; + } + goto exit; +} - if(tmp[0] < 0 || tmp[0] >= 360) { - fprintf(stderr, - "Invalid azimuth angle `%g'. Azimuth must be in [0, 360[ degrees.\n", +static res_T +parse_location + (const char* str, + double current_loc[2], + struct solstice_args* args) +{ + size_t len; + double tmp[2]; + res_T res = RES_OK; + ASSERT(str && current_loc && args); + + res = cstr_to_list_double(str, ',', tmp, &len, 2); + if(res != RES_OK || len != 2) { + if(res == RES_OK) res = RES_BAD_ARG; + fprintf(stderr, "Invalid location: `%s'.\n", str); + goto error; + } + + if(tmp[0] < -90 || tmp[0] > 90) { + fprintf(stderr, + "Invalid latitude `%g'. Must be in [-90, 90] degrees.\n", tmp[0]); - res = RES_BAD_ARG; - goto error; - } - if(tmp[1] < 0 || tmp[1] > 90) { - fprintf(stderr, - "Invalid elevation angle `%g'. Elevation must be in [0, 90] degrees.\n", + res = RES_BAD_ARG; + goto error; + } + if(tmp[1] < -180 || tmp[1] > 180) { + fprintf(stderr, + "Invalid longitude `%g'. Must be in [-180, 180] degrees.\n", tmp[1]); - res = RES_BAD_ARG; - goto error; - } + res = RES_BAD_ARG; + goto error; + } - spherical.azimuth = tmp[0]; - spherical.elevation = tmp[1]; - sa_push(args->sun_dirs, spherical); + d2_set(current_loc, tmp); +exit: + return res; +error: + goto exit; +} - tk = strtok_r(NULL, ":", &ctx); +static res_T +parse_time + (const char* str, + const double current_loc[2], + struct solstice_args* args) +{ + res_T res = RES_OK; + struct solstice_args_sun_dir sun_dir; + char* p; + ASSERT(str && current_loc && args); + + p = strptime(str, "%Y-%m-%dT%H:%M:%S", &sun_dir.location_time.time); + if(p == NULL || *p != '\0') { + res = RES_BAD_ARG; + fprintf(stderr, "Invalid time `%s'.\n", str); + goto error; } - args->nsun_dirs += sa_size(args->sun_dirs); + sun_dir.type = SOLSTICE_ARGS_LOCATION_TIME; + sun_dir.location_time.latitude = current_loc[0]; + sun_dir.location_time.longitude = current_loc[1]; + sa_push(args->sun_dirs, sun_dir); + args->nsun_dirs++; exit: return res; @@ -545,7 +611,8 @@ res_T solstice_args_init(struct solstice_args* args, const int argc, char** argv) { int opt; - int i; + int i, L_defined = 0; + double current_loc[2]; res_T res = RES_OK; ASSERT(args && argc && argv); @@ -563,10 +630,10 @@ solstice_args_init(struct solstice_args* args, const int argc, char** argv) } optind = 0; - while((opt = getopt(argc, argv, "D:fG:g:hn:o:p:qR:r:t:v")) != -1) { + while((opt = getopt(argc, argv, "D:fG:g:hn:L:o:p:qR:r:T:t:v")) != -1) { switch(opt) { - case 'D': /* Sun directions */ - res = parse_sun_dir_list(optarg, args); + case 'D': /* 1 sun direction */ + res = parse_sun_spherical(optarg, args); break; case 'f': args->force_overwriting = 1; break; case 'h': /* Print short help and exit */ @@ -574,7 +641,7 @@ solstice_args_init(struct solstice_args* args, const int argc, char** argv) solstice_args_release(args); args->quit = 1; goto exit; - case 'n': /* Define the number of experiments */ + case 'n': /* Define the number of MC samples */ res = cstr_to_ulong(optarg, &args->nexperiments); if(res == RES_OK && !args->nexperiments) res = RES_BAD_ARG; break; @@ -589,6 +656,10 @@ solstice_args_init(struct solstice_args* args, const int argc, char** argv) res = RES_BAD_ARG; } break; + case 'L': /* Plant location */ + L_defined = 1; + res = parse_location(optarg, current_loc, args); + break; case 'o': args->output_filename = optarg; break; case 'p': /* Switch in dump radiative paths mode and configure it */ args->dump_paths = 1; @@ -600,6 +671,16 @@ solstice_args_init(struct solstice_args* args, const int argc, char** argv) args->rendering = 1; res = parse_multiple_options(optarg, args, parse_rendering_option); break; + case 'T': /* 1 time */ + if(!L_defined) { + fprintf(stderr, + "%s: cannot use the -%c option with no current location defined '%s'\n", + argv[0], opt, optarg); + res = RES_BAD_ARG; + break; + } + res = parse_time(optarg, current_loc, args); + break; case 't': /* Submit an hint on the number of threads to use */ res = cstr_to_uint(optarg, &args->nthreads); if(res == RES_OK && !args->nthreads) res = RES_BAD_ARG; @@ -610,16 +691,17 @@ solstice_args_init(struct solstice_args* args, const int argc, char** argv) if(res != RES_OK) { if(optarg) { fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", - argv[0], optarg, opt); + argv[0], optarg, opt); } goto error; } } if(!args->rendering - && args->dump_format == SOLSTICE_ARGS_DUMP_NONE - && !args->nsun_dirs) { - fprintf(stderr, "Missing sun direction.\n"); + && args->dump_format == SOLSTICE_ARGS_DUMP_NONE + && args->nsun_dirs == 0) + { + fprintf(stderr, "No computation request.\n"); res = RES_BAD_ARG; goto error; } diff --git a/src/solstice_args.h.in b/src/solstice_args.h.in @@ -21,10 +21,34 @@ #include <solstice/ssol.h> #include <rsys/math.h> +#include <time.h> + +struct solstice_args_location { + double longiture; + double latitude; +}; + +#define SOLSTICE_ARGS_LOCATION_NONE { 0, 0 } + +enum solstice_args_sun_dir_type { + SOLSTICE_ARGS_SPHERICAL, + SOLSTICE_ARGS_LOCATION_TIME +}; struct solstice_args_spherical { - double azimuth; /* In radians */ - double elevation; /* In radians */ + double azimuth; + double elevation; +}; + +struct solstice_args_location_time { + struct tm time; + double latitude, longitude; +}; + +struct solstice_args_sun_dir { + enum solstice_args_sun_dir_type type; + struct solstice_args_location_time location_time; + struct solstice_args_spherical spherical; }; enum solstice_args_dump_format { @@ -51,7 +75,7 @@ struct solstice_args { unsigned nthreads; /* #threads */ /* List of sun directions */ - struct solstice_args_spherical* sun_dirs; + struct solstice_args_sun_dir* sun_dirs; size_t nsun_dirs; struct { diff --git a/src/solstice_sun.c b/src/solstice_sun.c @@ -43,7 +43,7 @@ create_sun_buie res = ssol_sun_set_buie_param(sun, solparser_sun->radang_distrib.buie.csr); if(res != RES_OK) { - fprintf(stderr, "Could setup the buie parameter of the solver sun.\n"); + fprintf(stderr, "Could not setup the buie parameter of the solver sun.\n"); goto error; } @@ -245,6 +245,12 @@ solstice_create_sun(struct solstice* solstice) } if(res != RES_OK) goto error; + res = ssol_sun_set_algorithm(sun, SSOL_SUN_DIRECTION_ALGORITHM_PSA); + if(res != RES_OK) { + fprintf(stderr, "Could not setup the sun direction algorithm.\n"); + goto error; + } + if(!SOLPARSER_ID_IS_VALID(solparser_sun->spectrum)) { res = create_default_sun_spectrum(solstice, solparser_sun, &spectrum); if(res != RES_OK) goto error; diff --git a/src/test_solstice_args.c b/src/test_solstice_args.c @@ -14,6 +14,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +#define _XOPEN_SOURCE /* strptime support */ + #include "solstice_args.h" #include "test_solstice_utils.h" @@ -23,6 +25,7 @@ #include <stdarg.h> #include <string.h> #include <limits.h> +#include <time.h> #ifdef COMPILER_CL #pragma warning(push) @@ -235,79 +238,205 @@ test_rendering(void) cmd_delete(cmd); } +static int +tm_equal(const struct tm* tm1, const struct tm* tm2) { + if(tm1->tm_year != tm2->tm_year) return 0; + if(tm1->tm_mon != tm2->tm_mon) return 0; + if(tm1->tm_mday != tm2->tm_mday) return 0; + if(tm1->tm_hour != tm2->tm_hour) return 0; + if(tm1->tm_min != tm2->tm_min) return 0; + if(tm1->tm_sec != tm2->tm_sec) return 0; + return 1; +} + static void test_sun_dirs(void) { + struct tm ref; struct solstice_args args = SOLSTICE_ARGS_NULL; + char *p; + char timeOK[] = "2020-10-20T10:10:10"; + char* dirsBAD[] = { + "1.2,3.4,5", + "1.2,3.4:1", + "1.2,3.4:5.2,A", + "-0.1,2", + "360,2", + "0,-1", + "0,91", + NULL + }; + char* locsBAD[] = { + "-100,0", + "100,0", + "0,-200", + "0,200", + ",0", + "0,", + "0;0", + "0,0,0", + "A,0", + "0,A", + NULL + }; + char* timeBAD[] = { + "2020-10-20t10:10:10", + "2020_10-20T10:10:10", + "2020-10_20T10:10:10", + "2020-10-20T10;10:10", + "2020-10-20T10:10;10", + "2020-19-20T10:10:10", + "2020-10-20T10:10:100", + "2020-10T10:10:10", + "2020-10-20T10:10", + "2020-A-20T10:10:10", + "2020-10-20T10:10:A", + "2020-10-20-5T10:10:10", + "2020-10-20T10:10:10:10", + "2020--20T10:10:10", + "2020-10-20T10::10", + "-10-20T10:10:10", + "2020-10-20T:10:10", + "2020-10-T10:10:10", + "2020-10-20T10:10:", + NULL + }; + size_t i; char** cmd = NULL; + p = strptime(timeOK, "%Y-%m-%dT%H:%M:%S", &ref); + CHK(p != NULL && *p == '\0'); + cmd = cmd_create(0, "test", "-D", "0,1", NULL); CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); CHK(args.nsun_dirs == 1); - CHK(args.sun_dirs[0].azimuth == 0); - CHK(eq_eps(args.sun_dirs[0].elevation, 1, 1.e-6) == 1); + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_SPHERICAL); + CHK(args.sun_dirs[0].spherical.azimuth == 0); + CHK(eq_eps(args.sun_dirs[0].spherical.elevation, 1, 1.e-6) == 1); solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", "1.2,3.4:3.14,0.123:", RES_OK); + cmd = cmd_create(0, "test", "-D", "1.2,3.4", "-D", "3.14,0.123", NULL); CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); CHK(args.nsun_dirs == 2); - CHK(eq_eps(args.sun_dirs[0].azimuth, 1.2, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[0].elevation, 3.4, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[1].azimuth, 3.14, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[1].elevation, 0.123, 1.e-6) == 1); + + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_SPHERICAL); + CHK(eq_eps(args.sun_dirs[0].spherical.azimuth, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[0].spherical.elevation, 3.4, 1.e-6) == 1); + CHK(args.sun_dirs[1].type == SOLSTICE_ARGS_SPHERICAL); + CHK(eq_eps(args.sun_dirs[1].spherical.azimuth, 3.14, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[1].spherical.elevation, 0.123, 1.e-6) == 1); solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", "1.2,3.4:3.14,0.123:2.01,23.1", RES_OK); + cmd = cmd_create(0, "test", "-L", "1.2,3.4", "-T", timeOK, NULL); CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); - CHK(args.nsun_dirs == 3); - CHK(eq_eps(args.sun_dirs[0].azimuth, 1.2, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[0].elevation, 3.4, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[1].azimuth, 3.14, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[1].elevation, 0.123, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[2].azimuth, 2.01, 1.e-6) == 1); - CHK(eq_eps(args.sun_dirs[2].elevation, 23.1, 1.e-6) == 1); + CHK(args.nsun_dirs == 1); + + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_LOCATION_TIME); + CHK(eq_eps(args.sun_dirs[0].location_time.latitude, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[0].location_time.longitude, 3.4, 1.e-6) == 1); + CHK(tm_equal(&args.sun_dirs[0].location_time.time, &ref) == 1); solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", "1.2,3.4,5", NULL); - CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); - cmd_delete(cmd); + cmd = cmd_create(0, "test", "-L", "1,4", "-L", "1.2,3.4", "-T", timeOK, NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); + CHK(args.nsun_dirs == 1); - cmd = cmd_create(0, "test", "-D", "1.2,3.4:1", NULL); - CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_LOCATION_TIME); + CHK(eq_eps(args.sun_dirs[0].location_time.latitude, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[0].location_time.longitude, 3.4, 1.e-6) == 1); + CHK(tm_equal(&args.sun_dirs[0].location_time.time, &ref) == 1); + solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", "1.2,3.4:5.2,A", NULL); - CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); - cmd_delete(cmd); + cmd = cmd_create(0, "test", "-L", "1.2,3.4", "-T", timeOK, "-T", timeOK, NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); + CHK(args.nsun_dirs == 2); - cmd = cmd_create(0, "test", "-D", "-0.1,2", NULL); - CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_LOCATION_TIME); + CHK(eq_eps(args.sun_dirs[0].location_time.latitude, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[0].location_time.longitude, 3.4, 1.e-6) == 1); + CHK(tm_equal(&args.sun_dirs[0].location_time.time, &ref) == 1); + CHK(args.sun_dirs[1].type == SOLSTICE_ARGS_LOCATION_TIME); + CHK(eq_eps(args.sun_dirs[1].location_time.latitude, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[1].location_time.longitude, 3.4, 1.e-6) == 1); + CHK(tm_equal(&args.sun_dirs[1].location_time.time, &ref) == 1); + solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", "360,2", NULL); - CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + cmd = cmd_create(0, "test", "-L", "1,3", "-T", timeOK, "-L", "2,4", "-T", timeOK, NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); + CHK(args.nsun_dirs == 2); + + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_LOCATION_TIME); + CHK(eq_eps(args.sun_dirs[0].location_time.latitude, 1, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[0].location_time.longitude, 3, 1.e-6) == 1); + CHK(tm_equal(&args.sun_dirs[0].location_time.time, &ref) == 1); + CHK(args.sun_dirs[1].type == SOLSTICE_ARGS_LOCATION_TIME); + CHK(eq_eps(args.sun_dirs[1].location_time.latitude, 2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[1].location_time.longitude, 4, 1.e-6) == 1); + CHK(tm_equal(&args.sun_dirs[1].location_time.time, &ref) == 1); + solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", "0,-1", NULL); - CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + cmd = cmd_create(0, "test", "-D", "1.2,3.4", "-L", "1.2,3.4", "-T", timeOK, NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); + CHK(args.nsun_dirs == 2); + + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_SPHERICAL); + CHK(eq_eps(args.sun_dirs[0].spherical.azimuth, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[0].spherical.elevation, 3.4, 1.e-6) == 1); + CHK(args.sun_dirs[1].type == SOLSTICE_ARGS_LOCATION_TIME); + CHK(eq_eps(args.sun_dirs[1].location_time.latitude, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[1].location_time.longitude, 3.4, 1.e-6) == 1); + CHK(tm_equal(&args.sun_dirs[1].location_time.time, &ref) == 1); + solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", "0,91", NULL); - CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + cmd = cmd_create(0, "test", "-D", "1.2,3.4", "-D", "3.14,0.123", "-D", "2.01,23.1", NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); + CHK(args.nsun_dirs == 3); + + CHK(args.sun_dirs[0].type == SOLSTICE_ARGS_SPHERICAL); + CHK(eq_eps(args.sun_dirs[0].spherical.azimuth, 1.2, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[0].spherical.elevation, 3.4, 1.e-6) == 1); + CHK(args.sun_dirs[1].type == SOLSTICE_ARGS_SPHERICAL); + CHK(eq_eps(args.sun_dirs[1].spherical.azimuth, 3.14, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[1].spherical.elevation, 0.123, 1.e-6) == 1); + CHK(args.sun_dirs[2].type == SOLSTICE_ARGS_SPHERICAL); + CHK(eq_eps(args.sun_dirs[2].spherical.azimuth, 2.01, 1.e-6) == 1); + CHK(eq_eps(args.sun_dirs[2].spherical.elevation, 23.1, 1.e-6) == 1); + solstice_args_release(&args); cmd_delete(cmd); - cmd = cmd_create(0, "test", "-D", NULL); + cmd = cmd_create(0, "test", "-T", timeOK, "-L", "1.2,3.4", NULL); CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); cmd_delete(cmd); + for(i = 0; i < sizeof(dirsBAD)/sizeof(*dirsBAD); i++) { + cmd = cmd_create(0, "test", "-D", dirsBAD[i], NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + cmd_delete(cmd); + } + + for(i = 0; i < sizeof(locsBAD)/sizeof(*locsBAD); i++) { + cmd = cmd_create(0, "test", "-L", locsBAD[i], NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + cmd_delete(cmd); + } + + for(i = 0; i < sizeof(timeBAD)/sizeof(*timeBAD); i++) { + cmd = cmd_create(0, "test", "-L", "1.2,3.4", "-T", timeBAD[i], NULL); + CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_BAD_ARG); + cmd_delete(cmd); + } + cmd = cmd_create(0, "test", "-D", "0,90", NULL); CHK(solstice_args_init(&args, cmd_size(cmd), cmd) == RES_OK); solstice_args_release(&args); cmd_delete(cmd); - } static void diff --git a/src/test_solstice_simulation_time.c b/src/test_solstice_simulation_time.c @@ -0,0 +1,613 @@ +/* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com) + * Copyright (C) 2016-2018 CNRS + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#define _XOPEN_SOURCE /* strptime support */ +#define _POSIX_C_SOURCE 200809L /* mkstemp support */ + +#include <star/scem.h> + +#include <rsys/rsys.h> +#include <rsys/math.h> +#include <rsys/double2.h> + +#include <time.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +#ifdef COMPILER_CL + /* Wrap POSIX functions and constants */ + #include <io.h> + #include <fcntl.h> + #include <sys/stat.h> + #define fdopen _fdopen + #define open _open +#endif + +enum side { + FRONT, + BACK +}; + +enum global_result_type { + GLOBAL_POTENTIAL, + GLOBAL_ABSORBED, + GLOBAL_COS, + GLOBAL_SHADOW, + GLOBAL_MISSING, + GLOBAL_ATMOSPHERE, + GLOBAL_REFLECTIVITY, + GLOBAL_RESULTS_COUNT__ +}; + +enum receiver_result_type { + FIRST_RECEIVER_RESULT, + FRONT_ABSORBED_FLUX = FIRST_RECEIVER_RESULT, + FRONT_INCOMING_FLUX, + FRONT_ABSORBED_FIELD_GAIN, + FRONT_ABSORBED_ATM_GAIN, + FRONT_EFFICIENCY, + BACK_ABSORBED_FLUX, + BACK_INCOMING_FLUX, + BACK_ABSORBED_FIELD_GAIN, + BACK_ABSORBED_ATM_GAIN, + BACK_EFFICIENCY, + RECEIVER_RESULTS_COUNT__ +}; + +enum primary_result_type { + FIRST_PRIMARY_RESULT, + PRIMARY_COS = FIRST_PRIMARY_RESULT, + PRIMARY_SHADOW, + PRIMARY_RESULTS_COUNT__ +}; + +struct counts { + unsigned long global, receiver, primary, realisation, failed; +}; + +static int +counts_ok(const struct counts* ref, const struct counts* c) +{ + CHK(ref->global == GLOBAL_RESULTS_COUNT__); + CHK(c->global == GLOBAL_RESULTS_COUNT__); + return ref->receiver == c->receiver + && ref->primary == c->primary + && ref->failed >= c->failed; +} + +#define MAX_LINE_LEN 2048 + +static const char +sundir_header [] = "#--- Sun direction:"; +static const char +suntime_header [] = "#--- Sun location and time:"; + + +#define IS_NEW_BLOCK(Line, Header) (!strncmp((Line), (Header), strlen(Header))) + + +#ifdef COMPILER_CL +/* mkstemp extracted from libc/sysdeps/posix/tempname.c. Copyright + * (C) 1991-1999, 2000, 2001, 2006 Free Software Foundation, Inc. + * + * The GNU C Library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. */ + +static const char letters [] = +"abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789"; + +/* Generate a temporary file name based on TMPL. TMPL must match the + * rules for mk[s]temp (i.e. end in "XXXXXX"). The name constructed + * does not exist at the time of the call to mkstemp. TMPL is + * overwritten with the result. */ +int +mkstemp(char *tmpl) +{ + size_t len; + char *XXXXXX; + static unsigned long long value; + unsigned long long random_time_bits; + unsigned int count; + int fd = -1; + int save_errno = errno; + + /* A lower bound on the number of temporary files to attempt to + * generate. The maximum total number of temporary file names that + * can exist for a given template is 62**6. It should never be + * necessary to try all these combinations. Instead if a reasonable + * number of names is tried (we define reasonable as 62**3) fail to + * give the system administrator the chance to remove the problems. */ + #define ATTEMPTS_MIN (62 * 62 * 62) + + /* The number of times to attempt to generate a temporary file. To + * conform to POSIX, this must be no smaller than TMP_MAX. */ + #if ATTEMPTS_MIN < TMP_MAX + unsigned int attempts = TMP_MAX; + #else + unsigned int attempts = ATTEMPTS_MIN; + #endif + + len = strlen(tmpl); + if (len < 6 || strcmp(&tmpl[len - 6], "XXXXXX")) { + errno = EINVAL; + return -1; + } + + /* This is where the Xs start. */ + XXXXXX = &tmpl[len - 6]; + + /* Get some more or less random data. */ + { + SYSTEMTIME stNow; + FILETIME ftNow; + + /* get system time */ + GetSystemTime(&stNow); + stNow.wMilliseconds = 500; + if (!SystemTimeToFileTime(&stNow, &ftNow)) { + errno = -1; + return -1; + } + + random_time_bits = (((unsigned long long)ftNow.dwHighDateTime << 32) + | (unsigned long long)ftNow.dwLowDateTime); + } + value += random_time_bits ^ (unsigned long long)GetCurrentThreadId(); + + for(count = 0; count < attempts; value += 7777, ++count) { + unsigned long long v = value; + + /* Fill in the random bits. */ + XXXXXX[0] = letters[v % 62]; + v /= 62; + XXXXXX[1] = letters[v % 62]; + v /= 62; + XXXXXX[2] = letters[v % 62]; + v /= 62; + XXXXXX[3] = letters[v % 62]; + v /= 62; + XXXXXX[4] = letters[v % 62]; + v /= 62; + XXXXXX[5] = letters[v % 62]; + + fd = open(tmpl, O_RDWR|O_CREAT|O_EXCL, S_IREAD|S_IWRITE); + if (fd >= 0) { + errno = save_errno; + return fd; + } + else if (errno != EEXIST) + return -1; + } + + /* We got out of the loop because we ran out of combinations to try. */ + errno = EEXIST; + return -1; +} +#endif + +static int +read_line(char* line, size_t max_line_len, FILE* stream) +{ + ASSERT(stream && line && max_line_len); + line = fgets(line, (int)max_line_len, stream); + if(!line) return 0; + CHK(strlen(line) + 1 < max_line_len); + return 1; +} + +static int +get_angles_and_counts + (FILE* file, + double angles[2], + struct counts* counts) +{ + char line[MAX_LINE_LEN]; + int r; + + CHK(file != NULL); + CHK(angles != NULL); + CHK(counts != NULL); + + /* Get sun dir */ + r = read_line(line, sizeof(line), file); + if (!r) { + CHK(feof(file) == 1); + return 0; + } + CHK(IS_NEW_BLOCK(line, sundir_header) == 1); + CHK(sscanf(line+strlen(sundir_header), "%lg%lg", &angles[0], &angles[1]) == 2); + + /* Get counts */ + CHK(read_line(line, sizeof(line), file) == 1); + CHK( + sscanf(line, + "%lu %lu %lu %lu %lu", + &counts->global, &counts->receiver, &counts->primary, + &counts->realisation, &counts->failed) == 5); + return 1; +} + +static int +get_dir_and_counts + (FILE* file, + double dir[3], + struct counts* counts) +{ + char line[MAX_LINE_LEN]; + int r, n; + + CHK(file != NULL); + CHK(dir != NULL); + CHK(counts != NULL); + + /* Get sun dir */ + r = read_line(line, sizeof(line), file); + if (!r) { + CHK(feof(file) == 1); + return 0; + } + CHK(IS_NEW_BLOCK(line, suntime_header) == 1); + n = sscanf(line+strlen(suntime_header), + "%*g %*g %*d-%*d-%*dT%*d:%*d:%*d (%lg %lg %lg)", + dir, dir+1, dir+2); + CHK(n == 3); + + /* Get counts */ + CHK(read_line(line, sizeof(line), file) == 1); + CHK( + sscanf(line, + "%lu %lu %lu %lu %lu", + &counts->global, &counts->receiver, &counts->primary, + &counts->realisation, &counts->failed) == 5); + return 1; +} + +static void +read_global(FILE* file, double* E, double* SE) +{ + char line[MAX_LINE_LEN]; + CHK(read_line(line, sizeof(line), file) == 1); + CHK(sscanf(line, "%lg %lg", E, SE) == 2); +} + +static void +read_recv(FILE* file, char name[], double E[], double SE[]) +{ + char line[MAX_LINE_LEN]; + + CHK(file != NULL); + CHK(name != NULL); + CHK(E != NULL); + CHK(SE != NULL); + + CHK(read_line(line, sizeof(line), file) == 1); + CHK( + sscanf(line, + "%s %*u %*g " + "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg " + "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", + name, /* ID, area */ + &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], + &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], + &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], + &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], + &E[FRONT_EFFICIENCY], &SE[FRONT_EFFICIENCY], + &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], + &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], + &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], + &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN], + &E[BACK_EFFICIENCY], &SE[BACK_EFFICIENCY]) == + 2 * RECEIVER_RESULTS_COUNT__ + 1); +} + +static void +read_primary + (FILE* file, char name[], double* area, double E[], double SE[]) +{ + char line[MAX_LINE_LEN]; + + CHK(file != NULL); + CHK(area != NULL); + CHK(E != NULL); + CHK(SE != NULL); + + CHK(read_line(line, sizeof(line), file) == 1); + CHK( + sscanf(line, + "%s %*u " + "%lg %*u " + "%lg %lg %lg %lg\n", + name, /* ID */ + area, /* count, */ + &E[PRIMARY_COS], &SE[PRIMARY_COS], + &E[PRIMARY_SHADOW], &SE[PRIMARY_SHADOW]) == + 2 * PRIMARY_RESULTS_COUNT__ + 2); +} + + +static void +read_recvXprim + (FILE* file, + unsigned long* rcv_id, + unsigned long* prim_id, + double E[], + double SE[]) +{ + char line[MAX_LINE_LEN]; + + CHK(file != NULL); + CHK(rcv_id != NULL); + CHK(prim_id != NULL); + CHK(E != NULL); + CHK(SE != NULL); + + CHK(read_line(line, sizeof(line), file) == 1); + CHK( + sscanf(line, + "%lu %lu " + "%lg %lg %lg %lg %lg %lg %lg %lg " + "%lg %lg %lg %lg %lg %lg %lg %lg", + rcv_id, prim_id, + &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], + &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], + &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], + &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], + &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], + &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], + &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], + &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN]) == + 2 * (RECEIVER_RESULTS_COUNT__ - 2 /* efficiencies not read */) + 2); +} + +static void +compute_estimate_intersection + (double intersection[2], + const double scale, + const double E0, + const double SE0, + const double E1, + const double SE1) +{ + double interval0[2], interval1[2]; + CHK(scale > 0); + interval0[0] = E0 - scale*SE0; + interval0[1] = E0 + scale*SE0; + interval1[0] = E1 - scale*SE1; + interval1[1] = E1 + scale*SE1; + intersection[0] = MMAX(interval0[0], interval1[0]); + intersection[1] = MMIN(interval0[1], interval1[1]); +} + +static void +check_estimate + (double ref_E, + double ref_SE, + double test_E, + double test_SE) +{ + if(ref_E == -1) { + CHK(ref_SE == -1); + CHK(test_E == -1); + CHK(test_SE == -1); + } else { + double interval[2]; + CHK(ref_SE >= 0); + CHK(test_E >= 0); + CHK(test_SE >= 0); + if(ref_SE == 0) ref_SE = ref_E / 1000.0; + if(test_SE == 0) test_SE = test_E / 1000.0; + compute_estimate_intersection(interval, 2, ref_E, ref_SE, test_E, test_SE); + CHK(interval[0] <= interval[1]); + } +} + +static void +check_1_reference + (FILE* ref_file, + FILE* test_file, + const struct counts* counts) +{ + unsigned n; + + CHK(ref_file != NULL); + CHK(test_file != NULL); + CHK(counts != NULL); + + /* both files' pointer are just past the new bloc header */ + + for(n = 0; n < counts->global; n++) { + double reference_E, reference_SE, test_E, test_SE; + read_global(ref_file, &reference_E, &reference_SE); + read_global(test_file, &test_E, &test_SE); + check_estimate(reference_E, reference_SE, test_E, test_SE); + } + for(n = 0; n < counts->receiver; n++) { + char ref_rcv_name[MAX_LINE_LEN], test_rcv_name[MAX_LINE_LEN]; + double reference_E[RECEIVER_RESULTS_COUNT__]; + double reference_SE[RECEIVER_RESULTS_COUNT__]; + double test_E[RECEIVER_RESULTS_COUNT__]; + double test_SE[RECEIVER_RESULTS_COUNT__]; + enum receiver_result_type r; + + read_recv(ref_file, ref_rcv_name, reference_E, reference_SE); + read_recv(test_file, test_rcv_name, test_E, test_SE); + CHK(strcmp(ref_rcv_name, test_rcv_name) == 0); + FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { + check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); + } + } + for(n = 0; n < counts->primary; n++) { + char ref_prim_name[MAX_LINE_LEN], test_prim_name[MAX_LINE_LEN]; + double reference_E[PRIMARY_RESULTS_COUNT__]; + double reference_SE[PRIMARY_RESULTS_COUNT__]; + double test_E[PRIMARY_RESULTS_COUNT__]; + double test_SE[PRIMARY_RESULTS_COUNT__]; + double ref_area, test_area; + enum primary_result_type r; + + read_primary(ref_file, ref_prim_name, &ref_area, reference_E, reference_SE); + read_primary(test_file, test_prim_name, &test_area, test_E, test_SE); + check_estimate(ref_area, 0, test_area, 0); + CHK(strcmp(ref_prim_name, test_prim_name) == 0); + FOR_EACH(r, FIRST_PRIMARY_RESULT, PRIMARY_RESULTS_COUNT__) { + check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); + } + } + for(n = 0; n < counts->receiver * counts->primary; n++) { + double reference_E[RECEIVER_RESULTS_COUNT__]; + double reference_SE[RECEIVER_RESULTS_COUNT__]; + double test_E[RECEIVER_RESULTS_COUNT__]; + double test_SE[RECEIVER_RESULTS_COUNT__]; + unsigned long ref_rcv_id, ref_prim_id; + unsigned long test_rcv_id, test_prim_id; + + enum receiver_result_type r; + read_recvXprim(ref_file, &ref_rcv_id, &ref_prim_id, reference_E, reference_SE); + read_recvXprim(test_file, &test_rcv_id, &test_prim_id, test_E, test_SE); + /* we rely on the order of outputs */ + CHK(ref_rcv_id == test_rcv_id); + CHK(ref_prim_id == test_prim_id); + FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { + if (r == FRONT_EFFICIENCY || r == BACK_EFFICIENCY) + continue; /* not read */ + check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); + } + } +} + +static FINLINE int +create_tmp_file(char* name, const size_t max_sizeof_name) +{ + const char* template = "solstice_tmp_file_XXXXXX"; + int fd; + CHK(name != NULL); + CHK(strlen(template)+1 <= max_sizeof_name-1); + strcpy(name, template); + fd = mkstemp(name); + CHK(fd != -1); + return fd; +} + +/* Scem an solstice agree on elevation. + * On the other hand, scem and solstice convention on azimuth differ: + * - scem CW with N=0 + * - solstice CCW with E=0 */ +static FINLINE void +scem_pos_to_solstice_pos + (const struct scem_sun_pos* scem, double solstice[2]) +{ + ASSERT(scem && solstice); + solstice[0] = fmod(360 + 90 - MRAD2DEG(scem->azimuth), 360); + solstice[1] = MRAD2DEG(scem->elevation); +} + +static void +do_check(const char* binary, const char* dir, const char* base_name) +{ + struct counts time_counts, dir_counts; + int n, h; + int err; + double lat, lon; + char cmd[512]; + char time_file_name[128]; + char dir_file_name[128]; + char sdate[64]; + char* p; + struct tm date; + struct scem_sun_pos scem_pos; + struct scem_location loc; + enum scem_sun_algo algorithm; + res_T res; + double sun_dir[3], dir_angles[2], solstice_pos[2]; + FILE* dir_file; + FILE* time_file; + int dir_fd, time_fd; + + lat = -90 + 180 * ((double)rand() / (double)RAND_MAX); + lon = 15 * ((double)rand() / (double)RAND_MAX); /* To allow to sample UTC h */ + h = (int)(9 + 10 * ((double)rand() / (double)RAND_MAX)); + algorithm = SCEM_SUN_PSA; + + snprintf(sdate, sizeof(sdate), "2020-05-05T%d:00:00", h); + p = strptime(sdate, "%Y-%m-%dT%H:%M:%S", &date); + CHK(p && *p == '\0'); + loc.latitude = lat; loc.longitude = lon; + + time_fd = create_tmp_file(time_file_name, sizeof(time_file_name)); + time_file = fdopen(time_fd, "r"); + CHK(time_file != NULL); + + n = snprintf(cmd, sizeof(cmd), + "%s -o %s -f -L %g,%g -T %s -n %lu -R %s%s_receiver.yaml %s%s.yaml", + binary, time_file_name, lat, lon, sdate, 10000L, + dir, base_name, dir, base_name); + CHK((unsigned)n < sizeof(cmd)); + + err = system(cmd); + CHK(err == 0); + + get_dir_and_counts(time_file, sun_dir, &time_counts); + + dir_fd = create_tmp_file(dir_file_name, sizeof(dir_file_name)); + dir_file = fdopen(dir_fd, "r"); + CHK(dir_file != NULL); + + res = scem_sun_position_from_earth(&date, &loc, algorithm, &scem_pos); + CHK(res == RES_OK); + scem_pos_to_solstice_pos(&scem_pos, solstice_pos); + n = snprintf(cmd, sizeof(cmd), + "%s -o %s -f -D %g,%g -n %lu -R %s%s_receiver.yaml %s%s.yaml", + binary, dir_file_name, + SPLIT2(solstice_pos), + 10000L, dir, base_name, dir, base_name); + CHK((unsigned)n < sizeof(cmd)); + + err = system(cmd); + CHK(err == 0); + + get_angles_and_counts(dir_file, dir_angles, &dir_counts); + + CHK(counts_ok(&time_counts, &dir_counts) == 1); + check_1_reference(time_file, dir_file, &time_counts); + + fclose(time_file); + fclose(dir_file); + remove(time_file_name); + remove(dir_file_name); +} + +int +main(int argc, char** argv) +{ + int err = 0; + + if(argc != 4) { + printf("Usage: %s <solstice-binary> <file-path> <file-base-name>\n", argv[0]); + goto error; + } + + do_check(argv[1], argv[2], argv[3]); + +exit: + return err; +error: + err = 1; + goto exit; +} +