solstice-solver

commit e5e169fc2de7a159c3d8d06063e38c1c896ce525
parent 581e30fff97dead628f48867b40f789f381c6f9e
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Fri, 29 May 2026 14:50:20 +0200

Add location and date as a way to set sun dir

A second way to set the sun direction through a location and a date.
The algorithm used to compute the direction can also be selected.

Diffstat:
M
Makefile
 | 
3
+++
M
config.mk
 | 
5
+++--
M
src/ssol.h
 | 
29
+++++++++++++++++++++++++++++
M
src/ssol_sun.c
 | 
98
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
M
src/ssol_sun_c.h
 | 
3
+++
M
src/test_ssol_sun.c
 | 
67
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--

6 files changed, 201 insertions(+), 4 deletions(-)
diff --git a/Makefile b/Makefile
@@ -89,6 +89,7 @@ src/.config: config.mk
 	$(PKG_CONFIG) --atleast-version $(SCPR_VERSION) scpr
 	$(PKG_CONFIG) --atleast-version $(SSF_VERSION) ssf
 	$(PKG_CONFIG) --atleast-version $(SSP_VERSION) star-sp
+	$(PKG_CONFIG) --atleast-version $(SCEM_VERSION) scem
 	@echo "config done" > $@
 
 $(DEP) : $(HDR)
@@ -111,6 +112,7 @@ src/ssol.pc: ssol.pc.in
 	    -e 's#@SCPR_VERSION@#$(SCPR_VERSION)#g'\
 	    -e 's#@SSF_VERSION@#$(SSF_VERSION)#g'\
 	    -e 's#@SSP_VERSION@#$(SSP_VERSION)#g'\
+	    -e 's#@SCEM_VERSION@#$(SCEM_VERSION)#g'\
 	    ssol.pc.in > $@
 
 src/ssol-local.pc: ssol.pc.in
@@ -123,6 +125,7 @@ src/ssol-local.pc: ssol.pc.in
 	    -e 's#@SCPR_VERSION@#$(SCPR_VERSION)#g'\
 	    -e 's#@SSF_VERSION@#$(SSF_VERSION)#g'\
 	    -e 's#@SSP_VERSION@#$(SSP_VERSION)#g'\
+	    -e 's#@SCEM_VERSION@#$(SCEM_VERSION)#g'\
 	    ssol.pc.in > $@
 
 src/ssol_version.h: src/ssol_version.h.in
diff --git a/config.mk b/config.mk
@@ -36,10 +36,11 @@ S3DUT_VERSION = 0.4
 SCPR_VERSION = 0.5
 SSF_VERSION = 0.10
 SSP_VERSION = 0.15
+SCEM_VERSION = 0.0
 
-INCS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags rsys s3d s3dut scpr ssf star-sp)\
+INCS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags rsys s3d s3dut scpr ssf star-sp scem)\
  -fopenmp
-LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs rsys s3d s3dut scpr ssf star-sp)\
+LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs rsys s3d s3dut scpr ssf star-sp scem)\
  -fopenmp -lm
 
 ################################################################################
diff --git a/src/ssol.h b/src/ssol.h
@@ -45,6 +45,7 @@
 struct logger;
 struct mem_allocator;
 struct ssp_rng;
+struct tm;
 
 /* Opaque Solstice solver types */
 struct ssol_atmosphere;
@@ -132,6 +133,12 @@ enum ssol_data_type {
   SSOL_DATA_SPECTRUM
 };
 
