scem_meeus.c (6716B)
1 /* Copyright (C) 2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "scem_c.h" 17 18 #include <rsys/math.h> 19 #include <rsys/rsys.h> 20 21 #include <math.h> 22 #include <time.h> 23 24 /******************************************************************************* 25 * Helper functions 26 ******************************************************************************/ 27 static INLINE double 28 julian_day(const struct tm* time) /* In UTC */ 29 { 30 double jd = 0; /* julian day */ 31 double d = 0; /* Decimal day */ 32 int m = 0; /* Month */ 33 int y = 0; /* Day */ 34 int a, b; 35 36 ASSERT(time); 37 38 d = (double)time->tm_mday 39 + (double)time->tm_hour/24.0 40 + (double)time->tm_min/1440.0 41 + (double)time->tm_sec/86400.0; 42 m = time->tm_mon + 1; /* tm_mon start at 0 */ 43 y = time->tm_year + 1900; /* tm_year is relative to year 1900 (see tm(3type)) */ 44 45 if(m < 3) { 46 y -= 1; 47 m += 12; 48 } 49 50 a = y / 100; 51 b = 2 - a + a/4; 52 53 /* Found in "Astronomical Algorithms" by Jean Meeus */ 54 jd = (int)(365.25 * (double)(y+4716)) 55 + (int)(30.6001 * (double)(m+1)) 56 + d + b - 1524.5; 57 58 return jd; 59 } 60 61 /* degrees minutes seconds to degrees conversion */ 62 static FINLINE double 63 dms2deg(int d, int m, double s) 64 { 65 return (double)d + (double)m/60.0 + s/3600.0; 66 } 67 68 static void 69 sun_pos_celestial_coords 70 (const double jd, 71 struct celestial_coords* celestial) 72 { 73 double e=0; /* excentricity of Earth's orbit */ 74 double t=0; /* time in Julian centuries */ 75 double v_rad=0; /* True anomaly of the sun */ 76 double R=0; /* Sun radius vector */ 77 double d1=0, d2=0, d3=0, d4=0; 78 double L0_deg=0; 79 double L_deg=0, L_rad=0; 80 double Lprime_deg=0, Lprime_rad=0; 81 double M_deg=0, M_rad=0; 82 double C_deg=0, C_rad=0; 83 double alpha_rad=0; 84 double delta_epsilon_deg=0; 85 double delta_rad=0; 86 double epsilon0_deg=0; 87 double epsilon_deg=0, epsilon_rad=0; 88 double lambda_deg=0, lambda_rad=0; 89 double omega_deg=0, omega_rad=0; 90 double theta_deg=0; 91 ASSERT(jd > 0 && celestial); 92 93 t = (jd - 2451545.0)/36525.0; /* time in julian centuries of 36525 days */ 94 95 /* L0: geometric mean longitude of the Sun, 96 * relative to the mean equinox of */ 97 L0_deg = 280.46645 + 36000.76983*t + 3.023e-4*t*t; 98 L0_deg = fmod(L0_deg, 360); 99 /* M: mean anomaly of the Sun */ 100 M_deg = 357.52910 + 35999.05030*t - 1.559e-4*t*t - 4.80e-7*t*t*t; 101 M_deg = fmod(M_deg, 360); 102 M_rad = MDEG2RAD(M_deg); 103 104 /* e: excentricity of Earth's orbit */ 105 e = 1.6708617e-2 - 4.2037e-5*t - 1.236e-7*t*t; 106 107 /* equation of the center of the Sun */ 108 C_deg = sin(1.0*M_rad) * (1.9146e-0 - 4.817e-3*t - 1.40e-5*t*t) 109 + sin(2.0*M_rad) * (1.9993e-2 - 1.010e-4*t) 110 + sin(3.0*M_rad) * (2.9000e-4); 111 C_rad = MDEG2RAD(C_deg); 112 113 /* theta: true longitude of the Sun */ 114 theta_deg = L0_deg + C_deg; 115 116 /* v: true anomaly of the Sun */ 117 v_rad = M_rad + C_rad; 118 119 /* Sun radius vector */ 120 R = 1.000001018*(1.0-e*e)/(1.0 + e*cos(v_rad)); /* [a.u] */ 121 (void)R; /* Not used */ 122 123 /* lambda: apparent longitude of the Sun, 124 * relative to the true equinox of the date */ 125 omega_deg = 125.04 - 1934.136*t; 126 omega_rad = MDEG2RAD(omega_deg); 127 lambda_deg = theta_deg - 5.69e-3 - 4.78e-3*sin(omega_rad); 128 lambda_rad = MDEG2RAD(lambda_deg); 129 130 /* epsilon0: mean obliquity of the ecliptic */ 131 d1 = dms2deg(23, 26, 21.448); 132 d2 = dms2deg( 0, 0, 46.8150); 133 d3 = dms2deg( 0, 0, 5.90e-4); 134 d4 = dms2deg( 0, 0, 1.813e-3); 135 epsilon0_deg = d1 - d2*t - d3*t*t + d4*t*t*t; 136 137 /* L: mean longitude of the Sun */ 138 L_deg = 280.4665 + 36000.7698*t; 139 L_deg = fmod(L_deg, 360); 140 L_rad = MDEG2RAD(L_deg); 141 142 /* Lprime: mean longitude of the Moon */ 143 Lprime_deg = 218.3165 + 481267.8813*t; 144 Lprime_deg = fmod(Lprime_deg, 360); 145 Lprime_rad = MDEG2RAD(Lprime_deg); 146 147 /* Delta_epsilon: nutation in obliquity */ 148 d1 = dms2deg(0, 0, 9.20); 149 d2 = dms2deg(0, 0, 5.70e-1); 150 d3 = dms2deg(0, 0, 1.0e-1); 151 d4 = dms2deg(0, 0, 9.0e-2); 152 delta_epsilon_deg = 153 d1*cos(omega_rad) 154 + d2*cos(2.0*L_rad) 155 + d3*cos(2.0*Lprime_rad) 156 - d4*cos(2.0*omega_rad); 157 158 /* epsilon: obliquity of the ecliptic */ 159 epsilon_deg = epsilon0_deg + delta_epsilon_deg + 2.56e-3*cos(omega_rad); 160 epsilon_rad = MDEG2RAD(epsilon_deg); 161 162 /* alpha: right ascension of the Sun */ 163 alpha_rad = atan2(cos(epsilon_rad)*sin(lambda_rad), cos(lambda_rad)); 164 165 /* delta: declination of the Sun */ 166 delta_rad = asin(sin(epsilon_rad)*sin(lambda_rad)); 167 168 /* Setup output data */ 169 celestial->right_ascension = alpha_rad; 170 celestial->declination = delta_rad; 171 } 172 173 /******************************************************************************* 174 * Local function 175 ******************************************************************************/ 176 res_T 177 sun_position_from_earth_meeus 178 (const struct tm* time, /* In UTC */ 179 const struct scem_location* pos, /* Local position */ 180 struct scem_sun_pos* sun) 181 { 182 struct celestial_coords celestial = CELESTIAL_COORDS_NULL; 183 double jd = 0; /* Julian Day */ 184 double theta_deg = 0; 185 double theta_rad = 0; 186 double H_rad = 0; 187 double latitude_rad = 0; 188 189 ASSERT(time && pos && sun); 190 191 jd = julian_day(time); 192 193 sun_pos_celestial_coords(jd, &celestial); 194 195 /* theta: true solar time (taking into account the longitude of the observer) */ 196 theta_deg = 280.1470 + 360.9856235*(jd-2451545.0) + pos->longitude /* [deg] */; 197 theta_deg = fmod(theta_deg, 360); 198 theta_rad = MDEG2RAD(theta_deg); 199 200 /* H: angular distance between the observer and the Sun */ 201 H_rad = theta_rad - celestial.right_ascension; 202 203 /* Zenital angle: angle between the horizontal plane and the Sun direction */ 204 latitude_rad = MDEG2RAD(pos->latitude /* [deg] */); 205 sun->zenith /* [rad] */ = asin 206 ( sin(latitude_rad)*sin(celestial.declination) 207 + cos(latitude_rad)*cos(celestial.declination)*cos(H_rad)); 208 209 /* Azimutal angle: angle between the North and the projection of the Sun 210 * direction on the horizontal plane, with a rotation direction from North to 211 * East. */ 212 sun->azimuth /* [rad] */ = PI + atan2 213 (sin(H_rad), 214 cos(H_rad)*sin(latitude_rad) - tan(celestial.declination)*cos(latitude_rad)); 215 216 return RES_OK; 217 }