mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-05-24 08:45:55 -04:00
Code Issues Packages Projects Releases Wiki Activity

Added Astronomy_HorizonFromVector, Astronomy_Rotation_EQD_HOR.

Browse Source 
Astronomy_HorizonFromVector is a specialized variant of
Astronomy_SphereFromVector that flips the orientation of the
azimuth angle to the more traditional clockwise-from-north
direction used in navigation and cartography.
It also allows the same optional refraction correction as
Astronomy_Horizon.

Astronomy_Rotation_EQD_HOR converts a equatorial-of-date vector
to a horizontal vector. The horizontal vector has the following
components:
    x = North
    y = West
    z = Zenith

Removed trailing whitespace in generate.c.
This commit is contained in:
Don Cross
2019-12-09 11:06:02 -05:00
parent 212c60b633
commit 4399e72334
6 changed files with 482 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files

60
generate/ctest.c
View File

@@ -2229,6 +2229,59 @@ fail:
return error;
}
static int Test_EQD_HOR(astro_body_t body)
{
astro_time_t time;
astro_observer_t observer;
astro_equatorial_t eqd;
astro_horizon_t hor;
astro_rotation_t rot;
astro_vector_t vec_eqd, vec_hor;
astro_spherical_t sphere;
double diff_az, diff_alt;
int error;
/* Use existing functions to calculate horizontal coordinates of the body for the time+observer. */
time = Astronomy_MakeTime(1970, 12, 13, 5, 15, 0.0);
observer = Astronomy_MakeObserver(-37.0, +45.0, 0.0);
CHECK_EQU(eqd, Astronomy_Equator(body, &time, observer, EQUATOR_OF_DATE, ABERRATION));
printf("Test_EQD_HOR %s: OFDATE ra=%0.3lf, dec=%0.3lf\n", Astronomy_BodyName(body), eqd.ra, eqd.dec);
hor = Astronomy_Horizon(&time, observer, eqd.ra, eqd.dec, REFRACTION_NORMAL);
/* Calculate the position of the body as an equatorial vector of date. */
sphere.status = ASTRO_SUCCESS;
sphere.dist = eqd.dist;
sphere.lat = eqd.dec;
sphere.lon = 15.0 * eqd.ra;
CHECK_VECTOR(vec_eqd, Astronomy_VectorFromSphere(sphere, time));
/* Calculate rotation matrix to convert equatorial J2000 vector to horizontal vector. */
rot = Astronomy_Rotation_EQD_HOR(time, observer);
/* Rotate the equator of date vector to a horizontal vector. */
CHECK_VECTOR(vec_hor, Astronomy_RotateVector(rot, vec_eqd));
/* Convert the horizontal vector to horizontal angular coordinates. */
sphere = Astronomy_HorizonFromVector(vec_hor, REFRACTION_NORMAL);
CHECK_STATUS(sphere);
diff_alt = fabs(sphere.lat - hor.altitude);
diff_az = fabs(sphere.lon - hor.azimuth);
printf("Test_EQD_HOR %s: trusted alt=%0.3lf, az=%0.3lf; test alt=%0.3lf, az=%0.3lf; diff_alt=%lg, diff_az=%lg\n",
Astronomy_BodyName(body), hor.altitude, hor.azimuth, sphere.lat, sphere.lon, diff_alt, diff_az);
if (diff_alt > 2.0e-14)
{
fprintf(stderr, "Test_EQD_HOR: EXCESSIVE HORIZONTAL ERROR.\n");
return 1;
}
error = 0;
fail:
return error;
}
static int RotationTest(void)
{
int error;
@@ -2262,12 +2315,19 @@ static int RotationTest(void)
CHECK(TestSpin(0.0, 0.0, -45.0, +1, 0, 0, +0.7071067811865476, -0.7071067811865476, 0));
CHECK(Test_EQJ_ECL());
CHECK(Test_EQJ_EQD(BODY_MERCURY));
CHECK(Test_EQJ_EQD(BODY_VENUS));
CHECK(Test_EQJ_EQD(BODY_MARS));
CHECK(Test_EQJ_EQD(BODY_JUPITER));
CHECK(Test_EQJ_EQD(BODY_SATURN));
CHECK(Test_EQD_HOR(BODY_MERCURY));
CHECK(Test_EQD_HOR(BODY_VENUS));
CHECK(Test_EQD_HOR(BODY_MARS));
CHECK(Test_EQD_HOR(BODY_JUPITER));
CHECK(Test_EQD_HOR(BODY_SATURN));
error = 0;
fail:
return error;

74
generate/generate.c
View File

@@ -97,14 +97,14 @@ static int UnitTestChebyshev(void);
#define NEPTUNE_PERIHELION 29.81
#define PLUTO_PERIHELION 29.658
/*
ErrorRadius[body] is the worst-case (smallest) distance
in AU over which a position error affects Earth-based observations.
/*
ErrorRadius[body] is the worst-case (smallest) distance
in AU over which a position error affects Earth-based observations.
For example, an error in the Earth/Moon Barycenter (EMB) affects
observations of Venus the most, because that planet comes closest
to Earth at 0.277 AU.
*/
const double ErrorRadius[] =
const double ErrorRadius[] =
{
EARTH_PERIHELION - MERCURY_APHELION, /* 0 = Mercury */
EARTH_PERIHELION - VENUS_APHELION, /* 1 = Venus */
@@ -138,7 +138,7 @@ int main(int argc, const char *argv[])
static int PrintUsage(void)
{
fprintf(stderr,
fprintf(stderr,
"\n"
"USAGE:\n"
"\n"
@@ -188,8 +188,8 @@ static int NovasBodyPos(double jd, int body, double pos[3])
{
double factor;
/*
The caller is asking for the Earth's position or the Moon's position.
/*
The caller is asking for the Earth's position or the Moon's position.
NOVAS does not directly represent either body.
Instead, we have to calculate the Earth or Moon using
the Earth/Moon Barycenter (EMB) and the Geocentric Moon (GM).
@@ -301,7 +301,7 @@ static int SearchVsop(int body)
winner_terms = trunc_terms;
winner_arcmin = max_arcmin;
winner_threshold = threshold;
error = SaveVsopFile(&model);
if (error) goto fail;
}
@@ -326,7 +326,7 @@ fail:
static int SaveVsopFile(const vsop_model_t *model)
{
char filename[100];
snprintf(filename, sizeof(filename), "output/vsop_%d.txt", (int)model->body);
return VsopWriteTrunc(model, filename);
}
@@ -436,13 +436,13 @@ static int PositionArcminError(int body, double jd, const double a[3], const dou
return 0;
}
if (body == BODY_MOON)
if (body == BODY_MOON)
{
/*
Measure the position of the moon as seen from the *SURFACE* of the Earth,
at the point along a line from the geocenter to the object.
This is closer than the geocenter or the barycenter, so the anglular error is larger.
For other bodies, there is no significant difference between
For other bodies, there is no significant difference between
topocentric error and geocentric error.
*/
const double scale = (ERAD / AU) / VectorLength(a);
@@ -465,7 +465,7 @@ static int PositionArcminError(int body, double jd, const double a[3], const dou
adiff.v[0] = a[0];
adiff.v[1] = a[1];
adiff.v[2] = a[2];
bdiff.t = jd; /* doesn't matter, but I don't like leaving undefined memory */
bdiff.v[0] = b[0];
bdiff.v[1] = b[1];
@@ -477,7 +477,7 @@ static int PositionArcminError(int body, double jd, const double a[3], const dou
if (body == BODY_EMB || body == BODY_EARTH)
{
/*
/*
For the Earth/Moon Barycenter (EMB=2), or Earth=11, we use a special estimate
of angular error. An error in the Earth's position affects the observed
position of all other bodies in the solar system.
@@ -528,7 +528,7 @@ calc:
adiff.v[0] = a[0] - epos[0];
adiff.v[1] = a[1] - epos[1];
adiff.v[2] = a[2] - epos[2];
bdiff.t = jd; /* doesn't matter, but I don't like leaving undefined memory */
bdiff.v[0] = b[0] - epos[0];
bdiff.v[1] = b[1] - epos[1];
@@ -578,7 +578,7 @@ static int MeasureError(const char *inFileName, int nsamples, error_stats_t *sta
memset(stats, 0, sizeof(error_stats_t));
error = EphFileOpen(&reader, inFileName);
if (error)
if (error)
{
fprintf(stderr, "MeasureError: Error %d trying to open file: %s\n", error, inFileName);
goto fail;
@@ -614,7 +614,7 @@ static int MeasureError(const char *inFileName, int nsamples, error_stats_t *sta
/* Calculate "correct" position at the time 'jd' (according to NOVAS). */
error = NovasBodyPos(jd, reader.body, npos);
if (error) goto fail;
/* Calculate our approximation of the position at the time 'jd'. */
x = ChebScale(record.jdStart, record.jdStart + record.jdDelta, jd);
ChebApprox(record.numpoly, 3, record.coeff, x, apos);
@@ -682,7 +682,7 @@ static int Resample(int body, const char *outFileName, int npoly, int startYear,
}
error = ChebInit(&encoder, npoly, VECTOR_DIM);
if (error)
if (error)
{
fprintf(stderr, "ERROR %d returned by ChebInit()\n", error);
goto fail;
@@ -745,7 +745,7 @@ static int SampleFunc(const void *context, double jd, double pos[CHEB_MAX_DIM])
double jed[2];
double sun_pos[3];
double vel[3]; /* we don't care about velocities... ignored */
jed[0] = jd;
jed[1] = 0.0;
error = state(jed, (short)c->body, pos, vel);
@@ -839,9 +839,9 @@ static int CheckEcliptic(const char *filename, int lnum, const char *line, doubl
double earth_ecl[3];
double geo[3];
double tt, jd, lon, dist, blon, elon, diff;
*arcmin = 99999.0;
*outbody = VSOP_INVALID_BODY;
*outbody = VSOP_INVALID_BODY;
/* Example line: "e Jupiter opp <tt> <au>" */
if (4 != sscanf(line, "e %9[A-Za-z] %9[a-z] %lf %lf", name, event, &tt, &dist))
@@ -883,7 +883,7 @@ static int CheckEcliptic(const char *filename, int lnum, const char *line, doubl
{
fprintf(stderr, "CheckEcliptic: Invalid distance %0.3lf AU for inferior conjunction; line %d file %s\n", dist, lnum, filename);
return 1;
}
}
}
else if (!strcmp(event, "sup"))
{
@@ -895,16 +895,16 @@ static int CheckEcliptic(const char *filename, int lnum, const char *line, doubl
{
fprintf(stderr, "CheckEcliptic: Invalid distance %0.3lf AU for superior conjunction; line %d file %s\n", dist, lnum, filename);
return 1;
}
}
}
else
{
fprintf(stderr, "CheckEcliptic: Invalid event code '%s' on line %d of file %s\n", event, lnum, filename);
return 1;
}
}
blon = EclipticLongitude(body_ecl);
elon = EclipticLongitude(earth_ecl);
elon = EclipticLongitude(earth_ecl);
diff = LongitudeOffset((elon - blon) - lon);
*arcmin = fabs(diff * 60.0);
@@ -943,14 +943,14 @@ static int CheckTestVector(const char *filename, int lnum, const char *line, dou
}
error = NovasBodyPos(jd, body, npos);
if (error)
if (error)
{
fprintf(stderr, "CheckTestVector: NovasBodyPos returned %d on line %d of file %s\n", error, lnum, filename);
return error;
}
error = PositionArcminError(body, jd, npos, xpos, arcmin);
if (error)
if (error)
{
fprintf(stderr, "CheckTestVector: PositionArcminError returned %d on line %d of file %s\n", error, lnum, filename);
return error;
@@ -996,7 +996,7 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
return 1;
}
if (body == BODY_EARTH || body == BODY_EMB)
if (body == BODY_EARTH || body == BODY_EMB)
{
fprintf(stderr, "CheckSkyPos: Cannot calculate sky position of body '%s'\n", name);
return 1;
@@ -1007,7 +1007,7 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
else if (body == BODY_SUN)
bodyIndex = 10;
else
bodyIndex = 1 + body;
bodyIndex = 1 + body;
error = make_object(0, (short)bodyIndex, name, NULL, &obj);
if (error)
@@ -1031,12 +1031,12 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
delta_ra = fabs(sky_ast.ra - ra);
if (delta_ra > 12.0)
{
/*
/*
Sometimes the two RA values can straddle the 24-hour mark.
For example, one of them is 0.001 and the other 23.999.
The actual error is then 0.002 hours, not 23.998.
In general, it is never "fair" to call the error greater than
12 hours or 180 degrees.
12 hours or 180 degrees.
*/
delta_ra = 24.0 - delta_ra;
}
@@ -1065,10 +1065,10 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
return error;
}
equ2hor(jd_utc, jd_tt-jd_utc, 1, 0.0, 0.0,
&location->on_surf, sky_dat.ra, sky_dat.dec, 0,
equ2hor(jd_utc, jd_tt-jd_utc, 1, 0.0, 0.0,
&location->on_surf, sky_dat.ra, sky_dat.dec, 0,
&check_zenith, &check_azimuth, &check_ra, &check_dec);
check_altitude = 90.0 - check_zenith;
delta_az = fabs(azimuth - check_azimuth);
@@ -1077,10 +1077,10 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
delta_az *= 60.0;
/* Just like RA, diminish the azimuth error as altitude approaches zenith/nadir... */
delta_az *= cos(check_altitude * DEG2RAD);
delta_alt = 60.0 * (altitude - check_altitude);
*arcmin_hor = sqrt(delta_az*delta_az + delta_alt*delta_alt);
*arcmin_hor = sqrt(delta_az*delta_az + delta_alt*delta_alt);
if (*arcmin_hor > 0.9)
{
@@ -1211,7 +1211,7 @@ static int CheckTestOutput(const char *filename)
CHECK(ParseObserver(filename, lnum, line, &location));
break;
case 'v': /* heliocentric cartesian vector represented in J2000 equatorial plane */
case 'v': /* heliocentric cartesian vector represented in J2000 equatorial plane */
CHECK(CheckTestVector(filename, lnum, line, &arcmin_helio, &body));
UpdateErrorStats(&bundle[body].helio, arcmin_helio);
break;
@@ -1226,7 +1226,7 @@ static int CheckTestOutput(const char *filename)
CHECK(CheckEcliptic(filename, lnum, line, &arcmin_eclip, &body));
UpdateErrorStats(&bundle[body].eclip, arcmin_eclip);
break;
default:
fprintf(stderr, "CheckTestOutput: Invalid first character on line %d of file %s\n", lnum, filename);
error = 1;

182
generate/template/astronomy.c
View File

@@ -4191,7 +4191,7 @@ astro_vector_t Astronomy_VectorFromSphere(astro_spherical_t sphere, astro_time_t
* Given a Cartesian vector, returns latitude, longitude, and distance.
*
* @param vector
* Cartesian vector to be converted to spherical coordinates
* Cartesian vector to be converted to spherical coordinates.
*
* @return
* Spherical coordinates that are equivalent to the given vector.
@@ -4235,6 +4235,100 @@ astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector)
return sphere;
}
/**
* @brief Converts Cartesian coordinates to horizontal coordinates.
*
* Given a horizontal Cartesian vector, returns horizontal azimuth and altitude.
*
* *IMPORTANT:* This function differs from #Astronomy_SphereFromVector in two ways:
* - `Astronomy_SphereFromVector` returns a `lon` value that represents azimuth defined counterclockwise
* from north (e.g., west = +90), but this function represents a clockwise rotation
* (e.g., east = +90). The difference is because `Astronomy_SphereFromVector` is intended
* to preserve the vector "right-hand rule", while this function defines azimuth in a more
* traditional way as used in navigation and cartography.
* - This function optionally corrects for atmospheric refraction, while `Astronomy_SphereFromVector`
* does not.
*
* The returned structure will contain the azimuth in `lon`.
* It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees.
*
* The altitude will be stored in `lat`.
*
* The distance to the observed object is stored in `dist`,
* and is expressed in astronomical units (AU).
*
* @param vector
* Cartesian vector to be converted to horizontal coordinates.
*
* @param refraction
* `REFRACTION_NORMAL`: correct altitude for atmospheric refraction (recommended).
* `REFRACTION_NONE`: no atmospheric refraction correction is performed.
* `REFRACTION_JPLHOR`: for JPL Horizons compatibility testing only; not recommended for normal use.
*
* @return
* If successful, `status` hold `ASTRO_SUCCESS` and the other fields are valid as described above.
* Otherwise `status` holds an error code and the other fields are undefined.
*/
astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refraction_t refraction)
{
astro_spherical_t sphere;
sphere = Astronomy_SphereFromVector(vector);
if (sphere.status == ASTRO_SUCCESS)
{
/* Convert azimuth from counterclockwise-from-north to clockwise-from-north. */
sphere.lon = 360.0 - sphere.lon;
if (sphere.lon >= 360.0)
sphere.lon -= 360.0;
else if (sphere.lon < 0.0)
sphere.lon += 360.0;
if (refraction == REFRACTION_NONE)
{
/* No atmospheric correction needed; do nothing. */
}
else if (refraction == REFRACTION_NORMAL || refraction == REFRACTION_JPLHOR)
{
double zd, refr, hd;
zd = 90.0 - sphere.lat;
// http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf
// JPL Horizons says it uses refraction algorithm from
// Meeus "Astronomical Algorithms", 1991, p. 101-102.
// I found the following Go implementation:
// https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go
// This is a translation from the function "Saemundsson" there.
// I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon.
// This is important because the 'refr' formula below goes crazy near hd = -5.11.
hd = sphere.lat;
if (hd < -1.0)
hd = -1.0;
refr = (1.02 / tan((hd+10.3/(hd+5.11))*DEG2RAD)) / 60.0;
if (refraction == REFRACTION_NORMAL && zd > 91.0)
{
// In "normal" mode we gradually reduce refraction toward the nadir
// so that we never get an altitude angle less than -90 degrees.
// When horizon angle is -1 degrees, zd = 91, and the factor is exactly 1.
// As zd approaches 180 (the nadir), the fraction approaches 0 linearly.
refr *= (180.0 - zd) / 89.0;
}
sphere.lat += refr;
}
else
{
/* Invalid refraction option. */
return SphereError(ASTRO_INVALID_PARAMETER);
}
}
return sphere;
}
/**
* @brief
* Applies a rotation to a vector, yielding a rotated vector.
@@ -4268,28 +4362,6 @@ astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t
}
#if 0
/**
* @brief
* Calculates a rotation matrix from equatorial of-date system to equatorial J2000.
*
* This is one of the family of MakeRotation functions that convert
* from one orientation system to another.
* Source: EQD = equatorial system, using equator of date.
* Target: EQJ = equatorial system, using equator at J2000 epoch.
*
* @param time
* The date and time defining the equator of the source orientation.
*
* @return
* A rotation matrix that converts EQD to EQJ.
*/
astro_rotation_t Astronomy_MakeRotation_EQD_EQJ(astro_time_t time)
{
}
#endif
/**
* @brief
* Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL).
@@ -4390,6 +4462,70 @@ astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time)
return Astronomy_CombineRotation(nut, prec);
}
/**
* @brief
* Calculates a rotation matrix from equatorial of-date (EQD) to horizontal (HOR).
*
* This is one of the family of functions that returns a rotation matrix
* for converting from one orientation to another.
* Source: EQD = equatorial system, using equator of the specified date/time.
* Target: HOR = horizontal system.
*
* Use #Astronomy_HorizonFromVector to convert the return value
* to a traditional altitude/azimuth pair.
*
* @param time
* The date and time at which the Earth's equator applies.
*
* @param observer
* A location near the Earth's mean sea level that defines the observer's horizon.
*
* @return
* A rotation matrix that converts EQD to HOR at `time` and for `observer`.
* The components of the horizontal vector are:
* x = north, y = west, z = zenith (straight up from the observer).
* These components are chosen so that the "right-hand rule" works for the vector
* and so that north represents the direction where azimuth = 0.
*/
astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer)
{
astro_rotation_t rot;
double uze[3], une[3], uwe[3];
double uz[3], un[3], uw[3];
double spin_angle;
double sinlat = sin(observer.latitude * DEG2RAD);
double coslat = cos(observer.latitude * DEG2RAD);
double sinlon = sin(observer.longitude * DEG2RAD);
double coslon = cos(observer.longitude * DEG2RAD);
uze[0] = coslat * coslon;
uze[1] = coslat * sinlon;
uze[2] = sinlat;
une[0] = -sinlat * coslon;
une[1] = -sinlat * sinlon;
une[2] = coslat;
uwe[0] = sinlon;
uwe[1] = -coslon;
uwe[2] = 0.0;
spin_angle = -15.0 * sidereal_time(&time);
spin(spin_angle, uze, uz);
spin(spin_angle, une, un);
spin(spin_angle, uwe, uw);
rot.rot[0][0] = un[0]; rot.rot[1][0] = un[1]; rot.rot[2][0] = un[2];
rot.rot[0][1] = uw[0]; rot.rot[1][1] = uw[1]; rot.rot[2][1] = uw[2];
rot.rot[0][2] = uz[0]; rot.rot[1][2] = uz[1]; rot.rot[2][2] = uz[2];
rot.status = ASTRO_SUCCESS;
return rot;
}
#ifdef __cplusplus
}
#endif

66
source/c/README.md
View File

@@ -473,6 +473,43 @@ This function optionally corrects for atmospheric refraction. For most uses, it
---
<a name="Astronomy_HorizonFromVector"></a>
### Astronomy_HorizonFromVector(vector, refraction) &#8658; [`astro_spherical_t`](#astro_spherical_t)
**Converts Cartesian coordinates to horizontal coordinates.**
Given a horizontal Cartesian vector, returns horizontal azimuth and altitude.
*IMPORTANT:* This function differs from [`Astronomy_SphereFromVector`](#Astronomy_SphereFromVector) in two ways:
- `Astronomy_SphereFromVector` returns a `lon` value that represents azimuth defined counterclockwise from north (e.g., west = +90), but this function represents a clockwise rotation (e.g., east = +90). The difference is because `Astronomy_SphereFromVector` is intended to preserve the vector "right-hand rule", while this function defines azimuth in a more traditional way as used in navigation and cartography.
- This function optionally corrects for atmospheric refraction, while `Astronomy_SphereFromVector` does not.
The returned structure will contain the azimuth in `lon`. It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees.
The altitude will be stored in `lat`.
The distance to the observed object is stored in `dist`, and is expressed in astronomical units (AU).
**Returns:** If successful, `status` hold `ASTRO_SUCCESS` and the other fields are valid as described above. Otherwise `status` holds an error code and the other fields are undefined.
| Type | Parameter | Description |
| --- | --- | --- |
| [`astro_vector_t`](#astro_vector_t) | `vector` | Cartesian vector to be converted to horizontal coordinates. |
| [`astro_refraction_t`](#astro_refraction_t) | `refraction` | `REFRACTION_NORMAL`: correct altitude for atmospheric refraction (recommended). `REFRACTION_NONE`: no atmospheric refraction correction is performed. `REFRACTION_JPLHOR`: for JPL Horizons compatibility testing only; not recommended for normal use. |
---
<a name="Astronomy_Illumination"></a>
@@ -764,6 +801,33 @@ This is one of the family of functions that returns a rotation matrix for conver
---
<a name="Astronomy_Rotation_EQD_HOR"></a>
### Astronomy_Rotation_EQD_HOR(time, observer) &#8658; [`astro_rotation_t`](#astro_rotation_t)
**Calculates a rotation matrix from equatorial of-date (EQD) to horizontal (HOR).**
This is one of the family of functions that returns a rotation matrix for converting from one orientation to another. Source: EQD = equatorial system, using equator of the specified date/time. Target: HOR = horizontal system.
Use [`Astronomy_HorizonFromVector`](#Astronomy_HorizonFromVector) to convert the return value to a traditional altitude/azimuth pair.
**Returns:** A rotation matrix that converts EQD to HOR at `time` and for `observer`. The components of the horizontal vector are: x = north, y = west, z = zenith (straight up from the observer). These components are chosen so that the "right-hand rule" works for the vector and so that north represents the direction where azimuth = 0.
| Type | Parameter | Description |
| --- | --- | --- |
| [`astro_time_t`](#astro_time_t) | `time` | The date and time at which the Earth's equator applies. |
| [`astro_observer_t`](#astro_observer_t) | `observer` | A location near the Earth's mean sea level that defines the observer's horizon. |
---
<a name="Astronomy_Rotation_EQJ_ECL"></a>
@@ -1162,7 +1226,7 @@ Given a Cartesian vector, returns latitude, longitude, and distance.
| Type | Parameter | Description |
| --- | --- | --- |
| [`astro_vector_t`](#astro_vector_t) | `vector` | Cartesian vector to be converted to spherical coordinates |
| [`astro_vector_t`](#astro_vector_t) | `vector` | Cartesian vector to be converted to spherical coordinates. |

182
source/c/astronomy.c
View File

@@ -5430,7 +5430,7 @@ astro_vector_t Astronomy_VectorFromSphere(astro_spherical_t sphere, astro_time_t
* Given a Cartesian vector, returns latitude, longitude, and distance.
*
* @param vector
* Cartesian vector to be converted to spherical coordinates
* Cartesian vector to be converted to spherical coordinates.
*
* @return
* Spherical coordinates that are equivalent to the given vector.
@@ -5474,6 +5474,100 @@ astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector)
return sphere;
}
/**
* @brief Converts Cartesian coordinates to horizontal coordinates.
*
* Given a horizontal Cartesian vector, returns horizontal azimuth and altitude.
*
* *IMPORTANT:* This function differs from #Astronomy_SphereFromVector in two ways:
* - `Astronomy_SphereFromVector` returns a `lon` value that represents azimuth defined counterclockwise
* from north (e.g., west = +90), but this function represents a clockwise rotation
* (e.g., east = +90). The difference is because `Astronomy_SphereFromVector` is intended
* to preserve the vector "right-hand rule", while this function defines azimuth in a more
* traditional way as used in navigation and cartography.
* - This function optionally corrects for atmospheric refraction, while `Astronomy_SphereFromVector`
* does not.
*
* The returned structure will contain the azimuth in `lon`.
* It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees.
*
* The altitude will be stored in `lat`.
*
* The distance to the observed object is stored in `dist`,
* and is expressed in astronomical units (AU).
*
* @param vector
* Cartesian vector to be converted to horizontal coordinates.
*
* @param refraction
* `REFRACTION_NORMAL`: correct altitude for atmospheric refraction (recommended).
* `REFRACTION_NONE`: no atmospheric refraction correction is performed.
* `REFRACTION_JPLHOR`: for JPL Horizons compatibility testing only; not recommended for normal use.
*
* @return
* If successful, `status` hold `ASTRO_SUCCESS` and the other fields are valid as described above.
* Otherwise `status` holds an error code and the other fields are undefined.
*/
astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refraction_t refraction)
{
astro_spherical_t sphere;
sphere = Astronomy_SphereFromVector(vector);
if (sphere.status == ASTRO_SUCCESS)
{
/* Convert azimuth from counterclockwise-from-north to clockwise-from-north. */
sphere.lon = 360.0 - sphere.lon;
if (sphere.lon >= 360.0)
sphere.lon -= 360.0;
else if (sphere.lon < 0.0)
sphere.lon += 360.0;
if (refraction == REFRACTION_NONE)
{
/* No atmospheric correction needed; do nothing. */
}
else if (refraction == REFRACTION_NORMAL || refraction == REFRACTION_JPLHOR)
{
double zd, refr, hd;
zd = 90.0 - sphere.lat;
// http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf
// JPL Horizons says it uses refraction algorithm from
// Meeus "Astronomical Algorithms", 1991, p. 101-102.
// I found the following Go implementation:
// https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go
// This is a translation from the function "Saemundsson" there.
// I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon.
// This is important because the 'refr' formula below goes crazy near hd = -5.11.
hd = sphere.lat;
if (hd < -1.0)
hd = -1.0;
refr = (1.02 / tan((hd+10.3/(hd+5.11))*DEG2RAD)) / 60.0;
if (refraction == REFRACTION_NORMAL && zd > 91.0)
{
// In "normal" mode we gradually reduce refraction toward the nadir
// so that we never get an altitude angle less than -90 degrees.
// When horizon angle is -1 degrees, zd = 91, and the factor is exactly 1.
// As zd approaches 180 (the nadir), the fraction approaches 0 linearly.
refr *= (180.0 - zd) / 89.0;
}
sphere.lat += refr;
}
else
{
/* Invalid refraction option. */
return SphereError(ASTRO_INVALID_PARAMETER);
}
}
return sphere;
}
/**
* @brief
* Applies a rotation to a vector, yielding a rotated vector.
@@ -5507,28 +5601,6 @@ astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t
}
#if 0
/**
* @brief
* Calculates a rotation matrix from equatorial of-date system to equatorial J2000.
*
* This is one of the family of MakeRotation functions that convert
* from one orientation system to another.
* Source: EQD = equatorial system, using equator of date.
* Target: EQJ = equatorial system, using equator at J2000 epoch.
*
* @param time
* The date and time defining the equator of the source orientation.
*
* @return
* A rotation matrix that converts EQD to EQJ.
*/
astro_rotation_t Astronomy_MakeRotation_EQD_EQJ(astro_time_t time)
{
}
#endif
/**
* @brief
* Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL).
@@ -5629,6 +5701,70 @@ astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time)
return Astronomy_CombineRotation(nut, prec);
}
/**
* @brief
* Calculates a rotation matrix from equatorial of-date (EQD) to horizontal (HOR).
*
* This is one of the family of functions that returns a rotation matrix
* for converting from one orientation to another.
* Source: EQD = equatorial system, using equator of the specified date/time.
* Target: HOR = horizontal system.
*
* Use #Astronomy_HorizonFromVector to convert the return value
* to a traditional altitude/azimuth pair.
*
* @param time
* The date and time at which the Earth's equator applies.
*
* @param observer
* A location near the Earth's mean sea level that defines the observer's horizon.
*
* @return
* A rotation matrix that converts EQD to HOR at `time` and for `observer`.
* The components of the horizontal vector are:
* x = north, y = west, z = zenith (straight up from the observer).
* These components are chosen so that the "right-hand rule" works for the vector
* and so that north represents the direction where azimuth = 0.
*/
astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer)
{
astro_rotation_t rot;
double uze[3], une[3], uwe[3];
double uz[3], un[3], uw[3];
double spin_angle;
double sinlat = sin(observer.latitude * DEG2RAD);
double coslat = cos(observer.latitude * DEG2RAD);
double sinlon = sin(observer.longitude * DEG2RAD);
double coslon = cos(observer.longitude * DEG2RAD);
uze[0] = coslat * coslon;
uze[1] = coslat * sinlon;
uze[2] = sinlat;
une[0] = -sinlat * coslon;
une[1] = -sinlat * sinlon;
une[2] = coslat;
uwe[0] = sinlon;
uwe[1] = -coslon;
uwe[2] = 0.0;
spin_angle = -15.0 * sidereal_time(&time);
spin(spin_angle, uze, uz);
spin(spin_angle, une, un);
spin(spin_angle, uwe, uw);
rot.rot[0][0] = un[0]; rot.rot[1][0] = un[1]; rot.rot[2][0] = un[2];
rot.rot[0][1] = uw[0]; rot.rot[1][1] = uw[1]; rot.rot[2][1] = uw[2];
rot.rot[0][2] = uz[0]; rot.rot[1][2] = uz[1]; rot.rot[2][2] = uz[2];
rot.status = ASTRO_SUCCESS;
return rot;
}
#ifdef __cplusplus
}
#endif

3
source/c/astronomy.h
View File

@@ -616,11 +616,12 @@ astro_rotation_t Astronomy_InverseRotation(astro_rotation_t rotation);
astro_rotation_t Astronomy_CombineRotation(astro_rotation_t a, astro_rotation_t b);
astro_vector_t Astronomy_VectorFromSphere(astro_spherical_t sphere, astro_time_t time);
astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector);
astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refraction_t refraction);
astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t vector);
astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time);
//astro_rotation_t Astronomy_Rotation_EQD_ECL(astro_time_t time);
//astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer);
astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer);
astro_rotation_t Astronomy_Rotation_EQJ_EQD(astro_time_t time);
astro_rotation_t Astronomy_Rotation_EQJ_ECL(void);
//astro_rotation_t Astronomy_Rotation_EQJ_HOR(astro_time_t time, astro_observer_t observer);