+/* The algorithm used to produce sun direction from location and date */
+enum ssol_sun_direction_algorithm {
+ SSOL_SUN_DIRECTION_ALGORITHM_MEEUS,
+ SSOL_SUN_DIRECTION_ALGORITHM_PSA
+};
+
 /* Describe a vertex data */
 struct ssol_vertex_data {
   enum ssol_attrib_usage usage; /* Semantic of the data */
@@ -1032,6 +1039,28 @@ ssol_sun_set_direction
   (struct ssol_sun* sun,
    const double direction[3]);
 
+/* Main sun direction defined by location and date.
+ * The algorithm must be set first */
+SSOL_API res_T
+ssol_sun_set_location_and_date
+  (struct ssol_sun* sun,
+   const double latitude,
+   const double longitude,
+   const struct tm* date);
+
+/* Set the algorithm used to produce sun direction from location and date.
+ * The algorithm must be set before calling ssol_sun_set_location_and_date */
+SSOL_API res_T
+ssol_sun_set_algorithm
+  (struct ssol_sun* sun,
+   const enum ssol_sun_direction_algorithm algorithm);
+
+SSOL_API res_T
+ssol_sun_get_algorithm
+  (const struct ssol_sun* sun,
+   enum ssol_sun_direction_algorithm* algorithm);
+
+/* Get direction from sun */
 SSOL_API res_T
 ssol_sun_get_direction
   (const struct ssol_sun* sun,
diff --git a/src/ssol_sun.c b/src/ssol_sun.c
@@ -21,6 +21,8 @@
 #include "ssol_ranst_sun_wl.h"
 #include "ssol_spectrum_c.h"
 
+#include <star/scem.h>
+
 #include <rsys/rsys.h>
 #include <rsys/mem_allocator.h>
 #include <rsys/ref_count.h>
@@ -65,6 +67,7 @@ sun_create
   SSOL(device_ref_get(dev));
   sun->dev = dev;
   sun->type = type;
+  sun->algorithm = SSOL_SUN_DIRECTION_ALGORITHM_MEEUS;
   ref_init(&sun->ref);
 
 exit:
@@ -78,6 +81,33 @@ error:
   goto exit;
 }
 
+static INLINE void
+scem_to_cartesian_sun_dir
+  (struct scem_sun_pos* scem, double sun_dir[3])
+{
+  double cos_azimuth;
+  double sin_azimuth;
+  double cos_zenith;
+  double sin_zenith;
+  double azimuth;
+  double zenith;
+  ASSERT(scem && sun_dir);
+
+  /* scem->azimuth: East = 0, CCW
+   * scem->zenith: XY plane = 0, UP = +PI/2, DOWN = -PI/2 */
+  azimuth = scem->azimuth;
+  zenith = scem->zenith;
+
+  cos_azimuth = cos(azimuth);
+  sin_azimuth = sin(azimuth);
+  cos_zenith = cos(zenith);
+  sin_zenith = sin(zenith);
+
+  sun_dir[0] = -(cos_zenith * cos_azimuth);
+  sun_dir[1] = -(cos_zenith * sin_azimuth);
+  sun_dir[2] = -(sin_zenith);
+}
+
 /*******************************************************************************
  * Exported ssol_image functions
  ******************************************************************************/
@@ -136,6 +166,74 @@ ssol_sun_set_direction(struct ssol_sun* sun, const double direction[3])
 }
 
 res_T
