scem_psa.c (5580B)
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 20 #include <math.h> 21 #include <time.h> 22 23 /******************************************************************************* 24 * Helper functions 25 ******************************************************************************/ 26 static double 27 julian_day(const struct tm* time) /* In UTC */ 28 { 29 double h = 0; /* Decimal hours */ 30 int month = 0; /* Month in [1, 12] */ 31 int year = 0; 32 long aux1 = 0; 33 long aux2 = 0; 34 35 ASSERT(time); 36 37 h = (double)time->tm_hour 38 + ((double)time->tm_min + (double)time->tm_sec/60.0) / 60.0; 39 40 month = time->tm_mon + 1; /* tm_mon start at 0 */ 41 year = time->tm_year + 1900; /* tm_year is relative to year 1900 */ 42 43 aux1 = (month - 14) / 12; 44 aux2 = (1461 * (year + 4800 + aux1)) / 4 45 + (367 * (month - 2 - 12*aux1)) / 12 46 - (3 * ((year + 4900 + aux1)/100)) / 4 47 + time->tm_mday - 32075; 48 49 return (double)aux2 - 0.5 + h/24.0; 50 } 51 52 /* Calculate ecliptic coordinates (ecliptic longitude and obliquity of the 53 * ecliptic in radians but without limiting the angle to be less than 2*Pi 54 * (i.e., the result may be greater than 2*Pi) */ 55 static void 56 sun_pos_ecliptic_coords(const double jd, struct ecliptic_coords* ecliptic) 57 { 58 double mean_lon_rad=0; 59 double mean_anomaly=0; 60 double omega=0; 61 double t=0; 62 63 ASSERT(jd >= 0 && ecliptic); 64 65 t = (jd - 2451545.0); 66 omega = 2.1429 - 0.0010394594*t; 67 mean_lon_rad = 4.8950630 + t*0.017202791698; 68 mean_anomaly = 6.2400600 + t*0.0172019699; 69 70 ecliptic->longitude = /* [rad] */ 71 mean_lon_rad 72 + 0.03341607*sin(mean_anomaly) 73 + 0.00034894*sin(mean_anomaly*2) - 0.0001134 74 - 0.0000203 *sin(omega); 75 76 ecliptic->obliquity = /* [rad] */ 77 0.4090928 - 6.2140e-9*t + 0.0000396*cos(omega); 78 } 79 80 /* Calculate celestial coordinates (right ascension and declination) in radians 81 * * but without limiting the angle to be less than 2*Pi (i.e., the result may 82 * be greater than 2*Pi) */ 83 static void 84 sun_pos_celestial_coords 85 (const struct ecliptic_coords* ecliptic, 86 struct celestial_coords* celestial) 87 { 88 double sin_ecliptic_lon=0; /* [rad] */ 89 double y=0; 90 double x=0; 91 92 sin_ecliptic_lon = sin(ecliptic->longitude); 93 y = cos(ecliptic->obliquity) * sin_ecliptic_lon; 94 x = cos(ecliptic->longitude); 95 96 celestial->right_ascension = atan2(y, x); 97 if(celestial->right_ascension < 0) celestial->right_ascension += 2*PI; 98 celestial->declination = asin(sin(ecliptic->obliquity) * sin_ecliptic_lon); 99 } 100 101 static void 102 sun_pos_local_coords 103 (const struct tm* time, /* In UTC */ 104 const struct scem_location* pos, 105 const double jd, /* Julian day */ 106 const struct celestial_coords* celestial, 107 struct scem_sun_pos* sun) 108 { 109 const double earth_radius = 6371.01; /* [km] */ 110 const double astro_unit = 149597890.0; /* [km] */ 111 112 double greenwitch_mean_side_real_time = 0; 113 double local_mean_side_real_time = 0; 114 double hour_angle = 0; 115 double lat_rad = 0; 116 double cos_lat = 0; 117 double sin_lat = 0; 118 double cos_hour_angle = 0; 119 double azimuth_rad = 0; 120 double zenith_rad = 0; 121 double parallax = 0; 122 double x = 0; 123 double y = 0; 124 double t = 0; 125 126 double h = 0; 127 128 ASSERT(time && pos && celestial && sun); 129 130 t = (jd - 2451545.0); 131 h = (double)time->tm_hour 132 + ((double)time->tm_min + (double)time->tm_sec/60.0) / 60.0; 133 134 greenwitch_mean_side_real_time = 6.6974243242 + 0.0657098283*t + h; 135 136 local_mean_side_real_time = 137 MDEG2RAD(greenwitch_mean_side_real_time*15 + pos->longitude); 138 139 hour_angle = local_mean_side_real_time - celestial->right_ascension; 140 lat_rad = MDEG2RAD(pos->latitude); 141 cos_lat = cos(lat_rad); 142 sin_lat = sin(lat_rad); 143 cos_hour_angle = cos(hour_angle); 144 145 zenith_rad = acos 146 ( cos_lat * cos(celestial->declination) * cos_hour_angle 147 + sin_lat * sin(celestial->declination)); 148 149 y = -sin(hour_angle); 150 x = tan(celestial->declination) * cos_lat - sin_lat*cos_hour_angle; 151 azimuth_rad = atan2(y, x); 152 153 if(azimuth_rad < 0.0) azimuth_rad += 2*PI; 154 155 parallax = (earth_radius / astro_unit) * sin(zenith_rad); 156 zenith_rad += parallax; 157 158 sun->azimuth = fmod(azimuth_rad, 2*PI); 159 sun->zenith = PI*0.5 - fmod(zenith_rad, 2*PI); 160 } 161 162 /******************************************************************************* 163 * Local symbols 164 ******************************************************************************/ 165 res_T 166 sun_position_from_earth_psa 167 (const struct tm* time, /* In UTC */ 168 const struct scem_location* pos, /* Local position */ 169 struct scem_sun_pos* sun) /* Sun position */ 170 { 171 struct ecliptic_coords ecliptic = ECLIPTIC_COORDS_NULL; 172 struct celestial_coords celestial = CELESTIAL_COORDS_NULL; 173 double jd = 0; /* Julian Day */ 174 175 ASSERT(time && pos && sun); 176 177 jd = julian_day(time); 178 sun_pos_ecliptic_coords(jd, &ecliptic); 179 sun_pos_celestial_coords(&ecliptic, &celestial); 180 sun_pos_local_coords(time, pos, jd, &celestial, sun); 181 182 return RES_OK; 183 }