From 528dad54cf85828e8cbcbcd94a573e34a2da9499 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 20 May 2019 14:02:29 -0400 Subject: [PATCH] Implemented C function SunPosition. Not yet tested. --- generate/template/astronomy.c | 101 ++++++++++++++++++++++++++++++++++ source/c/astronomy.c | 101 ++++++++++++++++++++++++++++++++++ source/c/astronomy.h | 12 ++++ 3 files changed, 214 insertions(+) diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index bee7071c..3ee6c641 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -52,6 +52,8 @@ static const double AU = 1.4959787069098932e+11; /* astronomical unit in static const double KM_PER_AU = 1.4959787069098932e+8; static const double ANGVEL = 7.2921150e-5; +static astro_ecliptic_t RotateEquatorialToEcliptic(const double pos[3], double obliq_radians); + double Astronomy_VectorLength(astro_vector_t vector) { return sqrt(vector.x*vector.x + vector.y*vector.y + vector.z*vector.z); @@ -99,6 +101,14 @@ static astro_equatorial_t EquError(astro_status_t status) return equ; } +static astro_ecliptic_t EclError(astro_status_t status) +{ + astro_ecliptic_t ecl; + ecl.status = status; + ecl.ex = ecl.ey = ecl.ez = ecl.elat = ecl.elon = NAN; + return ecl; +} + typedef struct { double mjd; @@ -1529,6 +1539,97 @@ astro_horizon_t Astronomy_Horizon( } + +/* + X = not yet implemented + T = still needs testing + + ------------------------------------------- + + X AngleFromSun + X Ecliptic + X EclipticLongitude + X Elongation + X GetPeformanceMetrics + X Illumination + X LongitudeFromSun + X MoonPhase + X NextLunarApsis + X NextMoonQuarter + X ResetPerformanceMetrics + X Search + X SearchHourAngle + X SearchLunarApsis + X SearchMaxElongation + X SearchMoonPhase + X SearchMoonQuarter + X SearchPeakMagnitude + X SearchRelativeLongitude + X SearchRiseSet + X SearchSunLongitude + X Seasons + T SunPosition +*/ + +astro_ecliptic_t Astronomy_SunPosition(astro_time_t observation_time) +{ + astro_time_t time; + astro_vector_t earth2000; + double sun2000[3]; + double stemp[3]; + double sun_ofdate[3]; + double true_obliq; + + /* Correct for light travel time from the Sun. */ + /* Otherwise season calculations (equinox, solstice) will all be early by about 8 minutes! */ + time = Astronomy_AddDays(observation_time, -1.0 / C_AUDAY); + + earth2000 = CalcEarth(time); + if (earth2000.status != ASTRO_SUCCESS) + return EclError(earth2000.status); + + /* Convert heliocentric location of Earth to geocentric location of Sun. */ + sun2000[0] = -earth2000.x; + sun2000[1] = -earth2000.y; + sun2000[2] = -earth2000.z; + + /* Convert to equatorial Cartesian coordinates of date. */ + precession(0.0, sun2000, time.tt, stemp); + nutation(time, 0, stemp, sun_ofdate); + + /* Convert equatorial coordinates to ecliptic coordinates. */ + true_obliq = DEG2RAD * e_tilt(time).tobl; + return RotateEquatorialToEcliptic(sun_ofdate, true_obliq); +} + +static astro_ecliptic_t RotateEquatorialToEcliptic(const double pos[3], double obliq_radians) +{ + astro_ecliptic_t ecl; + double cos_ob, sin_ob; + double xyproj; + + cos_ob = cos(obliq_radians); + sin_ob = sin(obliq_radians); + + ecl.ex = +pos[0]; + ecl.ey = +pos[1]*cos_ob + pos[2]*sin_ob; + ecl.ez = -pos[1]*sin_ob + pos[2]*cos_ob; + + xyproj = sqrt(ecl.ex*ecl.ex + ecl.ey*ecl.ey); + if (xyproj > 0.0) + { + ecl.elon = RAD2DEG * atan2(ecl.ey, ecl.ex); + if (ecl.elon < 0.0) + ecl.elon += 360.0; + } + else + ecl.elon = 0.0; + + ecl.elat = RAD2DEG * atan2(ecl.ez, xyproj); + ecl.status = ASTRO_SUCCESS; + return ecl; +} + #ifdef __cplusplus } #endif diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 340bc535..05c79dee 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -52,6 +52,8 @@ static const double AU = 1.4959787069098932e+11; /* astronomical unit in static const double KM_PER_AU = 1.4959787069098932e+8; static const double ANGVEL = 7.2921150e-5; +static astro_ecliptic_t RotateEquatorialToEcliptic(const double pos[3], double obliq_radians); + double Astronomy_VectorLength(astro_vector_t vector) { return sqrt(vector.x*vector.x + vector.y*vector.y + vector.z*vector.z); @@ -99,6 +101,14 @@ static astro_equatorial_t EquError(astro_status_t status) return equ; } +static astro_ecliptic_t EclError(astro_status_t status) +{ + astro_ecliptic_t ecl; + ecl.status = status; + ecl.ex = ecl.ey = ecl.ez = ecl.elat = ecl.elon = NAN; + return ecl; +} + typedef struct { double mjd; @@ -2585,6 +2595,97 @@ astro_horizon_t Astronomy_Horizon( } + +/* + X = not yet implemented + T = still needs testing + + ------------------------------------------- + + X AngleFromSun + X Ecliptic + X EclipticLongitude + X Elongation + X GetPeformanceMetrics + X Illumination + X LongitudeFromSun + X MoonPhase + X NextLunarApsis + X NextMoonQuarter + X ResetPerformanceMetrics + X Search + X SearchHourAngle + X SearchLunarApsis + X SearchMaxElongation + X SearchMoonPhase + X SearchMoonQuarter + X SearchPeakMagnitude + X SearchRelativeLongitude + X SearchRiseSet + X SearchSunLongitude + X Seasons + T SunPosition +*/ + +astro_ecliptic_t Astronomy_SunPosition(astro_time_t observation_time) +{ + astro_time_t time; + astro_vector_t earth2000; + double sun2000[3]; + double stemp[3]; + double sun_ofdate[3]; + double true_obliq; + + /* Correct for light travel time from the Sun. */ + /* Otherwise season calculations (equinox, solstice) will all be early by about 8 minutes! */ + time = Astronomy_AddDays(observation_time, -1.0 / C_AUDAY); + + earth2000 = CalcEarth(time); + if (earth2000.status != ASTRO_SUCCESS) + return EclError(earth2000.status); + + /* Convert heliocentric location of Earth to geocentric location of Sun. */ + sun2000[0] = -earth2000.x; + sun2000[1] = -earth2000.y; + sun2000[2] = -earth2000.z; + + /* Convert to equatorial Cartesian coordinates of date. */ + precession(0.0, sun2000, time.tt, stemp); + nutation(time, 0, stemp, sun_ofdate); + + /* Convert equatorial coordinates to ecliptic coordinates. */ + true_obliq = DEG2RAD * e_tilt(time).tobl; + return RotateEquatorialToEcliptic(sun_ofdate, true_obliq); +} + +static astro_ecliptic_t RotateEquatorialToEcliptic(const double pos[3], double obliq_radians) +{ + astro_ecliptic_t ecl; + double cos_ob, sin_ob; + double xyproj; + + cos_ob = cos(obliq_radians); + sin_ob = sin(obliq_radians); + + ecl.ex = +pos[0]; + ecl.ey = +pos[1]*cos_ob + pos[2]*sin_ob; + ecl.ez = -pos[1]*sin_ob + pos[2]*cos_ob; + + xyproj = sqrt(ecl.ex*ecl.ex + ecl.ey*ecl.ey); + if (xyproj > 0.0) + { + ecl.elon = RAD2DEG * atan2(ecl.ey, ecl.ex); + if (ecl.elon < 0.0) + ecl.elon += 360.0; + } + else + ecl.elon = 0.0; + + ecl.elat = RAD2DEG * atan2(ecl.ez, xyproj); + ecl.status = ASTRO_SUCCESS; + return ecl; +} + #ifdef __cplusplus } #endif diff --git a/source/c/astronomy.h b/source/c/astronomy.h index 41d669f5..3dbc21b0 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -97,6 +97,17 @@ typedef struct } astro_equatorial_t; +typedef struct +{ + astro_status_t status; + double ex; + double ey; + double ez; + double elat; + double elon; +} +astro_ecliptic_t; + typedef struct { double azimuth; @@ -132,6 +143,7 @@ astro_equatorial_t Astronomy_Equator( int ofdate, int aberration ); +astro_ecliptic_t Astronomy_SunPosition(astro_time_t time); astro_horizon_t Astronomy_Horizon(astro_time_t time, astro_observer_t observer, double ra, double dec, astro_refraction_t refraction); #ifdef __cplusplus