+ssol_sun_set_location_and_date
+  (struct ssol_sun* sun,
+   const double latitude,
+   const double longitude,
+   const struct tm* date)
+{
+  res_T res = RES_OK;
+  struct scem_sun_pos sun_pos;
+  struct scem_location loc;
+  enum scem_sun_algo algorithm;
+
+  if(!sun || !date) {
+    res = RES_BAD_ARG;
+    goto error;
+  }
+  if(latitude < -90 || latitude > 90 || longitude < -180 || longitude > 180) {
+    res = RES_BAD_ARG;
+    goto error;
+  }
+
+  loc.latitude = latitude;
+  loc.longitude = longitude;
+
+  switch(sun->algorithm ) {
+    case SSOL_SUN_DIRECTION_ALGORITHM_PSA :
+      algorithm = SCEM_SUN_PSA ;
+      break;
+    case SSOL_SUN_DIRECTION_ALGORITHM_MEEUS :
+      algorithm = SCEM_SUN_MEEUS ;
+      break;
+    default: FATAL("Unknown enum value.");
+  }
+  res = scem_sun_position_from_earth(date, &loc, algorithm, &sun_pos);
+  if(res != RES_OK) goto error;
+
+  scem_to_cartesian_sun_dir(&sun_pos, sun->direction);
+
+exit:
+  return res;
+error:
+  goto exit;
+}
+
+res_T
+ssol_sun_set_algorithm
+  (struct ssol_sun* sun,
+   const enum ssol_sun_direction_algorithm algorithm)
+{
+  if(!sun ||
+      (algorithm != SSOL_SUN_DIRECTION_ALGORITHM_MEEUS
+       && algorithm != SSOL_SUN_DIRECTION_ALGORITHM_PSA))
+    return RES_BAD_ARG;
+  sun->algorithm = algorithm;
+  return RES_OK;
+}
+
+res_T
+ssol_sun_get_algorithm
+  (const struct ssol_sun* sun,
+   enum ssol_sun_direction_algorithm* algorithm)
+{
+  if(!sun || !algorithm)
+    return RES_BAD_ARG;
+  *algorithm = sun->algorithm;
+  return RES_OK;
+}
+
+res_T
 ssol_sun_get_direction(const struct ssol_sun* sun, double direction[3])
 {
   if(!sun || !direction)
diff --git a/src/ssol_sun_c.h b/src/ssol_sun_c.h
@@ -17,6 +17,8 @@
 #ifndef SSOL_SUN_C_H
 #define SSOL_SUN_C_H
 
+#include "ssol.h"
+
 #include <rsys/ref_count.h>
 #include <rsys/list.h>
 
@@ -55,6 +57,7 @@ struct ssol_sun {
     struct gaussian gaussian;
     struct buie csr;
   } data;
+  enum ssol_sun_direction_algorithm algorithm;
 
   struct ssol_device* dev;
   ref_T ref;
diff --git a/src/test_ssol_sun.c b/src/test_ssol_sun.c
@@ -14,11 +14,15 @@
  * 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 "ssol.h"
 #include "test_ssol_utils.h"
 
 #include <rsys/double3.h>
 
+#include <time.h>
+
 int
 main(int argc, char** argv)
 {
@@ -29,8 +33,10 @@ main(int argc, char** argv)
   struct ssol_sun* sun;
   const double dir0[3] = { 0, 0, 0 };
   double dir[3] = { 1, 0, 0 };
-  double tmp[3];
+  double tmp[3], res[3];
   double dni;
+  struct tm date;
+  enum ssol_sun_direction_algorithm algo;
   (void) argc, (void) argv;
 
   mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -60,13 +66,70 @@ main(int argc, char** argv)
   CHK(ssol_sun_set_direction(sun, NULL) == RES_BAD_ARG);
   CHK(ssol_sun_set_direction(sun, dir0) == RES_BAD_ARG);
   CHK(ssol_sun_set_direction(sun, dir) == RES_OK);
-  CHK(ssol_sun_set_direction(sun, dir) == RES_OK);
 
   CHK(ssol_sun_get_direction(NULL, tmp) == RES_BAD_ARG);
   CHK(ssol_sun_get_direction(sun, NULL) == RES_BAD_ARG);
   CHK(ssol_sun_get_direction(sun, tmp) == RES_OK);
   CHK(d3_eq(dir, tmp) == 1);
 
+  CHK(ssol_sun_set_algorithm(NULL, 999) == RES_BAD_ARG);
+  CHK(ssol_sun_set_algorithm(NULL, SSOL_SUN_DIRECTION_ALGORITHM_MEEUS) == RES_BAD_ARG);
+  CHK(ssol_sun_set_algorithm(sun, 999) == RES_BAD_ARG);
+  CHK(ssol_sun_set_algorithm(sun, SSOL_SUN_DIRECTION_ALGORITHM_PSA) == RES_OK);
+
+  CHK(ssol_sun_get_algorithm(NULL, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_get_algorithm(NULL, &algo) == RES_BAD_ARG);
+  CHK(ssol_sun_get_algorithm(sun, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_get_algorithm(sun, &algo) == RES_OK);
+  CHK(algo == SSOL_SUN_DIRECTION_ALGORITHM_PSA);
+
+  CHK(ssol_sun_set_algorithm(sun, SSOL_SUN_DIRECTION_ALGORITHM_MEEUS) == RES_OK);
+  strptime("2026-6-21T12:4:0", "%Y-%m-%dT%H:%M:%S", &date);
+  CHK(ssol_sun_set_location_and_date(NULL, 100, 200, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(NULL, 100, 200, &date) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(NULL, 100, 0, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(NULL, 100, 0, &date) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(NULL, 0, 200, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(NULL, 0, 200, &date) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(NULL, 0, 0, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(NULL, 0, 0, &date) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 100, 200, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 100, 200, &date) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 100, 0, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 100, 0, &date) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 0, 200, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 0, 200, &date) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 0, 0, NULL) == RES_BAD_ARG);
+  CHK(ssol_sun_set_location_and_date(sun, 23.436, 0, &date) == RES_OK);
+
+  CHK(ssol_sun_get_direction(sun, tmp) == RES_OK);
+  d3(res, 0,0,-1);
+  CHK(d3_eq_eps(res, tmp, 1e-4) == 1);
+
+  CHK(ssol_sun_set_algorithm(sun, SSOL_SUN_DIRECTION_ALGORITHM_PSA) == RES_OK);
+  strptime("2026-6-21T12:1:49", "%Y-%m-%dT%H:%M:%S", &date);
+  CHK(ssol_sun_set_location_and_date(sun, 23.436, 0, &date) == RES_OK);
+
+  CHK(ssol_sun_get_direction(sun, tmp) == RES_OK);
+  d3(res, 0,0,-1);
+  CHK(d3_eq_eps(res, tmp, 1e-4) == 1);
+
+  CHK(ssol_sun_set_algorithm(sun, SSOL_SUN_DIRECTION_ALGORITHM_MEEUS) == RES_OK);
+  strptime("2026-6-21T10:0:0", "%Y-%m-%dT%H:%M:%S", &date);
+  CHK(ssol_sun_set_location_and_date(sun, 0, 0, &date) == RES_OK);
+
+  CHK(ssol_sun_get_direction(sun, tmp) == RES_OK);
+  d3(res, -0.397792,-0.47252,-0.78643);
+  CHK(d3_eq_eps(res, tmp, 1e-4) == 1);
+
+  CHK(ssol_sun_set_algorithm(sun, SSOL_SUN_DIRECTION_ALGORITHM_PSA) == RES_OK);
+  strptime("2026-6-21T10:0:0", "%Y-%m-%dT%H:%M:%S", &date);
+  CHK(ssol_sun_set_location_and_date(sun, 0, 0, &date) == RES_OK);
+
+  CHK(ssol_sun_get_direction(sun, tmp) == RES_OK);
+  d3(res, -0.397733,-0.465014,-0.79093);
+  CHK(d3_eq_eps(res, tmp, 1e-4) == 1);
+
   CHK(ssol_sun_set_dni(NULL, 1000) == RES_BAD_ARG);
   CHK(ssol_sun_set_dni(sun, 0) == RES_BAD_ARG);
   CHK(ssol_sun_set_dni(sun, 1000) == RES_OK);