commit 4063e991dac98879f7af669cac8807a266cfb914
parent 013e2e967e1e1dd72cb8b369e8bce3cc8799433c
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 19 Jun 2026 11:48:16 +0200
Merge branch 'release_0.1'
Diffstat:
10 files changed, 197 insertions(+), 108 deletions(-)
diff --git a/Makefile b/Makefile
@@ -22,7 +22,7 @@ LIBNAME_SHARED = libscem.so
LIBNAME = $(LIBNAME_$(LIB_TYPE))
default: library
-all: default tests util
+all: default tests util scem
################################################################################
# Library
@@ -75,7 +75,7 @@ UTIL_DEP = $(UTIL_SRC:.c=.d)
PKG_CONFIG_LOCAL = PKG_CONFIG_PATH="./:$${PKG_CONFIG_PATH}" $(PKG_CONFIG)
INCS_UTIL = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys scem-local)
-LIBS_UTIL = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys scem-local)
+LIBS_UTIL = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys scem-local) -lm
CFLAGS_UTIL = $(CFLAGS_EXE) $(INCS_UTIL)
LDFLAGS_UTIL = $(LDFLAGS_EXE) $(LIBS_UTIL)
diff --git a/README.md b/README.md
@@ -19,6 +19,15 @@ Edit config.mk as needed, then run:
## Release notes
+### Version 0.1
+
+- Rename zenith in elevation in man and code.
+ This breaks the API as the scem_sun_pos struct is modified.
+ However, the code still compute the same value: this is only a rename.
+- Add an utility function to convert a scem_sun_pos to a 3D vector giving the
+ direction from the sun to the observer.
+- Man wording improvement.
+
### Version 0.0
- Provide the scem library, which calculates the position of the sun at
diff --git a/config.mk b/config.mk
@@ -1,4 +1,8 @@
-VERSION = 0.0.0
+VERSION_MAJOR = 0
+VERSION_MINOR = 1
+VERSION_PATCH = 0
+VERSION = $(VERSION_MAJOR).$(VERSION_MINOR).$(VERSION_PATCH)
+
PREFIX = /usr/local
LIB_TYPE = SHARED
diff --git a/doc/scem.1 b/doc/scem.1
@@ -23,7 +23,7 @@
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
.Sh SYNOPSIS
.Nm
-.Op Fl hvr
+.Op Fl DhVvr
.Op Fl a Ar algo
.Op Fl d Ar utc_date
.Ar latitude
@@ -34,36 +34,43 @@
computes the position of the sun relative to a given point on the Earth's
surface, on a given date.
.Pp
-The longitude of the Earth's position is between
-.Bq -180, 180
-decimal degrees relative to Greenwich.
-The angle is positive towards the east.
-The latitude of the Earth's position is between
-.Bq -90, 90
-decimal degrees relative to the equator.
+The latitude must be in [-90, 90] degrees relative to the equator.
It is positive towards the north.
+The longitude must be in [-180, 180] degrees relative to Greenwich.
+The angle is positive towards the east.
.Pp
-On output,
-.Nm
-displays the position of the sun in zenith and azimuth on the standard
-output.
-The message is formatted as follows:
+The output depends on two options,
+.Fl D
+and
+.Fl V .
+If
+.Fl D
+is not provided,
+.No
+outputs the position of the sun in elevation and azimuth to the standard
+output with the following format:
.Bd -literal -offset Ds
-"%f %f\\n", zenith, azimuth
+"%f %f\\n", elevation, azimuth
.Ed
.Pp
-The zenith angle is between
-.Bq -90, 90
-degrees.
+The elevation angle is in [-90, 90] degrees.
It defines the elevation of the sun relative to the horizon.
-The azimuth angle ranges from
-.Bq 0, 360
-degrees, with 0° indicating north,
+The azimuth angle is in [0, 360] degrees, with 0° indicating north,
90° east,
180° south,
and 270° west.
.Pp
-The options are as follows
+If
+.Fl V
+is provided,
+.Nm
+outputs the direction from the sun as a 3D vector to the standard output
+with the following format:
+.Bd -literal -offset Ds
+"%f %f %f\\n", vec[0], vec[1], vec[2]
+.Ed
+.Pp
+The options are as follows:
.Bl -tag -width Ds
.\""""""""""""""""""""""""""""""""""""""
.It Fl a Ar algo
@@ -75,13 +82,15 @@ algorithm.
The solar position algorithms are as follows:
.Bl -tag -width Ds
.It Cm meeus
-The algorithm described in the Jean Meeus' book,
-.Dq Astronomical Algorithm .
+The algorithm described in Jean Meeus' book,
+.Dq Astronomical Algorithms .
.It Cm psa
The algorithm developed by the
.Dq Plataforma Solar de Almería .
.El
.\""""""""""""""""""""""""""""""""""""""
+.It Fl D
+Suppress the sun direction from the output.
.It Fl d Ar utc_date
Date on which the position of the sun is calculated at the specified
location.
@@ -99,29 +108,37 @@ date -u +"%Y-%m-%dT%H:%M:%S"
.It Fl h
Output short help and exit.
.\""""""""""""""""""""""""""""""""""""""
+.It Fl V
+Add the 3D vector from the sun to the observer to the output with the first,
+second and third components pointing east, north and up, respectively.
.It Fl v
Make
.Nm
verbose.
.\""""""""""""""""""""""""""""""""""""""
.It Fl r
-The solar output coordinates are defined in radians rather than degrees.
+Output the solar direction in radians rather than degrees.
+Meaningless if
+.Fl D
+is defined.
.El
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-.Sh EXIT
+.Sh EXIT STATUS
.Ex -std
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
.Sh EXAMPLES
-Calculate the position of the sun on the current date for a latitude of
-43.559962° east and a longitude of 1.468150° north:
+Compute the current position of the sun for a latitude of 43.559962° north and a
+longitude of 1.468150° east.
+Output both the direction in radians and the corresponding 3D vector:
.Bd -literal -offset Ds
-scem 43.559962 1.468150
+scem -V -r 43.559962 1.468150
.Ed
.Pp
-Same as above, but setting the observation date to October 21, 2015, at
-7:28:00 a.m.:
+Same location as above, setting the observation date to October 21, 2015, at
+11:28:00 AM UTC+2.
+Output the 3D vector only:
.Bd -literal -offset Ds
-scem -d 2015-10-21T7:28:00 43.559962 1.468150
+scem -D -V -d 2015-10-21T09:28:00 43.559962 1.468150
.Ed
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
.Sh SEE ALSO
@@ -131,14 +148,10 @@ scem -d 2015-10-21T7:28:00 43.559962 1.468150
.%D 1991
.Re
.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
+.%A Manuel Blanco
+.%A Kypros Milidonis
+.%A Aristides Bonanos
+.%T PSA+ updated PSA sun position algorithm
+.%U https://github.com/CST-Modelling-Tools/Updated-PSA-sun-position-algorithm
+.%D 2025
.Re
diff --git a/src/scem.c b/src/scem.c
@@ -16,6 +16,8 @@
#include "scem.h"
#include "scem_c.h"
+#include <math.h>
+
/*******************************************************************************
* Exported functions
******************************************************************************/
@@ -48,3 +50,31 @@ scem_sun_position_from_earth
}
return res;
}
+
+res_T
+scem_sun_position_to_sun_vector
+ (const struct scem_sun_pos* spherical,
+ double sun_dir[3])
+{
+ double cos_azimuth;
+ double sin_azimuth;
+ double cos_elevation;
+ double sin_elevation;
+ double azimuth;
+ double elevation;
+
+ if(!spherical || !sun_dir) return RES_BAD_ARG;
+
+ azimuth = spherical->azimuth;
+ elevation = spherical->elevation;
+
+ cos_azimuth = cos(azimuth);
+ sin_azimuth = sin(azimuth);
+ cos_elevation = cos(elevation);
+ sin_elevation = sin(elevation);
+
+ sun_dir[0] = -(cos_elevation * sin_azimuth);
+ sun_dir[1] = -(cos_elevation * cos_azimuth);
+ sun_dir[2] = -(sin_elevation);
+ return RES_OK;
+}
diff --git a/src/scem.h b/src/scem.h
@@ -39,7 +39,7 @@ enum scem_sun_algo {
/* Sun position */
struct scem_sun_pos {
/* In [-PI/2, PI/2]. Positive toward the sky */
- double zenith; /* [rad] */
+ double elevation; /* [rad] */
/* In [0, 2PI] with 0 aligned with north, PI/2 east, PI south, and 3PI/2 west */
double azimuth; /* [rad] */
@@ -73,10 +73,10 @@ scem_sun_algo_cstr(const enum scem_sun_algo algo)
BEGIN_DECLS
-/* Compute the position of the Sun in local coordinates (zenith and azimuth) for
- * a given observer at the surface of the Earth (can be extended to other
- * planets if required) i.e. the observer has to be located at the surface for
- * the required date */
+/* For a given location and time at the surface of the Earth (the observer has
+ * to be located at the surface), compute the position of the sun in local
+ * coordinates (elevation and azimuth) in radians.
+ * Could be extended to other planets if required. */
SCEM_API res_T
scem_sun_position_from_earth
(const struct tm* time, /* In UTC */
@@ -84,6 +84,14 @@ scem_sun_position_from_earth
const enum scem_sun_algo algo,
struct scem_sun_pos* sun);
+/* Convert the position of the sun (elevation and azimuth) in radians, as
+ * produced by scem_sun_position_from_earth, to a 3D vector (sun to observer).
+ * First component east, second component north, third component up. */
+SCEM_API res_T
+scem_sun_position_to_sun_vector
+ (const struct scem_sun_pos* sun,
+ double sun_vector[3]);
+
END_DECLS
#endif /* SCEM_H */
diff --git a/src/scem_main.c b/src/scem_main.c
@@ -18,6 +18,7 @@
#include "scem.h"
+#include <rsys/rsys.h>
#include <rsys/cstr.h>
#include <rsys/mem_allocator.h>
@@ -33,9 +34,11 @@ struct args {
int radian; /* Sun position in radians instead of degrees */
int verbose;
int quit;
+ int no_dir_output;
+ int vec_output;
};
static const struct args ARGS_DEFAULT = {
- SCEM_LOCATION_NULL__, {0}, SCEM_SUN_MEEUS, 0, 0, 0
+ SCEM_LOCATION_NULL__, {0}, SCEM_SUN_MEEUS, 0, 0, 0, 0, 0
};
/*******************************************************************************
@@ -45,7 +48,7 @@ static INLINE void
usage(FILE* stream)
{
fprintf(stream,
- "usage: scem [-hr] [-a algo] [-d utc_date] latitude longitude\n");
+ "usage: scem [-DhrVv] [-a algo] [-d utc_date] latitude longitude\n");
}
static res_T
@@ -66,9 +69,11 @@ parse_algo(const char* str, enum scem_sun_algo* algo)
static res_T
parse_date(const char* str, struct tm* time)
{
+ char *p;
ASSERT(str && time);
- if(!strptime(str, "%Y-%m-%dT%H:%M:%S\n", time)) {
+ p = strptime(str, "%Y-%m-%dT%H:%M:%S\n", time);
+ if(!p || *p != '\0') { /* strptime failed or trailing chars */
return RES_BAD_ARG;
} else {
return RES_OK;
@@ -87,7 +92,10 @@ parse_dbl
ASSERT(str && out_dbl && min <= max);
if((res = cstr_to_double(str, &dbl)) != RES_OK) return res;
- if(dbl < min || dbl > max) return RES_BAD_ARG;
+ if(dbl < min || dbl > max) {
+ fprintf(stderr, "scem: value not in range: %g [%g %g]\n", dbl, min, max);
+ return RES_BAD_ARG;
+ }
*out_dbl = dbl;
return RES_OK;
@@ -119,11 +127,14 @@ args_init(struct args* args, int argc, char** argv)
if((res = get_current_date(&args->time)) != RES_OK) goto error;
- while((opt = getopt(argc, argv, "a:d:hvr")) != -1) {
+ /* Don't process the last 2 args as options, they are lat and long */
+ while((opt = getopt(argc - 2, argv, "a:Dd:hVvr")) != -1) {
switch(opt) {
case 'a': res = parse_algo(optarg, &args->algo); break;
+ case 'D': args->no_dir_output = 1; break;
case 'd': res = parse_date(optarg, &args->time); break;
case 'h': usage(stdout); args->quit = 1; goto exit;
+ case 'V': args->vec_output = 1; break;
case 'v': args->verbose = 1; break;
case 'r': args->radian = 1; break;
default: res = RES_BAD_ARG; break;
@@ -159,6 +170,7 @@ main(int argc, char** argv)
{
struct args args = ARGS_DEFAULT;
struct scem_sun_pos pos = SCEM_SUN_POS_NULL;
+ double sun_dir[3];
int err = 0;
res_T res = RES_OK;
@@ -174,10 +186,17 @@ main(int argc, char** argv)
goto error;
}
- if(args.radian) {
- printf("%g %g\n", pos.zenith, pos.azimuth);
- } else {
- printf("%g %g\n", MRAD2DEG(pos.zenith), MRAD2DEG(pos.azimuth));
+ if(!args.no_dir_output) {
+ if(args.radian) {
+ printf("%g %g\n", pos.elevation, pos.azimuth);
+ } else {
+ printf("%g %g\n", MRAD2DEG(pos.elevation), MRAD2DEG(pos.azimuth));
+ }
+ }
+ if(args.vec_output) {
+ res = scem_sun_position_to_sun_vector(&pos, sun_dir);
+ if(res != RES_OK) goto error;
+ fprintf(stderr, "%g %g %g\n", SPLIT3(sun_dir));
}
exit:
diff --git a/src/scem_meeus.c b/src/scem_meeus.c
@@ -202,7 +202,7 @@ sun_position_from_earth_meeus
/* Zenital angle: angle between the horizontal plane and the Sun direction */
latitude_rad = MDEG2RAD(pos->latitude /* [deg] */);
- sun->zenith /* [rad] */ = asin
+ sun->elevation /* [rad] */ = asin
( sin(latitude_rad)*sin(celestial.declination)
+ cos(latitude_rad)*cos(celestial.declination)*cos(H_rad));
diff --git a/src/scem_psa.c b/src/scem_psa.c
@@ -26,11 +26,11 @@
static double
julian_day(const struct tm* time) /* In UTC */
{
- double h = 0; /* Decimal hours */
- int month = 0; /* Month in [1, 12] */
- int year = 0;
- long aux1 = 0;
- long aux2 = 0;
+ double h; /* Decimal hours */
+ int month; /* Month in [1, 12] */
+ int year;
+ long aux1;
+ long aux2;
ASSERT(time);
@@ -42,11 +42,11 @@ julian_day(const struct tm* time) /* In UTC */
aux1 = (month - 14) / 12;
aux2 = (1461 * (year + 4800 + aux1)) / 4
- + (367 * (month - 2 - 12*aux1)) / 12
- - (3 * ((year + 4900 + aux1)/100)) / 4
+ + (367 * (month - 2 - 12 * aux1)) / 12
+ - (3 * ((year + 4900 + aux1) / 100)) / 4
+ time->tm_mday - 32075;
- return (double)aux2 - 0.5 + h/24.0;
+ return (double)aux2 - 0.5 + h / 24.0;
}
/* Calculate ecliptic coordinates (ecliptic longitude and obliquity of the
@@ -55,26 +55,26 @@ julian_day(const struct tm* time) /* In UTC */
static void
sun_pos_ecliptic_coords(const double jd, struct ecliptic_coords* ecliptic)
{
- double mean_lon_rad=0;
- double mean_anomaly=0;
- double omega=0;
- double t=0;
+ double mean_lon_rad;
+ double mean_anomaly;
+ double omega;
+ double t;
ASSERT(jd >= 0 && ecliptic);
t = (jd - 2451545.0);
- omega = 2.1429 - 0.0010394594*t;
- mean_lon_rad = 4.8950630 + t*0.017202791698;
- mean_anomaly = 6.2400600 + t*0.0172019699;
+ omega = 2.267127827 - 0.00093003392670 * t;
+ mean_lon_rad = 4.895036035 + 0.01720279602 * t;
+ mean_anomaly = 6.239468336 + 0.01720200135 * t;
ecliptic->longitude = /* [rad] */
mean_lon_rad
- + 0.03341607*sin(mean_anomaly)
- + 0.00034894*sin(mean_anomaly*2) - 0.0001134
- - 0.0000203 *sin(omega);
+ + 0.03338320972 * sin(mean_anomaly)
+ + 0.0003497596876 * sin(mean_anomaly * 2) - 0.0001544353226
+ - 0.00000868972936 * sin(omega);
ecliptic->obliquity = /* [rad] */
- 0.4090928 - 6.2140e-9*t + 0.0000396*cos(omega);
+ 0.4090904909 - 6.213605399e-9 * t + 0.00004418094944 * cos(omega);
}
/* Calculate celestial coordinates (right ascension and declination) in radians
@@ -85,9 +85,9 @@ sun_pos_celestial_coords
(const struct ecliptic_coords* ecliptic,
struct celestial_coords* celestial)
{
- double sin_ecliptic_lon=0; /* [rad] */
- double y=0;
- double x=0;
+ double sin_ecliptic_lon; /* [rad] */
+ double y;
+ double x;
sin_ecliptic_lon = sin(ecliptic->longitude);
y = cos(ecliptic->obliquity) * sin_ecliptic_lon;
@@ -106,35 +106,34 @@ sun_pos_local_coords
const struct celestial_coords* celestial,
struct scem_sun_pos* sun)
{
- const double earth_radius = 6371.01; /* [km] */
- const double astro_unit = 149597890.0; /* [km] */
-
- double greenwitch_mean_side_real_time = 0;
- double local_mean_side_real_time = 0;
- double hour_angle = 0;
- double lat_rad = 0;
- double cos_lat = 0;
- double sin_lat = 0;
- double cos_hour_angle = 0;
- double azimuth_rad = 0;
- double zenith_rad = 0;
- double parallax = 0;
- double x = 0;
- double y = 0;
- double t = 0;
-
- double h = 0;
+ const double earth_radius = 6371.01; /* km */
+ const double astro_unit = 149597870.7; /* km */
+
+ double greenwich_mean_side_real_time;
+ double local_mean_side_real_time;
+ double hour_angle;
+ double lat_rad;
+ double cos_lat;
+ double sin_lat;
+ double cos_hour_angle;
+ double azimuth_rad;
+ double zenith_rad;
+ double parallax;
+ double x;
+ double y;
+ double t;
+ double h;
ASSERT(time && pos && celestial && sun);
t = (jd - 2451545.0);
h = (double)time->tm_hour
- + ((double)time->tm_min + (double)time->tm_sec/60.0) / 60.0;
+ + ((double)time->tm_min + (double)time->tm_sec / 60.0) / 60.0;
- greenwitch_mean_side_real_time = 6.6974243242 + 0.0657098283*t + h;
+ greenwich_mean_side_real_time = 6.697096103 + 0.06570984737 * t + h;
local_mean_side_real_time =
- MDEG2RAD(greenwitch_mean_side_real_time*15 + pos->longitude);
+ MDEG2RAD(greenwich_mean_side_real_time * 15 + pos->longitude);
hour_angle = local_mean_side_real_time - celestial->right_ascension;
lat_rad = MDEG2RAD(pos->latitude);
@@ -143,20 +142,20 @@ sun_pos_local_coords
cos_hour_angle = cos(hour_angle);
zenith_rad = acos
- ( cos_lat * cos(celestial->declination) * cos_hour_angle
+ (cos_lat * cos(celestial->declination) * cos_hour_angle
+ sin_lat * sin(celestial->declination));
y = -sin(hour_angle);
x = tan(celestial->declination) * cos_lat - sin_lat*cos_hour_angle;
azimuth_rad = atan2(y, x);
- if(azimuth_rad < 0.0) azimuth_rad += 2*PI;
+ if(azimuth_rad < 0.0) azimuth_rad += 2 * PI;
parallax = (earth_radius / astro_unit) * sin(zenith_rad);
zenith_rad += parallax;
- sun->azimuth = fmod(azimuth_rad, 2*PI);
- sun->zenith = PI*0.5 - fmod(zenith_rad, 2*PI);
+ sun->azimuth = fmod(azimuth_rad, 2 * PI);
+ sun->elevation = PI * 0.5 - fmod(zenith_rad, 2 * PI);
}
/*******************************************************************************
@@ -170,7 +169,7 @@ sun_position_from_earth_psa
{
struct ecliptic_coords ecliptic = ECLIPTIC_COORDS_NULL;
struct celestial_coords celestial = CELESTIAL_COORDS_NULL;
- double jd = 0; /* Julian Day */
+ double jd; /* Julian Day */
ASSERT(time && pos && sun);
diff --git a/src/test_scem_sun_position.c b/src/test_scem_sun_position.c
@@ -33,6 +33,7 @@ test_api(void)
struct tm utc = {0};
struct scem_location pos = SCEM_LOCATION_NULL;
struct scem_sun_pos sun = SCEM_SUN_POS_NULL;
+ double vec[3];
utc.tm_sec = 0;
utc.tm_min = 28;
@@ -47,6 +48,10 @@ test_api(void)
pos.longitude = 1.468150;
printf("Lat: %g°; lon: %g°\n", pos.latitude, pos.longitude);
+ CHK(scem_sun_position_to_sun_vector(NULL, NULL) == RES_BAD_ARG);
+ CHK(scem_sun_position_to_sun_vector(&sun, NULL) == RES_BAD_ARG);
+ CHK(scem_sun_position_to_sun_vector(NULL, vec) == RES_BAD_ARG);
+
CHK(scem_sun_position_from_earth(NULL, &pos, algo, &sun) == RES_BAD_ARG);
CHK(scem_sun_position_from_earth(&utc, NULL, algo, &sun) == RES_BAD_ARG);
CHK(scem_sun_position_from_earth(&utc, &pos, SCEM_SUN_ALGO_NONE__, &sun)
@@ -69,8 +74,10 @@ test_api(void)
CHK(scem_sun_position_from_earth(&utc, &pos, algo, &sun) == RES_OK);
printf("%s - zenith: %g°; azimuth: %g°\n",
scem_sun_algo_cstr(algo),
- MRAD2DEG(sun.zenith),
+ MRAD2DEG(sun.elevation),
MRAD2DEG(sun.azimuth));
+ CHK(scem_sun_position_to_sun_vector(&sun, vec) == RES_OK);
+ printf("Vector: %g %g %g\n", SPLIT3(vec));
}
}
