mirror of
https://github.com/cosinekitty/astronomy.git
synced 2026-05-24 08:45:55 -04:00
Implemented C version of Horizon function. Fixed bug in C GeoPos.
The C version of GeoPos was returning an ante-dated time value, not the time the caller asked about. In other words, it was returning the time when the observed body emitted the light the observer sees, not the time the observation was made on the Earth, like it was supposed to. Fortunately, my unit test caught that.
This commit is contained in:
@@ -6,9 +6,12 @@
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#include "astronomy.h"
|
||||
|
||||
#define PI 3.14159265358979323846
|
||||
|
||||
#define CHECK(x) do{if(0 != (error = (x))) goto fail;}while(0)
|
||||
|
||||
static int CheckVector(int lnum, astro_vector_t v)
|
||||
@@ -23,15 +26,23 @@ static int CheckVector(int lnum, astro_vector_t v)
|
||||
|
||||
#define CHECK_VECTOR(v,x) CHECK(CheckVector(__LINE__, ((v) = (x))))
|
||||
|
||||
static int Test_Geo(void);
|
||||
static int Test_AstroTime(void);
|
||||
static int AstroCheck(void);
|
||||
|
||||
int main()
|
||||
int main(int argc, const char *argv[])
|
||||
{
|
||||
int error;
|
||||
int error = 0;
|
||||
|
||||
CHECK(Test_AstroTime());
|
||||
CHECK(AstroCheck());
|
||||
if (argc == 2 && !strcmp(argv[1], "geo"))
|
||||
{
|
||||
CHECK(Test_Geo()); /* ad hoc test for debugging */
|
||||
}
|
||||
else
|
||||
{
|
||||
CHECK(Test_AstroTime());
|
||||
CHECK(AstroCheck());
|
||||
}
|
||||
|
||||
fail:
|
||||
printf("ctest exiting with %d\n", error);
|
||||
@@ -75,6 +86,7 @@ static int AstroCheck(void)
|
||||
astro_body_t body;
|
||||
astro_vector_t pos;
|
||||
astro_sky_t sky;
|
||||
astro_horizon_t hor;
|
||||
astro_observer_t observer = Astronomy_MakeObserver(28.0, -81.0, 10.0);
|
||||
|
||||
outfile = fopen(filename, "wt");
|
||||
@@ -101,9 +113,10 @@ static int AstroCheck(void)
|
||||
if (body != BODY_EARTH)
|
||||
{
|
||||
CHECK_VECTOR(pos, Astronomy_GeoVector(body, time));
|
||||
|
||||
/* FIXFIXFIX: add SkyPos, Horizon calls here; output as 's' record. */
|
||||
sky = Astronomy_SkyPos(pos, observer);
|
||||
hor = Astronomy_Horizon(sky.t, observer, sky.ofdate.ra, sky.ofdate.dec, REFRACTION_NONE);
|
||||
fprintf(outfile, "s %s %0.16lf %0.16lf %0.16lf %0.16lf %0.16lf %0.16lf %0.16lf\n",
|
||||
Astronomy_BodyName(body), pos.t.tt, pos.t.ut, sky.j2000.ra, sky.j2000.dec, sky.j2000.dist, hor.azimuth, hor.altitude);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -111,10 +124,12 @@ static int AstroCheck(void)
|
||||
CHECK_VECTOR(pos, Astronomy_GeoMoon(time));
|
||||
fprintf(outfile, "v GM %0.16lf %0.16lf %0.16lf %0.16lf\n", pos.t.tt, pos.x, pos.y, pos.z);
|
||||
|
||||
/* FIXFIXFIX: Test SkyPos, Horizon here; output GM 's' record. */
|
||||
sky = Astronomy_SkyPos(pos, observer);
|
||||
hor = Astronomy_Horizon(sky.t, observer, sky.ofdate.ra, sky.ofdate.dec, REFRACTION_NONE);
|
||||
fprintf(outfile, "s GM %0.16lf %0.16lf %0.16lf %0.16lf %0.16lf %0.16lf %0.16lf\n",
|
||||
pos.t.tt, pos.t.ut, sky.j2000.ra, sky.j2000.dec, sky.j2000.dist, hor.azimuth, hor.altitude);
|
||||
|
||||
time = Astronomy_AddDays(time, 10.03141592);
|
||||
time = Astronomy_AddDays(time, 10.0 + PI/100.0);
|
||||
}
|
||||
|
||||
(void)sky;
|
||||
@@ -124,3 +139,18 @@ fail:
|
||||
fclose(outfile);
|
||||
return error;
|
||||
}
|
||||
|
||||
static int Test_Geo(void)
|
||||
{
|
||||
int error = 0;
|
||||
astro_time_t time;
|
||||
astro_vector_t pos;
|
||||
|
||||
time.ut = -109492.2564886306936387;
|
||||
time.tt = -109492.2562455624283757;
|
||||
|
||||
CHECK_VECTOR(pos, Astronomy_GeoVector(BODY_MERCURY, time));
|
||||
printf("#(C ) GeoVector tt=%0.16lf ut=%0.16lf x=%0.16lf y=%0.16lf z=%0.16lf\n", time.tt, time.ut, pos.x, pos.y, pos.z);
|
||||
fail:
|
||||
return error;
|
||||
}
|
||||
|
||||
@@ -1015,7 +1015,6 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
|
||||
/* First call to place: ask for astrometric coordinates (J2000) : coord_sys=3 */
|
||||
error = place(jd_tt, &obj, location, jd_tt-jd_utc, 3, 1, &sky_ast);
|
||||
if (error)
|
||||
|
||||
{
|
||||
fprintf(stderr, "CheckSkyPos: place(3) returned %d on line %d of file %s\n", error, lnum, filename);
|
||||
return error;
|
||||
|
||||
@@ -639,6 +639,32 @@ static void geo_pos(astro_time_t time, astro_observer_t observer, double outpos[
|
||||
precession(time.tt, pos2, 0.0, outpos);
|
||||
}
|
||||
|
||||
static void spin(double angle, const double pos1[3], double vec2[3])
|
||||
{
|
||||
double angr = angle * DEG2RAD;
|
||||
double cosang = cos(angr);
|
||||
double sinang = sin(angr);
|
||||
double xx = cosang;
|
||||
double yx = sinang;
|
||||
double zx = 0;
|
||||
double xy = -sinang;
|
||||
double yy = cosang;
|
||||
double zy = 0;
|
||||
double xz = 0;
|
||||
double yz = 0;
|
||||
double zz = 1;
|
||||
|
||||
vec2[0] = xx*pos1[0] + yx*pos1[1] + zx*pos1[2];
|
||||
vec2[1] = xy*pos1[0] + yy*pos1[1] + zy*pos1[2];
|
||||
vec2[2] = xz*pos1[0] + yz*pos1[1] + zz*pos1[2];
|
||||
}
|
||||
|
||||
static void ter2cel(astro_time_t time, const double vec1[3], double vec2[3])
|
||||
{
|
||||
double gast = sidereal_time(time);
|
||||
spin(-15.0 * gast, vec1, vec2);
|
||||
}
|
||||
|
||||
/*------------------ CalcMoon ------------------*/
|
||||
|
||||
#define DECLARE_PASCAL_ARRAY_1(elemtype,name,xmin,xmax) \
|
||||
@@ -1191,8 +1217,6 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time)
|
||||
double dt;
|
||||
int iter;
|
||||
|
||||
vector.t = time;
|
||||
|
||||
switch (body)
|
||||
{
|
||||
case BODY_EARTH:
|
||||
@@ -1228,7 +1252,7 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time)
|
||||
ltime2 = Astronomy_AddDays(time, -Astronomy_VectorLength(vector) / C_AUDAY);
|
||||
dt = fabs(ltime2.tt - ltime.tt);
|
||||
if (dt < 1.0e-9)
|
||||
return vector;
|
||||
goto finished; /* Tricky: ensures we patch 'vector.t' with current time, not ante-dated time. */
|
||||
|
||||
ltime = ltime2;
|
||||
}
|
||||
@@ -1236,6 +1260,8 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time)
|
||||
break;
|
||||
}
|
||||
|
||||
finished:
|
||||
vector.t = time;
|
||||
return vector;
|
||||
}
|
||||
|
||||
@@ -1262,6 +1288,126 @@ astro_sky_t Astronomy_SkyPos(astro_vector_t gc_vector, astro_observer_t observer
|
||||
return sky;
|
||||
}
|
||||
|
||||
astro_horizon_t Astronomy_Horizon(
|
||||
astro_time_t time, astro_observer_t observer, double ra, double dec, astro_refraction_t refraction)
|
||||
{
|
||||
astro_horizon_t hor;
|
||||
double uze[3], une[3], uwe[3];
|
||||
double uz[3], un[3], uw[3];
|
||||
double p[3], pz, pn, pw, proj;
|
||||
double az, zd;
|
||||
|
||||
double sinlat = sin(observer.latitude * DEG2RAD);
|
||||
double coslat = cos(observer.latitude * DEG2RAD);
|
||||
double sinlon = sin(observer.longitude * DEG2RAD);
|
||||
double coslon = cos(observer.longitude * DEG2RAD);
|
||||
double sindc = sin(dec * DEG2RAD);
|
||||
double cosdc = cos(dec * DEG2RAD);
|
||||
double sinra = sin(ra * 15 * DEG2RAD);
|
||||
double cosra = cos(ra * 15 * 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;
|
||||
|
||||
ter2cel(time, uze, uz);
|
||||
ter2cel(time, une, un);
|
||||
ter2cel(time, uwe, uw);
|
||||
|
||||
p[0] = cosdc * cosra;
|
||||
p[1] = cosdc * sinra;
|
||||
p[2] = sindc;
|
||||
|
||||
pz = p[0]*uz[0] + p[1]*uz[1] + p[2]*uz[2];
|
||||
pn = p[0]*un[0] + p[1]*un[1] + p[2]*un[2];
|
||||
pw = p[0]*uw[0] + p[1]*uw[1] + p[2]*uw[2];
|
||||
|
||||
proj = sqrt(pn*pn + pw*pw);
|
||||
az = 0.0;
|
||||
if (proj > 0.0) {
|
||||
az = -atan2(pw, pn) * RAD2DEG;
|
||||
if (az < 0)
|
||||
az += 360;
|
||||
if (az >= 360)
|
||||
az -= 360;
|
||||
}
|
||||
zd = atan2(proj, pz) * RAD2DEG;
|
||||
hor.ra = ra;
|
||||
hor.dec = dec;
|
||||
|
||||
if (refraction == REFRACTION_NORMAL || refraction == REFRACTION_JPLHOR)
|
||||
{
|
||||
double zd0, refr, hd;
|
||||
int j;
|
||||
|
||||
zd0 = zd;
|
||||
|
||||
// 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 = 90.0 - zd;
|
||||
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;
|
||||
}
|
||||
|
||||
zd -= refr;
|
||||
|
||||
if (refr > 0.0 && zd > 3.0e-4)
|
||||
{
|
||||
double sinzd = sin(zd * DEG2RAD);
|
||||
double coszd = cos(zd * DEG2RAD);
|
||||
double sinzd0 = sin(zd0 * DEG2RAD);
|
||||
double coszd0 = cos(zd0 * DEG2RAD);
|
||||
double pr[3];
|
||||
|
||||
for (j=0; j<3; ++j)
|
||||
pr[j] = ((p[j] - coszd0 * uz[j]) / sinzd0)*sinzd + uz[j]*coszd;
|
||||
|
||||
proj = sqrt(pr[0]*pr[0] + pr[1]*pr[1]);
|
||||
if (proj > 0)
|
||||
{
|
||||
hor.ra = atan2(pr[1], pr[0]) * RAD2DEG / 15;
|
||||
if (hor.ra < 0)
|
||||
hor.ra += 24;
|
||||
if (hor.ra >= 24)
|
||||
hor.ra -= 24;
|
||||
}
|
||||
else
|
||||
{
|
||||
hor.ra = 0;
|
||||
}
|
||||
hor.dec = atan2(pr[2], proj) * RAD2DEG;
|
||||
}
|
||||
}
|
||||
|
||||
hor.azimuth = az;
|
||||
hor.altitude = 90.0 - zd;
|
||||
return hor;
|
||||
}
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
|
||||
@@ -1326,7 +1326,7 @@ function ter2cel(time, vec1) {
|
||||
* @param {(Date|number|Astronomy.AstroTime)} date
|
||||
* The date and time for which to find horizontal coordinates.
|
||||
*
|
||||
* @param {Astronomy.Observer} location
|
||||
* @param {Astronomy.Observer} observer
|
||||
* The location of the observer for which to find horizontal coordinates.
|
||||
*
|
||||
* @param {number} ra
|
||||
@@ -1355,13 +1355,13 @@ function ter2cel(time, vec1) {
|
||||
*
|
||||
* @returns {Astronomy.HorizontalCoordinates}
|
||||
*/
|
||||
Astronomy.Horizon = function(date, location, ra, dec, refraction) { // based on NOVAS equ2hor()
|
||||
Astronomy.Horizon = function(date, observer, ra, dec, refraction) { // based on NOVAS equ2hor()
|
||||
let time = Astronomy.MakeTime(date);
|
||||
|
||||
const sinlat = Math.sin(location.latitude * DEG2RAD);
|
||||
const coslat = Math.cos(location.latitude * DEG2RAD);
|
||||
const sinlon = Math.sin(location.longitude * DEG2RAD);
|
||||
const coslon = Math.cos(location.longitude * DEG2RAD);
|
||||
const sinlat = Math.sin(observer.latitude * DEG2RAD);
|
||||
const coslat = Math.cos(observer.latitude * DEG2RAD);
|
||||
const sinlon = Math.sin(observer.longitude * DEG2RAD);
|
||||
const coslon = Math.cos(observer.longitude * DEG2RAD);
|
||||
const sindc = Math.sin(dec * DEG2RAD);
|
||||
const cosdc = Math.cos(dec * DEG2RAD);
|
||||
const sinra = Math.sin(ra * 15 * DEG2RAD);
|
||||
@@ -1392,7 +1392,7 @@ Astronomy.Horizon = function(date, location, ra, dec, refraction) { // based
|
||||
let out_dec = dec;
|
||||
|
||||
if (refraction) {
|
||||
let zd1, refr, j;
|
||||
let refr, j;
|
||||
let zd0 = zd;
|
||||
|
||||
if (refraction === 'normal' || refraction === 'jplhor') {
|
||||
|
||||
@@ -730,6 +730,31 @@ static void geo_pos(astro_time_t time, astro_observer_t observer, double outpos[
|
||||
precession(time.tt, pos2, 0.0, outpos);
|
||||
}
|
||||
|
||||
static void spin(double angle, const double pos1[3], double vec2[3]) {
|
||||
double angr = angle * DEG2RAD;
|
||||
double cosang = cos(angr);
|
||||
double sinang = sin(angr);
|
||||
double xx = cosang;
|
||||
double yx = sinang;
|
||||
double zx = 0;
|
||||
double xy = -sinang;
|
||||
double yy = cosang;
|
||||
double zy = 0;
|
||||
double xz = 0;
|
||||
double yz = 0;
|
||||
double zz = 1;
|
||||
|
||||
vec2[0] = xx*pos1[0] + yx*pos1[1] + zx*pos1[2];
|
||||
vec2[1] = xy*pos1[0] + yy*pos1[1] + zy*pos1[2];
|
||||
vec2[2] = xz*pos1[0] + yz*pos1[1] + zz*pos1[2];
|
||||
}
|
||||
|
||||
static void ter2cel(astro_time_t time, const double vec1[3], double vec2[3])
|
||||
{
|
||||
double gast = sidereal_time(time);
|
||||
spin(-15.0 * gast, vec1, vec2);
|
||||
}
|
||||
|
||||
/*------------------ CalcMoon ------------------*/
|
||||
|
||||
#define DECLARE_PASCAL_ARRAY_1(elemtype,name,xmin,xmax) \
|
||||
@@ -2081,8 +2106,6 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time)
|
||||
double dt;
|
||||
int iter;
|
||||
|
||||
vector.t = time;
|
||||
|
||||
switch (body)
|
||||
{
|
||||
case BODY_EARTH:
|
||||
@@ -2118,7 +2141,7 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time)
|
||||
ltime2 = Astronomy_AddDays(time, -Astronomy_VectorLength(vector) / C_AUDAY);
|
||||
dt = fabs(ltime2.tt - ltime.tt);
|
||||
if (dt < 1.0e-9)
|
||||
return vector;
|
||||
goto finished; /* Tricky: ensures we patch 'vector.t' with current time, not ante-dated time. */
|
||||
|
||||
ltime = ltime2;
|
||||
}
|
||||
@@ -2126,6 +2149,8 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time)
|
||||
break;
|
||||
}
|
||||
|
||||
finished:
|
||||
vector.t = time;
|
||||
return vector;
|
||||
}
|
||||
|
||||
@@ -2152,6 +2177,126 @@ astro_sky_t Astronomy_SkyPos(astro_vector_t gc_vector, astro_observer_t observer
|
||||
return sky;
|
||||
}
|
||||
|
||||
astro_horizon_t Astronomy_Horizon(
|
||||
astro_time_t time, astro_observer_t observer, double ra, double dec, astro_refraction_t refraction)
|
||||
{
|
||||
astro_horizon_t hor;
|
||||
double uze[3], une[3], uwe[3];
|
||||
double uz[3], un[3], uw[3];
|
||||
double p[3], pz, pn, pw, proj;
|
||||
double az, zd;
|
||||
|
||||
double sinlat = sin(observer.latitude * DEG2RAD);
|
||||
double coslat = cos(observer.latitude * DEG2RAD);
|
||||
double sinlon = sin(observer.longitude * DEG2RAD);
|
||||
double coslon = cos(observer.longitude * DEG2RAD);
|
||||
double sindc = sin(dec * DEG2RAD);
|
||||
double cosdc = cos(dec * DEG2RAD);
|
||||
double sinra = sin(ra * 15 * DEG2RAD);
|
||||
double cosra = cos(ra * 15 * 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;
|
||||
|
||||
ter2cel(time, uze, uz);
|
||||
ter2cel(time, une, un);
|
||||
ter2cel(time, uwe, uw);
|
||||
|
||||
p[0] = cosdc * cosra;
|
||||
p[1] = cosdc * sinra;
|
||||
p[2] = sindc;
|
||||
|
||||
pz = p[0]*uz[0] + p[1]*uz[1] + p[2]*uz[2];
|
||||
pn = p[0]*un[0] + p[1]*un[1] + p[2]*un[2];
|
||||
pw = p[0]*uw[0] + p[1]*uw[1] + p[2]*uw[2];
|
||||
|
||||
proj = sqrt(pn*pn + pw*pw);
|
||||
az = 0.0;
|
||||
if (proj > 0.0) {
|
||||
az = -atan2(pw, pn) * RAD2DEG;
|
||||
if (az < 0)
|
||||
az += 360;
|
||||
if (az >= 360)
|
||||
az -= 360;
|
||||
}
|
||||
zd = atan2(proj, pz) * RAD2DEG;
|
||||
hor.ra = ra;
|
||||
hor.dec = dec;
|
||||
|
||||
if (refraction == REFRACTION_NORMAL || refraction == REFRACTION_JPLHOR)
|
||||
{
|
||||
double zd0, refr, hd;
|
||||
int j;
|
||||
|
||||
zd0 = zd;
|
||||
|
||||
// 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 = 90.0 - zd;
|
||||
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;
|
||||
}
|
||||
|
||||
zd -= refr;
|
||||
|
||||
if (refr > 0.0 && zd > 3.0e-4)
|
||||
{
|
||||
double sinzd = sin(zd * DEG2RAD);
|
||||
double coszd = cos(zd * DEG2RAD);
|
||||
double sinzd0 = sin(zd0 * DEG2RAD);
|
||||
double coszd0 = cos(zd0 * DEG2RAD);
|
||||
double pr[3];
|
||||
|
||||
for (j=0; j<3; ++j) {
|
||||
pr[j] = ((p[j] - coszd0 * uz[j]) / sinzd0)*sinzd + uz[j]*coszd;
|
||||
}
|
||||
proj = sqrt(pr[0]*pr[0] + pr[1]*pr[1]);
|
||||
if (proj > 0)
|
||||
{
|
||||
hor.ra = atan2(pr[1], pr[0]) * RAD2DEG / 15;
|
||||
if (hor.ra < 0)
|
||||
hor.ra += 24;
|
||||
if (hor.ra >= 24)
|
||||
hor.ra -= 24;
|
||||
}
|
||||
else
|
||||
{
|
||||
hor.ra = 0;
|
||||
}
|
||||
hor.dec = atan2(pr[2], proj) * RAD2DEG;
|
||||
}
|
||||
}
|
||||
|
||||
hor.azimuth = az;
|
||||
hor.altitude = 90.0 - zd;
|
||||
return hor;
|
||||
}
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
|
||||
@@ -102,9 +102,27 @@ typedef struct
|
||||
}
|
||||
astro_sky_t;
|
||||
|
||||
typedef struct
|
||||
{
|
||||
double azimuth;
|
||||
double altitude;
|
||||
double ra;
|
||||
double dec;
|
||||
}
|
||||
astro_horizon_t;
|
||||
|
||||
typedef enum
|
||||
{
|
||||
REFRACTION_NONE,
|
||||
REFRACTION_NORMAL,
|
||||
REFRACTION_JPLHOR
|
||||
}
|
||||
astro_refraction_t;
|
||||
|
||||
|
||||
/*---------- functions ----------*/
|
||||
|
||||
double Astronomy_VectorLength(astro_vector_t vector);
|
||||
const char *Astronomy_BodyName(astro_body_t body);
|
||||
astro_observer_t Astronomy_MakeObserver(double latitude, double longitude, double height);
|
||||
astro_time_t Astronomy_MakeTime(int year, int month, int day, int hour, int minute, double second);
|
||||
@@ -113,7 +131,7 @@ astro_vector_t Astronomy_HelioVector(astro_body_t body, astro_time_t time);
|
||||
astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time);
|
||||
astro_vector_t Astronomy_GeoMoon(astro_time_t time);
|
||||
astro_sky_t Astronomy_SkyPos(astro_vector_t gc_vector, astro_observer_t observer);
|
||||
double Astronomy_VectorLength(astro_vector_t vector);
|
||||
astro_horizon_t Astronomy_Horizon(astro_time_t time, astro_observer_t observer, double ra, double dec, astro_refraction_t refraction);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
|
||||
@@ -482,7 +482,7 @@ function calls may be more efficient than passing in native JavaScript Date obje
|
||||
|
||||
<a name="Astronomy.Horizon"></a>
|
||||
|
||||
### Astronomy.Horizon(date, location, ra, dec, refraction) ⇒ [<code>HorizontalCoordinates</code>](#Astronomy.HorizontalCoordinates)
|
||||
### Astronomy.Horizon(date, observer, ra, dec, refraction) ⇒ [<code>HorizontalCoordinates</code>](#Astronomy.HorizontalCoordinates)
|
||||
Given a date and time, a geographic location of an observer on the Earth, and
|
||||
equatorial coordinates (right ascension and declination) of a celestial object,
|
||||
returns horizontal coordinates (azimuth and altitude angles) for that object
|
||||
@@ -493,7 +493,7 @@ as seen by that observer. Allows optional correction for atmospheric refraction.
|
||||
| Param | Type | Description |
|
||||
| --- | --- | --- |
|
||||
| date | <code>Date</code> \| <code>number</code> \| [<code>AstroTime</code>](#Astronomy.AstroTime) | The date and time for which to find horizontal coordinates. |
|
||||
| location | [<code>Observer</code>](#Astronomy.Observer) | The location of the observer for which to find horizontal coordinates. |
|
||||
| observer | [<code>Observer</code>](#Astronomy.Observer) | The location of the observer for which to find horizontal coordinates. |
|
||||
| ra | <code>number</code> | Right ascension in sidereal hours of the celestial object, referred to the mean equinox of date for the J2000 epoch. |
|
||||
| dec | <code>number</code> | Declination in degrees of the celestial object, referred to the mean equator of date for the J2000 epoch. Positive values are north of the celestial equator and negative values are south. |
|
||||
| refraction | <code>string</code> | If omitted or has a false-like value (false, null, undefined, etc.) the calculations are performed without any correction for atmospheric refraction. If the value is the string <code>"normal"</code>, uses the recommended refraction correction based on Meeus "Astronomical Algorithms" with a linear taper more than 1 degree below the horizon. The linear taper causes the refraction to linearly approach 0 as the altitude of the body approaches the nadir (-90 degrees). If the value is the string <code>"jplhor"</code>, uses a JPL Horizons compatible formula. This is the same algorithm as <code>"normal"</code>, only without linear tapering; this can result in physically impossible altitudes of less than -90 degrees, which may cause problems for some applications. (The <code>"jplhor"</code> option was created for unit testing against data generated by JPL Horizons, and is otherwise not recommended for use.) |
|
||||
|
||||
@@ -2164,7 +2164,7 @@ function ter2cel(time, vec1) {
|
||||
* @param {(Date|number|Astronomy.AstroTime)} date
|
||||
* The date and time for which to find horizontal coordinates.
|
||||
*
|
||||
* @param {Astronomy.Observer} location
|
||||
* @param {Astronomy.Observer} observer
|
||||
* The location of the observer for which to find horizontal coordinates.
|
||||
*
|
||||
* @param {number} ra
|
||||
@@ -2193,13 +2193,13 @@ function ter2cel(time, vec1) {
|
||||
*
|
||||
* @returns {Astronomy.HorizontalCoordinates}
|
||||
*/
|
||||
Astronomy.Horizon = function(date, location, ra, dec, refraction) { // based on NOVAS equ2hor()
|
||||
Astronomy.Horizon = function(date, observer, ra, dec, refraction) { // based on NOVAS equ2hor()
|
||||
let time = Astronomy.MakeTime(date);
|
||||
|
||||
const sinlat = Math.sin(location.latitude * DEG2RAD);
|
||||
const coslat = Math.cos(location.latitude * DEG2RAD);
|
||||
const sinlon = Math.sin(location.longitude * DEG2RAD);
|
||||
const coslon = Math.cos(location.longitude * DEG2RAD);
|
||||
const sinlat = Math.sin(observer.latitude * DEG2RAD);
|
||||
const coslat = Math.cos(observer.latitude * DEG2RAD);
|
||||
const sinlon = Math.sin(observer.longitude * DEG2RAD);
|
||||
const coslon = Math.cos(observer.longitude * DEG2RAD);
|
||||
const sindc = Math.sin(dec * DEG2RAD);
|
||||
const cosdc = Math.cos(dec * DEG2RAD);
|
||||
const sinra = Math.sin(ra * 15 * DEG2RAD);
|
||||
@@ -2230,7 +2230,7 @@ Astronomy.Horizon = function(date, location, ra, dec, refraction) { // based
|
||||
let out_dec = dec;
|
||||
|
||||
if (refraction) {
|
||||
let zd1, refr, j;
|
||||
let refr, j;
|
||||
let zd0 = zd;
|
||||
|
||||
if (refraction === 'normal' || refraction === 'jplhor') {
|
||||
|
||||
Reference in New Issue
Block a user