/* * Partial C port of Greg Miller's public domain astro library (gmiller@gregmiller.net) 2019 * https://github.com/gmiller123456/astrogreg * * Ported by Joey Castillo for Sensor Watch * https://github.com/joeycastillo/Sensor-Watch/ * * Public Domain * * THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * SOFTWARE. */ #include #include #include #include #include "astrolib.h" #include "vsop87a_milli.h" double astro_convert_utc_to_tt(double jd) ; double astro_get_GMST(double ut1); astro_cartesian_coordinates_t astro_subtract_cartesian(astro_cartesian_coordinates_t a, astro_cartesian_coordinates_t b); astro_cartesian_coordinates_t astro_rotate_from_vsop_to_J2000(astro_cartesian_coordinates_t c); astro_matrix_t astro_get_x_rotation_matrix(double r); astro_matrix_t astro_get_y_rotation_matrix(double r); astro_matrix_t astro_get_z_rotation_matrix(double r); astro_matrix_t astro_transpose_matrix(astro_matrix_t m); astro_matrix_t astro_dot_product(astro_matrix_t a, astro_matrix_t b); astro_matrix_t astro_get_precession_matrix(double jd); astro_cartesian_coordinates_t astro_matrix_multiply(astro_cartesian_coordinates_t v, astro_matrix_t m); astro_cartesian_coordinates_t astro_convert_geodedic_latlon_to_ITRF_XYZ(double lat, double lon, double height); astro_cartesian_coordinates_t astro_convert_ITRF_to_GCRS(astro_cartesian_coordinates_t r, double ut1); astro_cartesian_coordinates_t astro_convert_coordinates_from_meters_to_AU(astro_cartesian_coordinates_t c); astro_cartesian_coordinates_t astro_get_observer_geocentric_coords(double jd, double lat, double lon); astro_cartesian_coordinates_t astro_get_body_coordinates(astro_body_t bodyNum, double et); astro_cartesian_coordinates_t astro_get_body_coordinates_light_time_adjusted(astro_body_t body, astro_cartesian_coordinates_t origin, double t); astro_equatorial_coordinates_t astro_convert_cartesian_to_polar(astro_cartesian_coordinates_t xyz); //Special "Math.floor()" function used by convertDateToJulianDate() static double _astro_special_floor(double d) { if(d > 0) { return floor(d); } return floor(d) - 1; } double astro_convert_date_to_julian_date(uint16_t year, uint8_t month, uint8_t day, uint8_t hour, uint8_t minute, uint8_t second) { if (month < 3){ year = year - 1; month = month + 12; } double b = 0; if (!(year < 1582 || (year == 1582 && (month < 10 || (month == 10 && day < 5))))) { double a = _astro_special_floor(year / 100.0); b = 2 - a + _astro_special_floor(a / 4.0); } double jd = _astro_special_floor(365.25 * (year + 4716)) + _astro_special_floor(30.6001 * (month + 1)) + day + b - 1524.5; jd += hour / 24.0; jd += minute / 24.0 / 60.0; jd += second / 24.0 / 60.0 / 60.0; return jd; } //Return all values in radians. //The positions are adjusted for the parallax of the Earth, and the offset of the observer from the Earth's center //All input and output angles are in radians! astro_equatorial_coordinates_t astro_get_ra_dec(double jd, astro_body_t body, double lat, double lon, bool calculate_precession) { double jdTT = astro_convert_utc_to_tt(jd); double t = astro_convert_jd_to_julian_millenia_since_j2000(jdTT); // Get current position of Earth and the target body astro_cartesian_coordinates_t earth_coords = astro_get_body_coordinates(ASTRO_BODY_EARTH, t); astro_cartesian_coordinates_t body_coords = astro_get_body_coordinates_light_time_adjusted(body, earth_coords, t); // Convert to Geocentric coordinate body_coords = astro_subtract_cartesian(body_coords, earth_coords); //Rotate ecliptic coordinates to J2000 coordinates body_coords = astro_rotate_from_vsop_to_J2000(body_coords); astro_matrix_t precession; // TODO: rotate body for precession, nutation and bias if(calculate_precession) { precession = astro_get_precession_matrix(jdTT); body_coords = astro_matrix_multiply(body_coords, precession); } //Convert to topocentric astro_cartesian_coordinates_t observerXYZ = astro_get_observer_geocentric_coords(jdTT, lat, lon); if(calculate_precession) { //TODO: rotate observerXYZ for precession, nutation and bias astro_matrix_t precessionInv = astro_transpose_matrix(precession); observerXYZ = astro_matrix_multiply(observerXYZ, precessionInv); } body_coords = astro_subtract_cartesian(body_coords, observerXYZ); //Convert to topocentric RA DEC by converting from cartesian coordinates to polar coordinates astro_equatorial_coordinates_t retval = astro_convert_cartesian_to_polar(body_coords); retval.declination = M_PI/2.0 - retval.declination; //Dec. Offset to make 0 the equator, and the poles +/-90 deg if(retval.right_ascension < 0) retval.right_ascension += 2*M_PI; //Ensure RA is positive return retval; } //Converts a Julian Date in UTC to Terrestrial Time (TT) double astro_convert_utc_to_tt(double jd) { //Leap seconds are hard coded, should be updated from the IERS website for other times //TAI = UTC + leap seconds (e.g. 32) //TT=TAI + 32.184 //return jd + (32.0 + 32.184) / 24.0 / 60.0 / 60.0; return jd + (37.0 + 32.184) / 24.0 / 60.0 / 60.0; /* https://data.iana.org/time-zones/tzdb-2018a/leap-seconds.list 2272060800 10 # 1 Jan 1972 2287785600 11 # 1 Jul 1972 2303683200 12 # 1 Jan 1973 2335219200 13 # 1 Jan 1974 2366755200 14 # 1 Jan 1975 2398291200 15 # 1 Jan 1976 2429913600 16 # 1 Jan 1977 2461449600 17 # 1 Jan 1978 2492985600 18 # 1 Jan 1979 2524521600 19 # 1 Jan 1980 2571782400 20 # 1 Jul 1981 2603318400 21 # 1 Jul 1982 2634854400 22 # 1 Jul 1983 2698012800 23 # 1 Jul 1985 2776982400 24 # 1 Jan 1988 2840140800 25 # 1 Jan 1990 2871676800 26 # 1 Jan 1991 2918937600 27 # 1 Jul 1992 2950473600 28 # 1 Jul 1993 2982009600 29 # 1 Jul 1994 3029443200 30 # 1 Jan 1996 3076704000 31 # 1 Jul 1997 3124137600 32 # 1 Jan 1999 3345062400 33 # 1 Jan 2006 3439756800 34 # 1 Jan 2009 3550089600 35 # 1 Jul 2012 3644697600 36 # 1 Jul 2015 3692217600 37 # 1 Jan 2017 */ } double astro_convert_jd_to_julian_millenia_since_j2000(double jd) { return (jd - 2451545.0) / 365250.0; } astro_cartesian_coordinates_t astro_subtract_cartesian(astro_cartesian_coordinates_t a, astro_cartesian_coordinates_t b) { astro_cartesian_coordinates_t retval; retval.x = a.x - b.x; retval.y = a.y - b.y; retval.z = a.z - b.z; return retval; } // Performs the rotation from ecliptic coordinates to J2000 coordinates for the given vector x astro_cartesian_coordinates_t astro_rotate_from_vsop_to_J2000(astro_cartesian_coordinates_t c) { /* From VSOP87.doc X +1.000000000000 +0.000000440360 -0.000000190919 X Y = -0.000000479966 +0.917482137087 -0.397776982902 Y Z FK5 0.000000000000 +0.397776982902 +0.917482137087 Z VSOP87A */ astro_cartesian_coordinates_t t; t.x = c.x + c.y * 0.000000440360 + c.z * -0.000000190919; t.y = c.x * -0.000000479966 + c.y * 0.917482137087 + c.z * -0.397776982902; t.z = c.y * 0.397776982902 + c.z * 0.917482137087; return t; } double astro_get_GMST(double ut1) { double D = ut1 - 2451545.0; double T = D/36525.0; double gmst = fmod((280.46061837 + 360.98564736629*D + 0.000387933*T*T - T*T*T/38710000.0), 360.0); if(gmst<0) { gmst+=360; } return gmst/15; } static astro_matrix_t _astro_get_empty_matrix() { astro_matrix_t t; for(uint8_t i = 0; i < 3 ; i++) { for(uint8_t j = 0 ; j < 3 ; j++) { t.elements[i][j] = 0; } } return t; } //Gets a rotation matrix about the x axis. Angle R is in radians astro_matrix_t astro_get_x_rotation_matrix(double r) { astro_matrix_t t = _astro_get_empty_matrix(); t.elements[0][0]=1; t.elements[0][1]=0; t.elements[0][2]=0; t.elements[1][0]=0; t.elements[1][1]=cos(r); t.elements[1][2]=sin(r); t.elements[2][0]=0; t.elements[2][1]=-sin(r); t.elements[2][2]=cos(r); return t; } //Gets a rotation matrix about the y axis. Angle R is in radians astro_matrix_t astro_get_y_rotation_matrix(double r) { astro_matrix_t t = _astro_get_empty_matrix(); t.elements[0][0]=cos(r); t.elements[0][1]=0; t.elements[0][2]=-sin(r); t.elements[1][0]=0; t.elements[1][1]=1; t.elements[1][2]=0; t.elements[2][0]=sin(r); t.elements[2][1]=0; t.elements[2][2]=cos(r); return t; } //Gets a rotation matrix about the z axis. Angle R is in radians astro_matrix_t astro_get_z_rotation_matrix(double r) { astro_matrix_t t = _astro_get_empty_matrix(); t.elements[0][0]=cos(r); t.elements[0][1]=sin(r); t.elements[0][2]=0; t.elements[1][0]=-sin(r); t.elements[1][1]=cos(r); t.elements[1][2]=0; t.elements[2][0]=0; t.elements[2][1]=0; t.elements[2][2]=1; return t; } void astro_print_matrix(char * title, astro_matrix_t matrix); void astro_print_matrix(char * title, astro_matrix_t matrix) { printf("%s\n", title); for(uint8_t i = 0; i < 3 ; i++) { printf("\t"); for(uint8_t j = 0 ; j < 3 ; j++) { printf("%12f", matrix.elements[i][j]); } printf("\n"); } printf("\n"); } astro_matrix_t astro_dot_product(astro_matrix_t a, astro_matrix_t b) { astro_matrix_t retval; for(uint8_t i = 0; i < 3 ; i++) { for(uint8_t j = 0 ; j < 3 ; j++) { double temp = 0; for(uint8_t k = 0; k < 3 ; k++) { temp += a.elements[i][k] * b.elements[k][j]; } retval.elements[i][j]=temp; } } return retval; } astro_matrix_t astro_transpose_matrix(astro_matrix_t m) { astro_matrix_t retval; for(uint8_t i = 0; i < 3 ; i++) { for(uint8_t j = 0 ; j < 3 ; j++) { retval.elements[i][j] = m.elements[j][i]; } } return retval; } astro_matrix_t astro_get_precession_matrix(double jd) { //2006 IAU Precession. Implemented from IERS Technical Note No 36 ch5. //https://www.iers.org/SharedDocs/Publikationen/EN/IERS/Publications/tn/TechnNote36/tn36_043.pdf?__blob=publicationFile&v=1 double t = (jd - 2451545.0) / 36525.0; //5.2 const double Arcsec2Radians = M_PI/180.0/60.0/60.0; //Converts arc seconds used in equations below to radians double e0 = 84381.406 * Arcsec2Radians; //5.6.4 double omegaA = e0 + ((-0.025754 + (0.0512623 + (-0.00772503 + (-0.000000467 + 0.0000003337*t) * t) * t) * t) * t) * Arcsec2Radians; //5.39 double psiA = ((5038.481507 + (-1.0790069 + (-0.00114045 + (0.000132851 - 0.0000000951*t) * t) * t) * t) * t) * Arcsec2Radians; //5.39 double chiA = ((10.556403 + (-2.3814292 + (-0.00121197 + (0.000170663 - 0.0000000560*t) * t) * t) * t) * t) * Arcsec2Radians; //5.40 //Rotation matrix from 5.4.5 //(R1(−e0) · R3(psiA) · R1(omegaA) · R3(−chiA)) //Above eq rotates from "of date" to J2000, so we reverse the signs to go from J2000 to "of date" astro_matrix_t m1 = astro_get_x_rotation_matrix(e0); astro_matrix_t m2 = astro_get_z_rotation_matrix(-psiA); astro_matrix_t m3 = astro_get_x_rotation_matrix(-omegaA); astro_matrix_t m4 = astro_get_z_rotation_matrix(chiA); astro_matrix_t m5 = astro_dot_product(m4, m3); astro_matrix_t m6 = astro_dot_product(m5, m2); astro_matrix_t precessionMatrix = astro_dot_product(m6, m1); return precessionMatrix; } astro_cartesian_coordinates_t astro_matrix_multiply(astro_cartesian_coordinates_t v, astro_matrix_t m) { astro_cartesian_coordinates_t t; t.x = v.x*m.elements[0][0] + v.y*m.elements[0][1] + v.z*m.elements[0][2]; t.y = v.x*m.elements[1][0] + v.y*m.elements[1][1] + v.z*m.elements[1][2]; t.z = v.x*m.elements[2][0] + v.y*m.elements[2][1] + v.z*m.elements[2][2]; return t; } //Converts cartesian XYZ coordinates to polar (e.g. J2000 xyz to Right Accention and Declication) astro_equatorial_coordinates_t astro_convert_cartesian_to_polar(astro_cartesian_coordinates_t xyz) { astro_equatorial_coordinates_t t; t.distance = sqrt(xyz.x * xyz.x + xyz.y * xyz.y + xyz.z * xyz.z); t.declination = acos(xyz.z / t.distance); t.right_ascension = atan2(xyz.y, xyz.x); if(t.declination < 0) t.declination += 2 * M_PI; if(t.right_ascension < 0) t.right_ascension += 2 * M_PI; return t; } //Convert Geodedic Lat Lon to geocentric XYZ position vector //All angles are input as radians astro_cartesian_coordinates_t astro_convert_geodedic_latlon_to_ITRF_XYZ(double lat, double lon, double height) { //Algorithm from Explanatory Supplement to the Astronomical Almanac 3rd ed. P294 const double a = 6378136.6; const double f = 1 / 298.25642; const double C = sqrt(((cos(lat)*cos(lat)) + (1.0-f)*(1.0-f) * (sin(lat)*sin(lat)))); const double S = (1-f)*(1-f)*C; double h = height; astro_cartesian_coordinates_t r; r.x = (a*C+h) * cos(lat) * cos(lon); r.y = (a*C+h) * cos(lat) * sin(lon); r.z = (a*S+h) * sin(lat); return r; } //Convert position vector to celestial "of date" system. //g(t)=R3(-GAST) r //(Remember to use UT1 for GAST, not ET) //All angles are input and output as radians astro_cartesian_coordinates_t astro_convert_ITRF_to_GCRS(astro_cartesian_coordinates_t r, double ut1) { //This is a simple rotation matrix implemenation about the Z axis, rotation angle is -GMST double GMST = astro_get_GMST(ut1); GMST =- GMST * 15.0 * M_PI / 180.0; astro_matrix_t m = astro_get_z_rotation_matrix(GMST); astro_cartesian_coordinates_t t = astro_matrix_multiply(r, m); return t; } astro_cartesian_coordinates_t astro_convert_coordinates_from_meters_to_AU(astro_cartesian_coordinates_t c) { astro_cartesian_coordinates_t t; t.x = c.x / 1.49597870691E+11; t.y = c.y / 1.49597870691E+11; t.z = c.z / 1.49597870691E+11; return t; } astro_cartesian_coordinates_t astro_get_observer_geocentric_coords(double jd, double lat, double lon) { astro_cartesian_coordinates_t r = astro_convert_geodedic_latlon_to_ITRF_XYZ(lat, lon,0); r = astro_convert_ITRF_to_GCRS(r, jd); r = astro_convert_coordinates_from_meters_to_AU(r); return r; } //Returns a body's cartesian coordinates centered on the Sun. //Requires vsop87a_milli_js, if you wish to use a different version of VSOP87, replace the class name vsop87a_milli below astro_cartesian_coordinates_t astro_get_body_coordinates(astro_body_t body, double et) { astro_cartesian_coordinates_t retval = {0}; double coords[3]; switch(body) { case ASTRO_BODY_SUN: return retval; //Sun is at the center for vsop87a case ASTRO_BODY_MERCURY: vsop87a_milli_getMercury(et, coords); break; case ASTRO_BODY_VENUS: vsop87a_milli_getVenus(et, coords); break; case ASTRO_BODY_EARTH: vsop87a_milli_getEarth(et, coords); break; case ASTRO_BODY_MARS: vsop87a_milli_getMars(et, coords); break; case ASTRO_BODY_JUPITER: vsop87a_milli_getJupiter(et, coords); break; case ASTRO_BODY_SATURN: vsop87a_milli_getSaturn(et, coords); break; case ASTRO_BODY_URANUS: vsop87a_milli_getUranus(et, coords); break; case ASTRO_BODY_NEPTUNE: vsop87a_milli_getNeptune(et, coords); break; case ASTRO_BODY_EMB: vsop87a_milli_getEmb(et, coords); break; case ASTRO_BODY_MOON: { double earth_coords[3]; double emb_coords[3]; vsop87a_milli_getEarth(et, earth_coords); vsop87a_milli_getEmb(et, emb_coords); vsop87a_milli_getMoon(earth_coords, emb_coords, coords); } break; } retval.x = coords[0]; retval.y = coords[1]; retval.z = coords[2]; return retval; } astro_cartesian_coordinates_t astro_get_body_coordinates_light_time_adjusted(astro_body_t body, astro_cartesian_coordinates_t origin, double t) { //Get current position of body astro_cartesian_coordinates_t body_coords = astro_get_body_coordinates(body, t); double newT = t; for(uint8_t i = 0 ; i < 2 ; i++) { //Calculate light time to body body_coords = astro_subtract_cartesian(body_coords, origin); double distance = sqrt(body_coords.x*body_coords.x + body_coords.y*body_coords.y + body_coords.z*body_coords.z); distance *= 1.496e+11; //Convert from AU to meters double lightTime = distance / 299792458.0; //Convert light time to Julian Millenia, and subtract it from the original value of t newT -= lightTime / 24.0 / 60.0 / 60.0 / 365250.0; //Recalculate body position adjusted for light time body_coords = astro_get_body_coordinates(body, newT); } return body_coords; } astro_horizontal_coordinates_t astro_ra_dec_to_alt_az(double jd, double lat, double lon, double ra, double dec) { double GMST = astro_get_GMST(jd) * M_PI/180.0 * 15.0; double h = GMST + lon - ra; double sina = sin(dec)*sin(lat) + cos(dec)*cos(h)*cos(lat); double a = asin(sina); double cosAz = (sin(dec)*cos(lat) - cos(dec)*cos(h)*sin(lat)) / cos(a); double Az = acos(cosAz); if(sin(h) > 0) Az = 2.0*M_PI - Az; astro_horizontal_coordinates_t retval; retval.altitude = a; retval.azimuth = Az; return retval; } double astro_degrees_to_radians(double degrees) { return degrees * M_PI / 180; } double astro_radians_to_degrees(double radians) { return radians * 180.0 / M_PI; } astro_angle_dms_t astro_radians_to_dms(double radians) { astro_angle_dms_t retval; int8_t sign = (radians < 0) ? -1 : 1; double degrees = fabs(astro_radians_to_degrees(radians)); retval.degrees = (uint16_t)degrees; double temp = 60.0 * (degrees - retval.degrees); retval.minutes = (uint8_t)temp; retval.seconds = (uint8_t)round(60.0 * (temp - retval.minutes)); if (retval.seconds > 59) { retval.seconds = 0.0; retval.minutes++; } if (retval.minutes > 59) { retval.minutes = 0; retval.degrees++; } degrees *= sign; return retval; } astro_angle_hms_t astro_radians_to_hms(double radians) { astro_angle_hms_t retval; double degrees = astro_radians_to_degrees(radians); double temp = degrees / 15.0; retval.hours = (uint8_t)temp; temp = 60.0 * (temp - retval.hours); retval.minutes = (uint8_t)temp; retval.seconds = (uint8_t)round(60.0 * (temp - retval.minutes)); if (retval.seconds > 59) { retval.seconds = 0; retval.minutes++; } if (retval.minutes > 59) { retval.minutes = 0; retval.hours++; } return retval; }