commit 66589442b2056899e5abf99d4a8712544df9ebd9
parent 8a8e9dbbd2c0129ce966cc5d37a87684de508834
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 10 Oct 2025 12:09:45 +0200
stardis: handle sun as an external source
Implement the plugin functions to calculate solar flux from an external
source rather than an imposed flux. To do this, at a given time, the
position of the sun is calculated dynamically at the location where the
meteorological data is provided. Solar power is calculated from the DNI
retrieved from the meteorological. Finally, the solar luminance
scattered at least once in the atmosphere is calculated from the diffuse
component of the shortwave downward flux.
Diffstat:
7 files changed, 193 insertions(+), 2 deletions(-)
diff --git a/Makefile b/Makefile
@@ -101,8 +101,8 @@ PLUGIN_SRC =\
PLUGIN_OBJ = $(PLUGIN_SRC:.c=.o)
PLUGIN_DEP = $(PLUGIN_SRC:.c=.d)
-INCS_PLUGIN = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys smeteo-local stardis)
-LIBS_PLUGIN = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys smeteo-local)
+INCS_PLUGIN = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys smeteo-local scem stardis)
+LIBS_PLUGIN = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys smeteo-local scem) -lm
CFLAGS_PLUGIN = -std=c89 $(CFLAGS_SO) $(INCS_PLUGIN)
LDFLAGS_PLUGIN = $(LDFLAGS_SO) $(LIBS_PLUGIN)
@@ -115,6 +115,7 @@ libstardis_smeteo.so $(PLUGIN_DEP) $(PLUGIN_OBJ): \
config.mk smeteo-local.pc .config_plugin
.config_plugin:
+ $(PKG_CONFIG) --atleast-version $(SCEM_VERSION) scem
$(PKG_CONFIG) --atleast-version $(STARDIS_VERSION) stardis
@echo 'config done' > $@
diff --git a/README.md b/README.md
@@ -11,6 +11,11 @@ See smeteo.5 for format specification.
- [RSys](https://gitlab.com/vaplv/rsys)
- [mandoc](https://mandoc.bsd.lv)
+Additional requirements for the optional Stardis plugin:
+
+- [Stardis](https://gitlab.com/meso-star/stardis)
+- [Star-CelestialMechanics](https://gitlab.com/meso-star/star-cem)
+
## Installation
Edit config.mk as needed, then run:
diff --git a/config.mk b/config.mk
@@ -38,6 +38,7 @@ LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs rsys)
# For Stardis plugin
STARDIS_VERSION = 0.11
+SCEM_VERSION = 0.0
################################################################################
# Compilation options
diff --git a/src/stardis_smeteo.c b/src/stardis_smeteo.c
@@ -13,10 +13,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 _POSIX_C_SOURCE 199606L /* localtime_r support */
+
#include "smeteo.h"
#include "stardis_smeteo.h"
#include "stardis_smeteo_library.h"
+#include <star/scem.h>
+
#include <rsys/algorithm.h>
#include <rsys/math.h>
#include <rsys/mem_allocator.h>
@@ -25,6 +29,10 @@
* header. It should by stardis, not by the caller. */
#include <limits.h>
+#include <time.h> /* localtime_r */
+
+#define EARTH_TO_SUN_DST 149.5978707e9 /* [m] */
+
struct boundary_condition {
struct stardis_smeteo_lib* lib;
struct stardis_smeteo_lib_desc lib_desc;
@@ -116,6 +124,78 @@ get_meteo_entry_id
}
}
+static void
+spherical_to_cartesian_dir
+ (const double azimuth, /* [In radians */
+ const double elevation, /* In radians */
+ double dir[3])
+{
+ double cos_azimuth;
+ double sin_azimuth;
+ double cos_elevation;
+ double sin_elevation;
+ ASSERT(azimuth >= 0 && azimuth < 2*PI);
+ ASSERT(elevation >= -PI/2.0 && elevation <= PI/2.0);
+ ASSERT(dir);
+
+ cos_azimuth = cos(azimuth);
+ sin_azimuth = sin(azimuth);
+ cos_elevation = cos(elevation);
+ sin_elevation = sin(elevation);
+
+ dir[0] = cos_elevation * cos_azimuth;
+ dir[1] = cos_elevation * sin_azimuth;
+ dir[2] = sin_elevation;
+}
+
+static void
+compute_sun_position
+ (const struct boundary_condition* bcond,
+ const double time, /* [s] */
+ struct scem_sun_pos* sun)
+{
+/* Star-CElestial-Mechanics */
+ struct scem_location pos = SCEM_LOCATION_NULL;
+
+ /* Time and date */
+ struct tm date_utc0 = {0};
+ time_t time_since_jan_1_1850_utc0 = 0; /* [s] */
+ time_t time_since_jan_1_1850_local = 0; /* [s] */
+
+ /* Miscellaneous */
+ size_t i = 0; /* Index of the meteo entry including "time" */
+ res_T res = RES_OK;
+
+ ASSERT(bcond && sun);
+
+ i = get_meteo_entry_id(bcond, time);
+
+ /* Calculate the number of seconds elapsed since
+ * January 1, 1850, local time */
+ time_since_jan_1_1850_local = /* [s] */
+ (time_t)(bcond->lib_desc.smeteo_desc.entries[i].day_1850 * 24 * 3600);
+
+ /* The limit condition stores the number of seconds elapsed since the epoch
+ * until January 1, 1850, local time. In other words, since the epoch is
+ * expressed in UTC+00:00, it defines an offset in seconds between local
+ * time and UTC+00:00 time coordinates. Thus, add this offset to the
+ * number of seconds elapsed since January 1, 1850, local time, to convert
+ * this number of seconds to UTC+00:00. */
+ time_since_jan_1_1850_utc0 = /* [s] */
+ bcond->lib_desc.jan_1_1850 + time_since_jan_1_1850_local;
+
+ /* Convert time in number of seconds into a broken-down time as expected by
+ * the Star-CElestial-Mechanics library */
+ CHK(localtime_r(&time_since_jan_1_1850_utc0, &date_utc0) != 0);
+
+ /* Compute the solar position at the specified location and date */
+ pos.latitude = bcond->lib_desc.smeteo_desc.latitude;
+ pos.longitude = bcond->lib_desc.smeteo_desc.longitude;
+ res = scem_sun_position_from_earth
+ (&date_utc0, &pos, bcond->lib_desc.algo, sun);
+ CHK(res == RES_OK);
+}
+
/*******************************************************************************
* Exported symbols
******************************************************************************/
@@ -301,6 +381,60 @@ stardis_t_range(void* data, double range[2] /* [K] */)
return range;
}
+double*
+stardis_spherical_source_position
+ (const double time, /* [s] */
+ double sun_pos[3], /* [m] */
+ void* data)
+{
+ struct scem_sun_pos sun = SCEM_SUN_POS_NULL;
+ double sun_dir[3] = {0,0,0};
+ ASSERT(sun_pos);
+
+ compute_sun_position(data, time, &sun);
+
+ /* Convert the sun's position to Cartesian coords as expected by Stardis */
+ spherical_to_cartesian_dir(sun.azimuth, sun.zenith, sun_dir);
+ sun_pos[0] = sun_dir[0] * EARTH_TO_SUN_DST;
+ sun_pos[1] = sun_dir[1] * EARTH_TO_SUN_DST;
+ sun_pos[2] = sun_dir[2] * EARTH_TO_SUN_DST;
+ return sun_pos;
+}
+
+double /* [W] */
+stardis_spherical_source_power(const double time /* [s] */, void* data)
+{
+ struct boundary_condition* bcond = data;
+ struct scem_sun_pos sun = SCEM_SUN_POS_NULL;
+ double power = 0; /* [W] */
+ size_t i = 0; /* Index of the meteo entry including "time" */
+ ASSERT(data);
+
+ compute_sun_position(data, time, &sun);
+
+ i = get_meteo_entry_id(bcond, time);
+ power = bcond->lib_desc.smeteo_desc.entries[i].SWdn_direct / sin(sun.zenith);
+ return power;
+}
+
+double /* [W/m^2/sr] */
+stardis_spherical_source_diffuse_radiance
+ (const double time /* [s] */,
+ const double dir[3],
+ void* data)
+{
+ double diffuse_radiance = 0; /* [W/m^2/sr] */
+ struct boundary_condition* bcond = data;
+ size_t i = 0; /* Index of the meteo entry including "time" */
+
+ ASSERT(data);
+ (void)dir; /* Avoid "unused variable" warning */
+
+ i = get_meteo_entry_id(bcond, time);
+ diffuse_radiance = bcond->lib_desc.smeteo_desc.entries[i].SWdn_diffuse / PI;
+ return diffuse_radiance;
+}
+
const char*
get_copyright_notice(void* data)
{
diff --git a/src/stardis_smeteo.h b/src/stardis_smeteo.h
@@ -161,6 +161,35 @@ stardis_t_range
double range[2]); /* >0 [K] */
/*******************************************************************************
+ * Functions used to treat solar flux not as a flux imposed on the ground
+ * surface, but as a flux coming from an external source, i.e., the sun.
+ ******************************************************************************/
+/* Position of the sun based on the time, as well as the location and
+ * time zone associated with the meteorological data */
+STARDIS_API double*
+stardis_spherical_source_position
+ (const double time, /* [s] */
+ double position[3], /* [m] */
+ void* data);
+
+/* Solar power at a given time. It is computed from the direct component of
+ * the shortwave downward flux provided by meteorological data */
+STARDIS_API double /* [W] */
+stardis_spherical_source_power
+ (const double time,
+ void* data);
+
+/* Describes the diffuse part of the sun's contribution at a given time, i.e.
+ * the radiance emitted by the sun and scattered at least once in the
+ * atmosphere. It is computed from the diffuse component of the shortwave
+ * downward flux provided by meteorological data */
+STARDIS_API double /* [W/m^2/sr] */
+stardis_spherical_source_diffuse_radiance
+ (const double time,
+ const double dir[3],
+ void* data);
+
+/*******************************************************************************
* Various mandatory legal functions
******************************************************************************/
STARDIS_API const char*
diff --git a/src/stardis_smeteo_library.c b/src/stardis_smeteo_library.c
@@ -33,6 +33,9 @@ struct stardis_smeteo_lib {
double Tsrf_range[2]; /* Range of the surface temperatures [K] */
double Trad_range[2]; /* Range of the radiative temperatures [K] */
+ /* Algorithm for computing the solar position */
+ enum scem_sun_algo algo;
+
/* Number of seconds elapsed since the epoch until January 1, 1850, local time
* at the location of the smeteo file. Since the epoch is defined in
* UTC+00:00, this number of seconds is also defined in UTC+00:00. The
@@ -132,6 +135,10 @@ setup_smeteo
if((res = setup_utc_reference(lib)) != RES_OK) goto error;
+ /* TODO Add an argument to the library creation function to allow the caller
+ * to define the algorithm to be used */
+ lib->algo = SCEM_SUN_MEEUS;
+
/* Retrieve the maximum convection coefficient from meteorological data */
FOR_EACH(i, 0, desc.nentries) max_H = MMAX(desc.entries[i].H, max_H);
lib->max_convection_coef = max_H;
@@ -243,4 +250,6 @@ stardis_smeteo_lib_get_desc
desc->Tsrf_range[1] = lib->Tsrf_range[1];
desc->Trad_range[0] = lib->Trad_range[0];
desc->Trad_range[1] = lib->Trad_range[1];
+ desc->jan_1_1850 = lib->jan_1_1850;
+ desc->algo = lib->algo;
}
diff --git a/src/stardis_smeteo_library.h b/src/stardis_smeteo_library.h
@@ -16,6 +16,7 @@
#ifndef STARDIS_SMETEO_LIBRARY_H
#define STARDIS_SMETEO_LIBRARY_H
+#include <star/scem.h>
#include <rsys/rsys.h>
struct stardis_smeteo_lib_desc {
@@ -24,6 +25,17 @@ struct stardis_smeteo_lib_desc {
double max_convection_coef;
double Tsrf_range[2]; /* Range of surface temperatures [K] */
double Trad_range[2]; /* Range of radiative temperatures [K] */
+
+ /* Algorithm for computing the solar position */
+ enum scem_sun_algo algo;
+
+ /* Number of seconds elapsed since the epoch until January 1, 1850, local time
+ * at the location of the smeteo file. Since the epoch is defined in
+ * UTC+00:00, this number of seconds is also defined in UTC+00:00. The
+ * day_1850 field in the smeteo file can therefore be used to calculate the
+ * number of seconds to add to this member variable in order to convert the
+ * smeteo time to UTC+00:00 */
+ time_t jan_1_1850;
};
#define STARDIS_SMETEO_LIB_DESC_NULL__ {0}
static const struct stardis_smeteo_lib_desc STARDIS_SMETEO_LIB_DESC_NULL =
