commit c640fcdcf2bd40bd0f1c361b50636dffcb816534
parent 6abf1f34d9a3c22ced471e13c5e5a6c001b6a8bc
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 7 Oct 2025 15:36:16 +0200
Add the PSA algorithm as an alternative
This was a near-literal translation of the PSA files from
http://www.psa.es/sdg/sunpos.htm; only the layout was modified to
conform to the coding style.
The position of the sun can thus be calculated using two different
algorithms: the one from Jean Meeus' book, "Astronomical Algorithms",
and the PSA algorithm. The caller chooses which algorithm to use when
invoking the function.
Please note that the original sources of the PSA algorithm contain no
licensing information or copyright notice. They are simply assumed to be
distributed under a free license, with no guarantee that this is the
case. Please contact us if this is not true and the code cannot be
distributed under the GPL license.
Diffstat:
| M | Makefile | | | 2 | +- |
| M | src/scem.c | | | 218 | +++++++------------------------------------------------------------------------ |
| M | src/scem.h | | | 43 | ++++++++++++++++++++++++++++++++----------- |
| A | src/scem_c.h | | | 46 | ++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/scem_meeus.c | | | 217 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/scem_psa.c | | | 183 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | src/test_scem_sun_position.c | | | 33 | ++++++++++++++++++++------------- |
7 files changed, 518 insertions(+), 224 deletions(-)
diff --git a/Makefile b/Makefile
@@ -27,7 +27,7 @@ all: default tests
################################################################################
# Library
################################################################################
-SRC = src/scem.c
+SRC = src/scem.c src/scem_meeus.c src/scem_psa.c
OBJ = $(SRC:.c=.o)
DEP = $(SRC:.c=.d)
diff --git a/src/scem.c b/src/scem.c
@@ -14,166 +14,7 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "scem.h"
-
-#include <rsys/math.h>
-#include <rsys/rsys.h>
-
-#include <math.h>
-#include <time.h>
-
-struct radr { /* right-ascension, declination, radius */
- double right_ascension; /* [rad] */
- double declination; /* [rad] */
- double radius; /* [m] */
-};
-static const struct radr RADR_NULL = {0};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static double
-julian_day(const struct tm* time) /* In UTC */
-{
- double jd = 0; /* julian day */
- double d = 0; /* Decimal day */
- int m = 0; /* Month */
- int y = 0; /* Day */
- int a, b;
-
- ASSERT(time);
-
- d = (double)time->tm_mday
- + (double)time->tm_hour/24.0
- + (double)time->tm_min/1440.0
- + (double)time->tm_sec/86400.0;
- m = time->tm_mon + 1; /* The months start at 0 */
- y = time->tm_year + 1900; /* tm_year is relative to year 1900 (see tm(3type)) */
-
- if((unsigned)time->tm_mon < 3) {
- y -= 1;
- m += 12;
- }
-
- a = y / 100;
- b = 2 - a + a/4;
-
- /* Found in "Astronomical Algorithms" by Jean Meeus */
- jd = (int)(365.25 * (double)(y+4715))
- + (int)(30.6001 * (double)(m+1))
- + d + b - 1524.5;
-
- return jd;
-}
-
-/* degrees minutes seconds to degrees conversion */
-static FINLINE double
-dms2deg(int d, int m, double s)
-{
- return (double)d + (double)m/60.0 + s/3600.0;
-}
-
-static void
-sun_pos_equatorial_coords(const double jd, struct radr* radr)
-{
- double e=0; /* excentricity of Earth's orbit */
- double t=0; /* time in Julian centuries */
- double v_rad=0; /* True anomaly of the sun */
- double R=0; /* Sun radisu vector */
- double d1=0, d2=0, d3=0, d4=0;
- double L0_deg=0;
- double L_deg=0, L_rad=0;
- double Lprime_deg=0, Lprime_rad=0;
- double M_deg=0, M_rad=0;
- double C_deg=0, C_rad=0;
- double alpha_rad=0;
- double delta_epsilon_deg=0;
- double delta_rad=0;
- double epsilon0_deg=0;
- double epsilon_deg=0, epsilon_rad=0;
- double lambda_deg=0, lambda_rad=0;
- double omega_deg=0, omega_rad=0;
- double theta_deg=0;
- ASSERT(jd > 0 && radr);
-
- t = (jd - 2451545.0)/36525.0; /* time in julian centuries of 36525 days */
-
- /* L0: geometric mean longitude of the Sun,
- * relative to the mean equinox of */
- L0_deg = 280.46645 + 36000.76983*t + 3.023e-4*t*t;
- L0_deg = fmod(L0_deg, 360);
- /* M: mean anomaly of the Sun */
- M_deg = 357.52910 + 35999.05030*t - 1.559e-4*t*t - 4.80e-7*t*t*t;
- M_deg = fmod(M_deg, 360);
- M_rad = MDEG2RAD(M_deg);
-
- /* e: excentricity of Earth's orbit */
- e = 1.6708617e-2 - 4.2037e-5*t - 1.236e-7*t*t;
-
- /* equation of the center of the Sun */
- C_deg = sin(1.0*M_rad) * (1.9146e-0 - 4.817e-3*t - 1.40e-5*t*t)
- + sin(2.0*M_rad) * (1.9993e-2 - 1.010e-4*t)
- + sin(3.0*M_rad) * (2.9000e-4);
- C_rad = MDEG2RAD(C_deg);
-
- /* theta: true longitude of the Sun */
- theta_deg = L0_deg + C_deg;
-
- /* v: true anomaly of the Sun */
- v_rad = M_rad + C_rad;
-
- /* Sun radius vector */
- R = 1.000001018*(1.0-e*e)/(1.0 + e*cos(v_rad)); /* [a.u] */
-
- /* lambda: apparent longitude of the Sun,
- * relative to the true equinox of the date */
- omega_deg = 125.04 - 1934.136*t;
- omega_rad = MDEG2RAD(omega_deg);
- lambda_deg = theta_deg - 5.69e-3 - 4.78e-3*sin(omega_rad);
- lambda_rad = MDEG2RAD(lambda_deg);
-
- /* epsilon0: mean obliquity of the ecliptic */
- d1 = dms2deg(23, 26, 21.448);
- d2 = dms2deg( 0, 0, 46.8150);
- d3 = dms2deg( 0, 0, 5.90e-4);
- d4 = dms2deg( 0, 0, 1.813e-3);
- epsilon0_deg = d1 - d2*t - d3*t*t + d4*t*t*t;
-
- /* L: mean longitude of the Sun */
- L_deg = 280.4665 + 36000.7698*t;
- L_deg = fmod(L_deg, 360);
- L_rad = MDEG2RAD(L_deg);
-
- /* Lprime: mean longitude of the Moon */
- Lprime_deg = 218.3165 + 481267.8813*t;
- Lprime_deg = fmod(Lprime_deg, 360);
- Lprime_rad = MDEG2RAD(Lprime_deg);
-
- /* Delta_epsilon: nutation in obliquity */
- d1 = dms2deg(0, 0, 9.20);
- d2 = dms2deg(0, 0, 5.70e-1);
- d3 = dms2deg(0, 0, 1.0e-1);
- d4 = dms2deg(0, 0, 9.0e-2);
- delta_epsilon_deg =
- d1*cos(omega_rad)
- + d2*cos(2.0*L_rad)
- + d3*cos(2.0*Lprime_rad)
- - d4*cos(2.0*omega_rad);
-
- /* epsilon: obliquity of the ecliptic */
- epsilon_deg = epsilon0_deg + delta_epsilon_deg + 2.56e-3*cos(omega_rad);
- epsilon_rad = MDEG2RAD(epsilon_deg);
-
- /* alpha: right ascension of the Sun */
- alpha_rad = atan2(cos(epsilon_rad)*sin(lambda_rad), cos(lambda_rad));
-
- /* delta: declination of the Sun */
- delta_rad = asin(sin(epsilon_rad)*sin(lambda_rad));
-
- /* Setup output data */
- radr->right_ascension = alpha_rad;
- radr->declination = delta_rad;
- radr->radius = R;
-}
+#include "scem_c.h"
/*******************************************************************************
* Exported functions
@@ -181,44 +22,23 @@ sun_pos_equatorial_coords(const double jd, struct radr* radr)
res_T
scem_sun_position_from_earth
(const struct tm* time, /* In UTC */
- const struct scem_gcs* pos, /* Local position */
- struct scem_zad* sun) /* Sun position */
+ const struct scem_location* pos, /* Local position */
+ const enum scem_sun_algo algo,
+ struct scem_sun_pos* sun)
{
- struct radr radr = RADR_NULL;
- double jd = 0; /* Julian Day */
- double theta_deg = 0;
- double theta_rad = 0;
- double H_rad = 0;
- double latitude_rad = 0;
-
- if(!time || !pos || !sun) return RES_BAD_ARG;
-
- jd = julian_day(time);
- sun_pos_equatorial_coords(jd, &radr);
-
- /* theta: true solar time (taking into account the longitude of the observer) */
- theta_deg = 280.1470 + 360.9856235*(jd-2451545.0) + pos->longitude /* [deg] */;
- theta_deg = fmod(theta_deg, 360);
- theta_rad = MDEG2RAD(theta_deg);
-
- /* H: angular distance between the observer and the Sun */
- H_rad = theta_rad - radr.right_ascension;
-
- /* Zenital angle: angle between the horizontal plane and the Sun direction */
- latitude_rad = MDEG2RAD(pos->latitude /* [deg] */);
- sun->zenith /* [rad] */ = asin
- ( sin(latitude_rad)*sin(radr.declination)
- + cos(latitude_rad)*cos(radr.declination)*cos(H_rad));
-
- /* Azimutal angle: angle between the North and the projection of the Sun
- * direction on the horizontal plane, with a rotation direction from North to
- * East. */
- sun->azimuth /* [rad] */ = PI + atan2
- (sin(H_rad), cos(H_rad)*sin(latitude_rad) - tan(radr.declination)*cos(latitude_rad));
-
- /* Do not take into account the impact of the transformation in the local
- * reference frame on the calculation of the distance between the observer and
- * the sun */
- sun->distance = radr.radius;
- return RES_OK;
+ res_T res = RES_OK;
+
+ if(!time || !pos || !sun || (unsigned)algo >= SCEM_SUN_ALGO_NONE__)
+ return RES_BAD_ARG;
+
+ switch(algo) {
+ case SCEM_SUN_MEEUS:
+ res = sun_position_from_earth_meeus(time, pos, sun);
+ break;
+ case SCEM_SUN_PSA:
+ res = sun_position_from_earth_psa(time, pos, sun);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ return res;
}
diff --git a/src/scem.h b/src/scem.h
@@ -27,17 +27,25 @@
/* Forward declaration */
struct tm;
-/* Position in zenith, azimuth, distance */
-struct scem_zad {
- double zenith; /* [-PI/2, PI/2] */
+enum scem_sun_algo {
+ /* Algorithm taken from Jean Meeus' book "Astronomical Algorithm" */
+ SCEM_SUN_MEEUS,
+ /* Algorithm from the "Plataforma Solar de Almería" */
+ SCEM_SUN_PSA,
+
+ SCEM_SUN_ALGO_NONE__
+};
+
+/* Sun position */
+struct scem_sun_pos {
+ double zenith; /* [-PI/2, PI/2]. Positive toward the sky */
double azimuth; /* [0, 2/PI] */
- double distance; /* [m] */
};
-#define SCEM_ZAD_NULL__ {0, 0, 0}
-static const struct scem_zad SCEM_ZAD_NULL = SCEM_ZAD_NULL__;
+#define SCEM_SUN_POS_NULL__ {0, 0}
+static const struct scem_sun_pos SCEM_SUN_POS_NULL = SCEM_SUN_POS_NULL__;
/* Position defined in geographic coordinate system */
-struct scem_gcs {
+struct scem_location {
/* In [-180,180] decimal degree relative to Greenwitch.
* Positive toward the east */
double longitude;
@@ -45,8 +53,20 @@ struct scem_gcs {
* Positive toward the north */
double latitude;
};
-#define SCEM_GCS_NULL__ {0,0}
-static const struct scem_gcs SCEM_GCS_NULL = SCEM_GCS_NULL__;
+#define SCEM_LOCATION_NULL__ {0,0}
+static const struct scem_location SCEM_LOCATION_NULL = SCEM_LOCATION_NULL__;
+
+static INLINE const char*
+scem_sun_algo_cstr(const enum scem_sun_algo algo)
+{
+ const char* cstr = "None";
+ switch(algo) {
+ case SCEM_SUN_MEEUS: cstr = "Meeus"; break;
+ case SCEM_SUN_PSA: cstr = "PSA"; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ return cstr;
+}
BEGIN_DECLS
@@ -57,8 +77,9 @@ BEGIN_DECLS
SCEM_API res_T
scem_sun_position_from_earth
(const struct tm* time, /* In UTC */
- const struct scem_gcs* position, /* Local position */
- struct scem_zad* sun);
+ const struct scem_location* position, /* Local position */
+ const enum scem_sun_algo algo,
+ struct scem_sun_pos* sun);
END_DECLS
diff --git a/src/scem_c.h b/src/scem_c.h
@@ -0,0 +1,46 @@
+/* Copyright (C) 2025 |Méso|Star> (contact@meso-star.com)
+ *
+ * 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/>. */
+
+#ifndef SCEM_C_H
+#define SCEM_C_H
+
+#include "scem.h"
+#include <rsys/rsys.h>
+
+struct celestial_coords { /* right-ascension, declination */
+ double right_ascension; /* [rad] */
+ double declination; /* [rad] */
+};
+static const struct celestial_coords CELESTIAL_COORDS_NULL = {0};
+
+struct ecliptic_coords {
+ double longitude; /* [rad] */
+ double obliquity; /* [rad] */
+};
+static const struct ecliptic_coords ECLIPTIC_COORDS_NULL = {0};
+
+extern LOCAL_SYM res_T
+sun_position_from_earth_psa
+ (const struct tm* time, /* In UTC */
+ const struct scem_location* pos,
+ struct scem_sun_pos* sun);
+
+extern LOCAL_SYM res_T
+sun_position_from_earth_meeus
+ (const struct tm* time, /* In UTC */
+ const struct scem_location* pos,
+ struct scem_sun_pos* sun);
+
+#endif /* SCEM_C_H */
diff --git a/src/scem_meeus.c b/src/scem_meeus.c
@@ -0,0 +1,217 @@
+/* Copyright (C) 2025 |Méso|Star> (contact@meso-star.com)
+ *
+ * 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/>. */
+
+#include "scem_c.h"
+
+#include <rsys/math.h>
+#include <rsys/rsys.h>
+
+#include <math.h>
+#include <time.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE double
+julian_day(const struct tm* time) /* In UTC */
+{
+ double jd = 0; /* julian day */
+ double d = 0; /* Decimal day */
+ int m = 0; /* Month */
+ int y = 0; /* Day */
+ int a, b;
+
+ ASSERT(time);
+
+ d = (double)time->tm_mday
+ + (double)time->tm_hour/24.0
+ + (double)time->tm_min/1440.0
+ + (double)time->tm_sec/86400.0;
+ m = time->tm_mon + 1; /* The months start at 0 */
+ y = time->tm_year + 1900; /* tm_year is relative to year 1900 (see tm(3type)) */
+
+ if((unsigned)time->tm_mon < 3) {
+ y -= 1;
+ m += 12;
+ }
+
+ a = y / 100;
+ b = 2 - a + a/4;
+
+ /* Found in "Astronomical Algorithms" by Jean Meeus */
+ jd = (int)(365.25 * (double)(y+4715))
+ + (int)(30.6001 * (double)(m+1))
+ + d + b - 1524.5;
+
+ return jd;
+}
+
+/* degrees minutes seconds to degrees conversion */
+static FINLINE double
+dms2deg(int d, int m, double s)
+{
+ return (double)d + (double)m/60.0 + s/3600.0;
+}
+
+static void
+sun_pos_celestial_coords
+ (const double jd,
+ struct celestial_coords* celestial)
+{
+ double e=0; /* excentricity of Earth's orbit */
+ double t=0; /* time in Julian centuries */
+ double v_rad=0; /* True anomaly of the sun */
+ double R=0; /* Sun radius vector */
+ double d1=0, d2=0, d3=0, d4=0;
+ double L0_deg=0;
+ double L_deg=0, L_rad=0;
+ double Lprime_deg=0, Lprime_rad=0;
+ double M_deg=0, M_rad=0;
+ double C_deg=0, C_rad=0;
+ double alpha_rad=0;
+ double delta_epsilon_deg=0;
+ double delta_rad=0;
+ double epsilon0_deg=0;
+ double epsilon_deg=0, epsilon_rad=0;
+ double lambda_deg=0, lambda_rad=0;
+ double omega_deg=0, omega_rad=0;
+ double theta_deg=0;
+ ASSERT(jd > 0 && celestial);
+
+ t = (jd - 2451545.0)/36525.0; /* time in julian centuries of 36525 days */
+
+ /* L0: geometric mean longitude of the Sun,
+ * relative to the mean equinox of */
+ L0_deg = 280.46645 + 36000.76983*t + 3.023e-4*t*t;
+ L0_deg = fmod(L0_deg, 360);
+ /* M: mean anomaly of the Sun */
+ M_deg = 357.52910 + 35999.05030*t - 1.559e-4*t*t - 4.80e-7*t*t*t;
+ M_deg = fmod(M_deg, 360);
+ M_rad = MDEG2RAD(M_deg);
+
+ /* e: excentricity of Earth's orbit */
+ e = 1.6708617e-2 - 4.2037e-5*t - 1.236e-7*t*t;
+
+ /* equation of the center of the Sun */
+ C_deg = sin(1.0*M_rad) * (1.9146e-0 - 4.817e-3*t - 1.40e-5*t*t)
+ + sin(2.0*M_rad) * (1.9993e-2 - 1.010e-4*t)
+ + sin(3.0*M_rad) * (2.9000e-4);
+ C_rad = MDEG2RAD(C_deg);
+
+ /* theta: true longitude of the Sun */
+ theta_deg = L0_deg + C_deg;
+
+ /* v: true anomaly of the Sun */
+ v_rad = M_rad + C_rad;
+
+ /* Sun radius vector */
+ R = 1.000001018*(1.0-e*e)/(1.0 + e*cos(v_rad)); /* [a.u] */
+ (void)R; /* Not used */
+
+ /* lambda: apparent longitude of the Sun,
+ * relative to the true equinox of the date */
+ omega_deg = 125.04 - 1934.136*t;
+ omega_rad = MDEG2RAD(omega_deg);
+ lambda_deg = theta_deg - 5.69e-3 - 4.78e-3*sin(omega_rad);
+ lambda_rad = MDEG2RAD(lambda_deg);
+
+ /* epsilon0: mean obliquity of the ecliptic */
+ d1 = dms2deg(23, 26, 21.448);
+ d2 = dms2deg( 0, 0, 46.8150);
+ d3 = dms2deg( 0, 0, 5.90e-4);
+ d4 = dms2deg( 0, 0, 1.813e-3);
+ epsilon0_deg = d1 - d2*t - d3*t*t + d4*t*t*t;
+
+ /* L: mean longitude of the Sun */
+ L_deg = 280.4665 + 36000.7698*t;
+ L_deg = fmod(L_deg, 360);
+ L_rad = MDEG2RAD(L_deg);
+
+ /* Lprime: mean longitude of the Moon */
+ Lprime_deg = 218.3165 + 481267.8813*t;
+ Lprime_deg = fmod(Lprime_deg, 360);
+ Lprime_rad = MDEG2RAD(Lprime_deg);
+
+ /* Delta_epsilon: nutation in obliquity */
+ d1 = dms2deg(0, 0, 9.20);
+ d2 = dms2deg(0, 0, 5.70e-1);
+ d3 = dms2deg(0, 0, 1.0e-1);
+ d4 = dms2deg(0, 0, 9.0e-2);
+ delta_epsilon_deg =
+ d1*cos(omega_rad)
+ + d2*cos(2.0*L_rad)
+ + d3*cos(2.0*Lprime_rad)
+ - d4*cos(2.0*omega_rad);
+
+ /* epsilon: obliquity of the ecliptic */
+ epsilon_deg = epsilon0_deg + delta_epsilon_deg + 2.56e-3*cos(omega_rad);
+ epsilon_rad = MDEG2RAD(epsilon_deg);
+
+ /* alpha: right ascension of the Sun */
+ alpha_rad = atan2(cos(epsilon_rad)*sin(lambda_rad), cos(lambda_rad));
+
+ /* delta: declination of the Sun */
+ delta_rad = asin(sin(epsilon_rad)*sin(lambda_rad));
+
+ /* Setup output data */
+ celestial->right_ascension = alpha_rad;
+ celestial->declination = delta_rad;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+sun_position_from_earth_meeus
+ (const struct tm* time, /* In UTC */
+ const struct scem_location* pos, /* Local position */
+ struct scem_sun_pos* sun)
+{
+ struct celestial_coords celestial = CELESTIAL_COORDS_NULL;
+ double jd = 0; /* Julian Day */
+ double theta_deg = 0;
+ double theta_rad = 0;
+ double H_rad = 0;
+ double latitude_rad = 0;
+
+ ASSERT(time && pos && sun);
+
+ jd = julian_day(time);
+
+ sun_pos_celestial_coords(jd, &celestial);
+
+ /* theta: true solar time (taking into account the longitude of the observer) */
+ theta_deg = 280.1470 + 360.9856235*(jd-2451545.0) + pos->longitude /* [deg] */;
+ theta_deg = fmod(theta_deg, 360);
+ theta_rad = MDEG2RAD(theta_deg);
+
+ /* H: angular distance between the observer and the Sun */
+ H_rad = theta_rad - celestial.right_ascension;
+
+ /* Zenital angle: angle between the horizontal plane and the Sun direction */
+ latitude_rad = MDEG2RAD(pos->latitude /* [deg] */);
+ sun->zenith /* [rad] */ = asin
+ ( sin(latitude_rad)*sin(celestial.declination)
+ + cos(latitude_rad)*cos(celestial.declination)*cos(H_rad));
+
+ /* Azimutal angle: angle between the North and the projection of the Sun
+ * direction on the horizontal plane, with a rotation direction from North to
+ * East. */
+ sun->azimuth /* [rad] */ = PI + atan2
+ (sin(H_rad),
+ cos(H_rad)*sin(latitude_rad) - tan(celestial.declination)*cos(latitude_rad));
+
+ return RES_OK;
+}
diff --git a/src/scem_psa.c b/src/scem_psa.c
@@ -0,0 +1,183 @@
+/* Copyright (C) 2025 |Méso|Star> (contact@meso-star.com)
+ *
+ * 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/>. */
+
+#include "scem_c.h"
+
+#include <rsys/math.h>
+
+#include <math.h>
+#include <time.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+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;
+
+ ASSERT(time);
+
+ h = (double)time->tm_hour
+ + ((double)time->tm_min + (double)time->tm_sec/60.0) / 60.0;
+
+ month = time->tm_mon + 1; /* tm_mon start at 0 */
+ year = time->tm_year + 1900; /* tm_year is relative to year 1900 */
+
+ aux1 = (month - 14) / 12;
+ aux2 = (1461 * (year + 4800 + aux1)) / 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;
+}
+
+/* Calculate ecliptic coordinates (ecliptic longitude and obliquity of the
+ * ecliptic in radians but without limiting the angle to be less than 2*Pi
+ * (i.e., the result may be greater than 2*Pi) */
+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;
+
+ ASSERT(jd >= 2451545.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;
+
+ ecliptic->longitude = /* [rad] */
+ mean_lon_rad
+ + 0.03341607*sin(mean_anomaly)
+ + 0.00034894*sin(mean_anomaly*2) - 0.0001134
+ - 0.0000203 *sin(omega);
+
+ ecliptic->obliquity = /* [rad] */
+ 0.4090928 - 6.2140e-9*t + 0.0000396*cos(omega);
+}
+
+/* Calculate celestial coordinates (right ascension and declination) in radians
+ * * but without limiting the angle to be less than 2*Pi (i.e., the result may
+ * be greater than 2*Pi) */
+static void
+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;
+
+ sin_ecliptic_lon = sin(ecliptic->longitude);
+ y = cos(ecliptic->obliquity) * sin_ecliptic_lon;
+ x = cos(ecliptic->longitude);
+
+ celestial->right_ascension = atan2(y, x);
+ if(celestial->right_ascension < 0) celestial->right_ascension += 2*PI;
+ celestial->declination = asin(sin(ecliptic->obliquity) * sin_ecliptic_lon);
+}
+
+static void
+sun_pos_local_coords
+ (const struct tm* time, /* In UTC */
+ const struct scem_location* pos,
+ const double jd, /* Julian day */
+ 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;
+
+ 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;
+
+ greenwitch_mean_side_real_time = 6.6974243242 + 0.0657098283*t + h;
+
+ local_mean_side_real_time =
+ MDEG2RAD(greenwitch_mean_side_real_time*15 + pos->longitude);
+
+ hour_angle = local_mean_side_real_time - celestial->right_ascension;
+ lat_rad = MDEG2RAD(pos->latitude);
+ cos_lat = cos(lat_rad);
+ sin_lat = sin(lat_rad);
+ cos_hour_angle = cos(hour_angle);
+
+ zenith_rad = acos
+ ( 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;
+
+ 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);
+}
+
+/*******************************************************************************
+ * Local symbols
+ ******************************************************************************/
+res_T
+sun_position_from_earth_psa
+ (const struct tm* time, /* In UTC */
+ const struct scem_location* pos, /* Local position */
+ struct scem_sun_pos* sun) /* Sun position */
+{
+ struct ecliptic_coords ecliptic = ECLIPTIC_COORDS_NULL;
+ struct celestial_coords celestial = CELESTIAL_COORDS_NULL;
+ double jd = 0; /* Julian Day */
+
+ ASSERT(time && pos && sun);
+
+ jd = julian_day(time);
+ sun_pos_ecliptic_coords(jd, &ecliptic);
+ sun_pos_celestial_coords(&ecliptic, &celestial);
+ sun_pos_local_coords(time, pos, jd, &celestial, sun);
+
+ return RES_OK;
+}
diff --git a/src/test_scem_sun_position.c b/src/test_scem_sun_position.c
@@ -27,28 +27,35 @@ static void
test_api(void)
{
time_t epoch;
- struct tm* utc = NULL;
- struct scem_gcs pos = SCEM_GCS_NULL;
- struct scem_zad sun = SCEM_ZAD_NULL;
+ enum scem_sun_algo algo = SCEM_SUN_MEEUS;
+ struct tm utc = {0};
+ struct scem_location pos = SCEM_LOCATION_NULL;
+ struct scem_sun_pos sun = SCEM_SUN_POS_NULL;
epoch = time(NULL);
CHK(epoch != ((time_t)-1));
- utc = gmtime(&epoch);
- CHK(utc != NULL);
+ utc.tm_sec = 0;
+ utc.tm_min = 28;
+ utc.tm_hour = 17;
+ utc.tm_mday = 6;
+ utc.tm_mon = 9;
+ utc.tm_year = 125;
- printf("%s", asctime(utc));
+ printf("%s", asctime(&utc));
- pos.latitude = 43.604600;
+ pos.latitude = 43.60460;
pos.longitude = 1.44422;
printf("Lat: %g°; lon: %g°\n", pos.latitude, pos.longitude);
- CHK(scem_sun_position_from_earth(NULL, &pos, &sun) == RES_BAD_ARG);
- CHK(scem_sun_position_from_earth(utc, NULL, &sun) == RES_BAD_ARG);
- CHK(scem_sun_position_from_earth(utc, &pos, NULL) == RES_BAD_ARG);
- CHK(scem_sun_position_from_earth(utc, &pos, &sun) == RES_OK);
-
- printf("azimuth: %g°; elevation %g°\n",
+ 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)
+ == RES_BAD_ARG);
+ CHK(scem_sun_position_from_earth(&utc, &pos, algo, NULL) == RES_BAD_ARG);
+ CHK(scem_sun_position_from_earth(&utc, &pos, algo, &sun) == RES_OK);
+ printf("%s - azimuth: %g°; elevation: %g°\n",
+ scem_sun_algo_cstr(algo),
MRAD2DEG(sun.azimuth),
MRAD2DEG(sun.zenith));
}
