From c008dcdff6feb17abe6248c0ffcd03c15defed77 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Thu, 12 Dec 2019 12:35:00 -0500 Subject: [PATCH] Finished C demo program horizon.c. This program demonstrates converting ecliptic coordinates to horizontal coordinates at a given time for a given observer. It searches for the two locations on the horizon where the ecliptic plane intersects it. --- demo/c/.gitignore | 1 + demo/c/README.md | 5 + demo/c/astro_demo_common.c | 4 +- demo/c/build | 12 +- demo/c/horizon.c | 190 ++++++++++++++++++++++++++++++++ demo/c/test/.gitignore | 1 + demo/c/test/horizon_correct.txt | 2 + demo/c/test/test | 4 + 8 files changed, 215 insertions(+), 4 deletions(-) create mode 100644 demo/c/horizon.c create mode 100644 demo/c/test/horizon_correct.txt diff --git a/demo/c/.gitignore b/demo/c/.gitignore index 13072f5e..7a03e711 100644 --- a/demo/c/.gitignore +++ b/demo/c/.gitignore @@ -3,3 +3,4 @@ positions riseset seasons culminate +horizon diff --git a/demo/c/README.md b/demo/c/README.md index 988af27a..b5ed688f 100644 --- a/demo/c/README.md +++ b/demo/c/README.md @@ -22,6 +22,11 @@ Culmination is also the moment a body crosses the *meridian*, the imaginary semi in the sky that passes from due north on the horizon, through the zenith (straight up), and then toward due south on the horizon. +### [Horizon Intersection](horizon.c) +This is a more advanced example. It shows how to use coordinate +transforms and the search mechanism to find where the +ecliptic intersects with an observer's horizon at a given date and time. + --- # [API Reference](../../source/c/) diff --git a/demo/c/astro_demo_common.c b/demo/c/astro_demo_common.c index e43af8cf..c71b6318 100644 --- a/demo/c/astro_demo_common.c +++ b/demo/c/astro_demo_common.c @@ -8,8 +8,8 @@ int ParseArgs(int argc, const char *argv[], astro_observer_t *observer, astro_ti { observer->height = 0.0; - if (1 != sscanf(argv[1], "%lf", &observer->latitude) || - observer->latitude < -90.0 || + if (1 != sscanf(argv[1], "%lf", &observer->latitude) || + observer->latitude < -90.0 || observer->latitude > +90.0) { fprintf(stderr, "ERROR: Invalid latitude '%s' on command line\n", argv[1]); diff --git a/demo/c/build b/demo/c/build index 1def3f88..ec7465d2 100755 --- a/demo/c/build +++ b/demo/c/build @@ -5,8 +5,16 @@ Fail() exit 1 } -for name in moonphase positions riseset seasons culminate; do - gcc -Wall -Werror -o ${name} -I../../source/c ../../source/c/astronomy.c astro_demo_common.c ${name}.c -lm || +if [[ "$1" == "debug" ]]; then + BUILDOPT='-g -O0' +elif [[ -z "$1" ]]; then + BUILDOPT='-O3' +else + Fail "unrecognized command line option" +fi + +for name in moonphase positions riseset seasons culminate horizon; do + gcc ${BUILDOPT} -Wall -Werror -o ${name} -I../../source/c ../../source/c/astronomy.c astro_demo_common.c ${name}.c -lm || Fail "Error building ${name}.c" done diff --git a/demo/c/horizon.c b/demo/c/horizon.c new file mode 100644 index 00000000..e98d273f --- /dev/null +++ b/demo/c/horizon.c @@ -0,0 +1,190 @@ +/* + horizon.c - Don Cross - 2019-12-11 + + Example C program for Astronomy Engine: + https://github.com/cosinekitty/astronomy + + This is a more advanced example. It shows how to use coordinate + transforms and a binary search to find the two azimuths where the + ecliptic intersects with an observer's horizon at a given date and time. +*/ + +#include +#include +#include "astro_demo_common.h" + +static int FindEclipticCrossings(astro_observer_t observer, astro_time_t time); + + +int main(int argc, const char *argv[]) +{ + int error; + astro_observer_t observer; + astro_time_t time; + + if (argc != 4) + { + fprintf(stderr, "USAGE: horizon latitude longitude yyyy-mm-ddThh:mm:ssZ\n"); + return 1; + } + + error = ParseArgs(argc, argv, &observer, &time); + if (error) + return error; + + return FindEclipticCrossings(observer, time); +} + +static int HorizontalCoords( + astro_spherical_t *hor, + double ecliptic_longitude, + astro_time_t time, + astro_rotation_t rot_ecl_hor) +{ + astro_vector_t ecl_vec; + astro_vector_t hor_vec; + astro_spherical_t eclip; + + eclip.status = ASTRO_SUCCESS; + eclip.lat = 0.0; /* being "on the ecliptic plane" means ecliptic latitude is zero. */ + eclip.lon = ecliptic_longitude; + eclip.dist = 1.0; /* any positive distance value will work fine. */ + + /* Convert ecliptic angular coordinates to ecliptic vector. */ + ecl_vec = Astronomy_VectorFromSphere(eclip, time); + if (ecl_vec.status != ASTRO_SUCCESS) + { + fprintf(stderr, "Altitude: error %d converting ecliptic angles to vector.\n", ecl_vec.status); + return 1; + } + + /* Use the rotation matrix to convert ecliptic vector to horizontal vector. */ + hor_vec = Astronomy_RotateVector(rot_ecl_hor, ecl_vec); + if (hor_vec.status != ASTRO_SUCCESS) + { + fprintf(stderr, "Altitude: error %d rotating ecliptic vector to horizontal vector.\n", hor_vec.status); + return 1; + } + + /* Find horizontal angular coordinates, correcting for atmospheric refraction. */ + *hor = Astronomy_HorizonFromVector(hor_vec, REFRACTION_NORMAL); + if (hor->status != ASTRO_SUCCESS) + { + fprintf(stderr, "Altitude: error %d converting horizontal vector to angles.\n", hor->status); + return 1; + } + return 0; /* success */ +} + +static int Search( + double *ecliptic_longitude_crossing, + astro_spherical_t *hor_crossing, + astro_time_t time, + astro_rotation_t rot_ecl_hor, + double e1, double e2) +{ + int error; + double e3; + astro_spherical_t h3; + const double tolerance = 1.0e-6; /* one-millionth of a degree is close enough! */ + + /* + Binary search: find the ecliptic longitude such that the horizontal altitude + ascends through a zero value. The caller must pass e1, e2 such that the altitudes + bound zero in ascending order. + */ + + for(;;) + { + e3 = (e1 + e2) / 2.0; + error = HorizontalCoords(&h3, e3, time, rot_ecl_hor); + if (error) + return error; + + if (fabs(e2-e1) < tolerance) + { + /* We have found the horizon crossing within tolerable limits. */ + *ecliptic_longitude_crossing = e3; + *hor_crossing = h3; + return 0; + } + + if (h3.lat < 0.0) + e1 = e3; + else + e2 = e3; + } +} + +#define NUM_SAMPLES 4 +#define ECLIPLON(i) ((360.0 * (i)) / NUM_SAMPLES) + +static int FindEclipticCrossings(astro_observer_t observer, astro_time_t time) +{ + int i; + astro_rotation_t rot; + astro_spherical_t hor[NUM_SAMPLES]; + + /* + The ecliptic is a celestial circle that describes the mean plane of + the Earth's orbit around the Sun. We use J2000 ecliptic coordinates, + meaning the x-axis is defined to where the plane of the Earth's + equator on January 1, 2000 at noon UTC intersects the ecliptic plane. + The positive x-axis points toward the March equinox. + Calculate a rotation matrix that converts J2000 ecliptic vectors + to horizontal vectors for this observer and time. + */ + rot = Astronomy_Rotation_ECL_HOR(time, observer); + if (rot.status != ASTRO_SUCCESS) + { + fprintf(stderr, "FindEclipticCrossings: error %d trying to find rotation matrix.\n", rot.status); + return 1; + } + + /* + Sample several points around the ecliptic. + Remember the horizontal coordinates for each sample. + */ + for (i=0; i a1) + { + /* Search for the ecliptic longitude and azimuth where altitude ascends through zero. */ + error = Search(&ex, &hx, time, rot, e1, e2); + } + else + { + /* Search for the ecliptic longitude and azimuth where altitude descends through zero. */ + error = Search(&ex, &hx, time, rot, e2, e1); + } + + if (error) + return error; + + if (hx.lon > 0.0 && hx.lon < 180.0) + direction = "ascends "; /* azimuth is more toward the east than the west */ + else + direction = "descends"; /* azimuth is more toward the west than the east */ + + printf("Ecliptic longitude %9.4lf %s through horizon az %9.4lf, alt %12lg\n", ex, direction, hx.lon, hx.lat); + } + } + return 0; +} diff --git a/demo/c/test/.gitignore b/demo/c/test/.gitignore index 91e3269e..d84be8c9 100644 --- a/demo/c/test/.gitignore +++ b/demo/c/test/.gitignore @@ -3,3 +3,4 @@ positions.txt riseset.txt seasons.txt culminate.txt +horizon.txt diff --git a/demo/c/test/horizon_correct.txt b/demo/c/test/horizon_correct.txt new file mode 100644 index 00000000..1023b917 --- /dev/null +++ b/demo/c/test/horizon_correct.txt @@ -0,0 +1,2 @@ +Ecliptic longitude 93.6460 descends through horizon az 296.3883, alt -1.73303e-07 +Ecliptic longitude 274.9704 ascends through horizon az 115.7274, alt 2.10448e-07 diff --git a/demo/c/test/test b/demo/c/test/test index 02a79701..2f52a043 100755 --- a/demo/c/test/test +++ b/demo/c/test/test @@ -28,5 +28,9 @@ echo "Testing example: culminate.c" ./culminate +30 -90 2015-02-28T00:00:00Z > test/culminate.txt || Fail "Error testing culminate.c." diff test/culminate.txt test/culminate_correct.txt || Fail "Error comparing culminate.c output." +echo "Testing example: horizon.c" +./horizon +25.5 -85.3 2016-12-25T12:30:45Z > test/horizon.txt || Fail "Error testing horizon.c." +diff test/horizon.txt test/horizon_correct.txt || Fail "Error comparing horizon.c output." + echo "PASS: C examples" exit 0