scem_psa.c (5562B)
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; /* Decimal hours */ 30 int month; /* Month in [1, 12] */ 31 int year; 32 long aux1; 33 long aux2; 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; 59 double mean_anomaly; 60 double omega; 61 double t; 62 63 ASSERT(jd >= 0 && ecliptic); 64 65 t = (jd - 2451545.0); 66 omega = 2.267127827 - 0.00093003392670 * t; 67 mean_lon_rad = 4.895036035 + 0.01720279602 * t; 68 mean_anomaly = 6.239468336 + 0.01720200135 * t; 69 70 ecliptic->longitude = /* [rad] */ 71 mean_lon_rad 72 + 0.03338320972 * sin(mean_anomaly) 73 + 0.0003497596876 * sin(mean_anomaly * 2) - 0.0001544353226 74 - 0.00000868972936 * sin(omega); 75 76 ecliptic->obliquity = /* [rad] */ 77 0.4090904909 - 6.213605399e-9 * t + 0.00004418094944 * 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; /* [rad] */ 89 double y; 90 double x; 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 = 149597870.7; /* km */ 111 112 double greenwich_mean_side_real_time; 113 double local_mean_side_real_time; 114 double hour_angle; 115 double lat_rad; 116 double cos_lat; 117 double sin_lat; 118 double cos_hour_angle; 119 double azimuth_rad; 120 double zenith_rad; 121 double parallax; 122 double x; 123 double y; 124 double t; 125 double h; 126 127 ASSERT(time && pos && celestial && sun); 128 129 t = (jd - 2451545.0); 130 h = (double)time->tm_hour 131 + ((double)time->tm_min + (double)time->tm_sec / 60.0) / 60.0; 132 133 greenwich_mean_side_real_time = 6.697096103 + 0.06570984737 * t + h; 134 135 local_mean_side_real_time = 136 MDEG2RAD(greenwich_mean_side_real_time * 15 + pos->longitude); 137 138 hour_angle = local_mean_side_real_time - celestial->right_ascension; 139 lat_rad = MDEG2RAD(pos->latitude); 140 cos_lat = cos(lat_rad); 141 sin_lat = sin(lat_rad); 142 cos_hour_angle = cos(hour_angle); 143 144 zenith_rad = acos 145 (cos_lat * cos(celestial->declination) * cos_hour_angle 146 + sin_lat * sin(celestial->declination)); 147 148 y = -sin(hour_angle); 149 x = tan(celestial->declination) * cos_lat - sin_lat*cos_hour_angle; 150 azimuth_rad = atan2(y, x); 151 152 if(azimuth_rad < 0.0) azimuth_rad += 2 * PI; 153 154 parallax = (earth_radius / astro_unit) * sin(zenith_rad); 155 zenith_rad += parallax; 156 157 sun->azimuth = fmod(azimuth_rad, 2 * PI); 158 sun->elevation = PI * 0.5 - fmod(zenith_rad, 2 * PI); 159 } 160 161 /******************************************************************************* 162 * Local symbols 163 ******************************************************************************/ 164 res_T 165 sun_position_from_earth_psa 166 (const struct tm* time, /* In UTC */ 167 const struct scem_location* pos, /* Local position */ 168 struct scem_sun_pos* sun) /* Sun position */ 169 { 170 struct ecliptic_coords ecliptic = ECLIPTIC_COORDS_NULL; 171 struct celestial_coords celestial = CELESTIAL_COORDS_NULL; 172 double jd; /* Julian Day */ 173 174 ASSERT(time && pos && sun); 175 176 jd = julian_day(time); 177 sun_pos_ecliptic_coords(jd, &ecliptic); 178 sun_pos_celestial_coords(&ecliptic, &celestial); 179 sun_pos_local_coords(time, pos, jd, &celestial, sun); 180 181 return RES_OK; 182 }