commit 4442b161fc9b3bba67374294469e2ccbaaa21a1e
parent d32534a583aca49f84f1a951af0d4c00820cee7e
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 4 Jun 2026 17:15:23 +0200
Update de PSA algorithm to the latest published
From the github official repository.
Diffstat:
| M | doc/scem.1 | | | 16 | ++++++---------- |
| M | src/scem_psa.c | | | 95 | +++++++++++++++++++++++++++++++++++++++---------------------------------------- |
2 files changed, 53 insertions(+), 58 deletions(-)
diff --git a/doc/scem.1 b/doc/scem.1
@@ -131,14 +131,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_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->zenith = 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);
