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

Unconfirmed calculation of horizontal coordinates.

Browse Source 
I have horizontal coordinates calculated, but they might
be wrong (both in how I call NOVAS functions and the JS code itself)
because I think I'm mixing up equinox of date coordinates with
J2000 coordinates for (RA,DEC).

Fixed bug that caused excessive estimate of angular error:
right ascension and azimuth are like longitudes -- they matter
less as an object approaches the poles.  Scale such longitudinal
errors by the cosine of the latitudinal counterpart.
This commit is contained in:
Don Cross
2019-04-11 14:15:31 -04:00
parent 8c392a9cc8
commit 005edb555b
5 changed files with 249 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
generate/astro_check.js
View File

@@ -2,7 +2,7 @@ var Astronomy = require('../source/js/astronomy.js');
var date = new Date('1900-01-01T00:00:00Z');
var stop = new Date('2100-01-01T00:00:00Z');
var body, pos, sky;
var body, pos, sky, hor;
const observer = Astronomy.MakeObserver(29, -81, 10);
console.log(`o ${observer.latitude} ${observer.longitude} ${observer.height}`);
@@ -16,7 +16,8 @@ while (date < stop) {
if (body !== 'Earth') {
pos = Astronomy.GeoVector(body, date);
sky = Astronomy.SkyPos(pos, observer);
console.log(`s ${body} ${pos.t.jd_tt} ${pos.t.jd_utc} ${sky.ra} ${sky.dec} ${sky.dist}`);
hor = Astronomy.Horizon(sky.t, observer, sky.ra, sky.dec);
console.log(`s ${body} ${pos.t.jd_tt} ${pos.t.jd_utc} ${sky.ra} ${sky.dec} ${sky.dist} ${hor.azimuth} ${hor.altitude}`);
}
}
}
@@ -24,7 +25,8 @@ while (date < stop) {
console.log(`v GM ${pos.t.jd_tt} ${pos.x} ${pos.y} ${pos.z}`);
sky = Astronomy.SkyPos(pos, observer);
console.log(`s GM ${pos.t.jd_tt} ${pos.t.jd_utc} ${sky.ra} ${sky.dec} ${sky.dist}`);
hor = Astronomy.Horizon(sky.t, observer, sky.ra, sky.dec);
console.log(`s GM ${pos.t.jd_tt} ${pos.t.jd_utc} ${sky.ra} ${sky.dec} ${sky.dist} ${hor.azimuth} ${hor.altitude}`);
date = new Date(date.getTime() + (24*3600*1000));
}

22
generate/debug.html Normal file
View File

@@ -0,0 +1,22 @@
<!DOCTYPE html>
<html>
<head>
<title>Astronomy debug</title>
<script src="../source/js/astronomy.js"></script>
<script>
function Run() {
var observer = Astronomy.MakeObserver(29, -81, 10);
var date = 2481207.5;
var pos = Astronomy.GeoMoon(date);
console.log(pos);
var sky = Astronomy.SkyPos(pos, observer);
console.log(sky);
var hor = Astronomy.Horizon(sky.t, observer, sky.ra, sky.dec);
console.log(hor);
}
</script>
</head>
<body>
<input type="button" onclick="Run()" value="Run" >
</body>
</html>

42
generate/generate.c
View File

@@ -139,6 +139,9 @@ static int PrintUsage(void)
"generate all\n"
" Generate astronomy source code for all target languages.\n"
"\n"
"generate check testfile\n"
" Verify the calculations in the testfile generated by a unit test.\n"
"\n"
);
return 1;
@@ -823,14 +826,15 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
{
int body, error, bodyIndex;
char name[10];
double jd_tt, jd_utc, ra, dec, dist;
double jd_tt, jd_utc, ra, dec, dist, azimuth, altitude;
double check_zenith, check_azimuth, check_altitude, check_ra, check_dec;
object obj;
sky_pos sky;
double delta_ra, delta_dec;
double delta_ra, delta_dec, delta_az, delta_alt;
*arcmin = 99999.0;
if (6 != sscanf(line, "s %9[A-Za-z] %lf %lf %lf %lf %lf", name, &jd_tt, &jd_utc, &ra, &dec, &dist))
if (8 != sscanf(line, "s %9[A-Za-z] %lf %lf %lf %lf %lf %lf %lf", name, &jd_tt, &jd_utc, &ra, &dec, &dist, &azimuth, &altitude))
{
fprintf(stderr, "CheckSkyPos: Invalid format on line %d of file %s\n", lnum, filename);
return 1;
@@ -888,12 +892,42 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const
}
delta_ra *= (15.0 * 60.0);
/* Right Ascension errors are less significant as the declination approaches the poles. */
/* For example, a 12-hour RA error for Polaris would not matter very much to its observed position. */
/* So diminish the error measurement as appropriate for this declination. */
delta_ra *= cos(sky.dec * DEG2RAD);
/* Calculate pythagorean error as if both were planar coordinates. */
*arcmin = sqrt(delta_ra*delta_ra + delta_dec*delta_dec);
if (*arcmin > 1.0)
{
fprintf(stderr, "CheckSkyPos: excessive angular error = %lf arcmin at line %d of file %s\n", *arcmin, lnum, filename);
fprintf(stderr, "CheckSkyPos: excessive (RA,DEC) error = %lf arcmin at line %d of file %s\n", *arcmin, lnum, filename);
return 1;
}
/* We have tested RA,DEC. Now measure the error in horizontal coordinates. */
equ2hor(jd_utc, jd_tt-jd_utc, 1, 0.0, 0.0,
&location->on_surf, sky.ra, sky.dec, 0,
&check_zenith, &check_azimuth, &check_ra, &check_dec);
check_altitude = 90.0 - check_zenith;
delta_az = fabs(azimuth - check_azimuth);
if (delta_az > 180.0)
delta_az = 360.0 - delta_az;
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);
/* Replace the error calculation with the sky coordinate calculation. */
*arcmin = sqrt(delta_az*delta_az + delta_alt*delta_alt);
if (*arcmin > 1.0)
{
fprintf(stderr, "CheckSkyPos: excessive (alt,az) error = %lf arcmin (alt=%lf, check=%lf) for body %d at line %d of file %s\n", *arcmin, altitude, check_altitude, body, lnum, filename);
return 1;
}

105
generate/template/astronomy.js
View File

@@ -11,6 +11,7 @@ const ERAD = 6378136.6; // mean earth radius
const AU = 1.4959787069098932e+11; // astronomical unit in meters
const ASEC2RAD = 4.848136811095359935899141e-6;
const DEG2RAD = 0.017453292519943296;
const RAD2DEG = 57.295779513082321;
const ASEC360 = 1296000;
const ANGVEL = 7.2921150e-5;
const AU_KM = 1.4959787069098932e+8;
@@ -98,26 +99,36 @@ function LeapSeconds(jd_utc) {
}
function Time(date) {
// Keep the original Date object.
this.date = date;
const MillisPerDay = 1000 * 3600 * 24;
// Convert JavaScript Date object to Julian UTC value.
this.jd_utc = (date - j2000) / (1000 * 3600 * 24) + T0;
if (date instanceof Date) {
// Keep the original Date object.
this.date = date;
// Convert UTC to Terrestrial Time (TT) by adjusting for leap seconds.
this.jd_tt = this.jd_utc + (LeapSeconds(this.jd_utc) + 32.184) / 86400.0;
// Convert JavaScript Date object to Julian UTC value.
this.jd_utc = (date - j2000) / MillisPerDay + T0;
// Convert UTC to Terrestrial Time (TT) by adjusting for leap seconds.
this.jd_tt = this.jd_utc + (LeapSeconds(this.jd_utc) + 32.184) / 86400.0;
return;
}
if (typeof date === 'number') {
this.date = new Date(j2000 - (T0 - date)*MillisPerDay);
this.jd_utc = date;
this.jd_tt = this.jd_utc + (LeapSeconds(this.jd_utc) + 32.184) / 86400.0;
return;
}
throw 'AstroTime() argument must be a Date object, a Time object, or a numeric UTC Julian date.';
}
function AstroTime(date) {
if (date instanceof Time) {
return date;
}
if (date instanceof Date) {
return new Time(date);
}
throw 'Argument must be a Date object or a Time object.';
return new Time(date);
}
var nals_t = [
@@ -831,6 +842,72 @@ function vector2radec(pos)
return { ra:ra, dec:dec, dist:dist };
}
function spin(angle, pos1) {
const angr = angle * DEG2RAD;
const cosang = Math.cos(angr);
const sinang = Math.sin(angr);
const xx = cosang;
const yx = sinang;
const zx = 0;
const xy = -sinang;
const yy = cosang;
const zy = 0;
const xz = 0;
const yz = 0;
const zz = 1;
let pos2 = [
xx*pos1[0] + yx*pos1[1] + zx*pos1[2],
xy*pos1[0] + yy*pos1[1] + zy*pos1[2],
xz*pos1[0] + yz*pos1[1] + zz*pos1[2]
];
return pos2;
}
function ter2cel(time, vec1) {
const gast = sidereal_time(time, 'gast');
let vec2 = spin(-15 * gast, vec1);
return vec2;
}
Astronomy.Horizon = function(date, location, ra, dec) { // based on NOVAS equ2hor()
let time = AstroTime(date);
// FIXFIXFIX: this code expects (ra,dec) to be in equator and equinox of date, but we are giving J2000?
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 sindc = Math.sin(dec * DEG2RAD);
const cosdc = Math.cos(dec * DEG2RAD);
const sinra = Math.sin(ra * 15 * DEG2RAD);
const cosra = Math.cos(ra * 15 * DEG2RAD);
let uze = [coslat*coslon, coslat*sinlon, sinlat];
let une = [-sinlat*coslon, -sinlat*sinlon, coslat];
let uwe = [sinlon, -coslon, 0];
let uz = ter2cel(time, uze);
let un = ter2cel(time, une);
let uw = ter2cel(time, uwe);
let p = [cosdc*cosra, cosdc*sinra, sindc];
const pz = p[0]*uz[0] + p[1]*uz[1] + p[2]*uz[2];
const pn = p[0]*un[0] + p[1]*un[1] + p[2]*un[2];
const pw = p[0]*uw[0] + p[1]*uw[1] + p[2]*uw[2];
const proj = Math.sqrt(pn*pn + pw*pw);
let az = 0;
if (proj > 0) {
az = -Math.atan2(pw, pn) * RAD2DEG;
if (az < 0) az += 360;
if (az >= 360) az -= 360;
}
let zd = Math.atan2(proj, pz) * RAD2DEG;
return { azimuth:az, altitude:90-zd };
}
function Observer(latitude_degrees, longitude_degrees, height_in_meters) {
this.latitude = latitude_degrees;
this.longitude = longitude_degrees;
@@ -848,7 +925,9 @@ Astronomy.SkyPos = function(gc_vector, observer) { // based on NOVAS place()
gc_vector.y - gc_observer[1],
gc_vector.z - gc_observer[2] ];
return vector2radec(topo_vector);
let sky = vector2radec(topo_vector);
sky.t = gc_vector.t; // patch in the Time object for convenience
return sky;
}
Astronomy.GeoMoon = function(date) {

105
source/js/astronomy.js
View File

@@ -11,6 +11,7 @@ const ERAD = 6378136.6; // mean earth radius
const AU = 1.4959787069098932e+11; // astronomical unit in meters
const ASEC2RAD = 4.848136811095359935899141e-6;
const DEG2RAD = 0.017453292519943296;
const RAD2DEG = 57.295779513082321;
const ASEC360 = 1296000;
const ANGVEL = 7.2921150e-5;
const AU_KM = 1.4959787069098932e+8;
@@ -698,26 +699,36 @@ function LeapSeconds(jd_utc) {
}
function Time(date) {
// Keep the original Date object.
this.date = date;
const MillisPerDay = 1000 * 3600 * 24;
// Convert JavaScript Date object to Julian UTC value.
this.jd_utc = (date - j2000) / (1000 * 3600 * 24) + T0;
if (date instanceof Date) {
// Keep the original Date object.
this.date = date;
// Convert UTC to Terrestrial Time (TT) by adjusting for leap seconds.
this.jd_tt = this.jd_utc + (LeapSeconds(this.jd_utc) + 32.184) / 86400.0;
// Convert JavaScript Date object to Julian UTC value.
this.jd_utc = (date - j2000) / MillisPerDay + T0;
// Convert UTC to Terrestrial Time (TT) by adjusting for leap seconds.
this.jd_tt = this.jd_utc + (LeapSeconds(this.jd_utc) + 32.184) / 86400.0;
return;
}
if (typeof date === 'number') {
this.date = new Date(j2000 - (T0 - date)*MillisPerDay);
this.jd_utc = date;
this.jd_tt = this.jd_utc + (LeapSeconds(this.jd_utc) + 32.184) / 86400.0;
return;
}
throw 'AstroTime() argument must be a Date object, a Time object, or a numeric UTC Julian date.';
}
function AstroTime(date) {
if (date instanceof Time) {
return date;
}
if (date instanceof Date) {
return new Time(date);
}
throw 'Argument must be a Date object or a Time object.';
return new Time(date);
}
var nals_t = [
@@ -1431,6 +1442,72 @@ function vector2radec(pos)
return { ra:ra, dec:dec, dist:dist };
}
function spin(angle, pos1) {
const angr = angle * DEG2RAD;
const cosang = Math.cos(angr);
const sinang = Math.sin(angr);
const xx = cosang;
const yx = sinang;
const zx = 0;
const xy = -sinang;
const yy = cosang;
const zy = 0;
const xz = 0;
const yz = 0;
const zz = 1;
let pos2 = [
xx*pos1[0] + yx*pos1[1] + zx*pos1[2],
xy*pos1[0] + yy*pos1[1] + zy*pos1[2],
xz*pos1[0] + yz*pos1[1] + zz*pos1[2]
];
return pos2;
}
function ter2cel(time, vec1) {
const gast = sidereal_time(time, 'gast');
let vec2 = spin(-15 * gast, vec1);
return vec2;
}
Astronomy.Horizon = function(date, location, ra, dec) { // based on NOVAS equ2hor()
let time = AstroTime(date);
// FIXFIXFIX: this code expects (ra,dec) to be in equator and equinox of date, but we are giving J2000?
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 sindc = Math.sin(dec * DEG2RAD);
const cosdc = Math.cos(dec * DEG2RAD);
const sinra = Math.sin(ra * 15 * DEG2RAD);
const cosra = Math.cos(ra * 15 * DEG2RAD);
let uze = [coslat*coslon, coslat*sinlon, sinlat];
let une = [-sinlat*coslon, -sinlat*sinlon, coslat];
let uwe = [sinlon, -coslon, 0];
let uz = ter2cel(time, uze);
let un = ter2cel(time, une);
let uw = ter2cel(time, uwe);
let p = [cosdc*cosra, cosdc*sinra, sindc];
const pz = p[0]*uz[0] + p[1]*uz[1] + p[2]*uz[2];
const pn = p[0]*un[0] + p[1]*un[1] + p[2]*un[2];
const pw = p[0]*uw[0] + p[1]*uw[1] + p[2]*uw[2];
const proj = Math.sqrt(pn*pn + pw*pw);
let az = 0;
if (proj > 0) {
az = -Math.atan2(pw, pn) * RAD2DEG;
if (az < 0) az += 360;
if (az >= 360) az -= 360;
}
let zd = Math.atan2(proj, pz) * RAD2DEG;
return { azimuth:az, altitude:90-zd };
}
function Observer(latitude_degrees, longitude_degrees, height_in_meters) {
this.latitude = latitude_degrees;
this.longitude = longitude_degrees;
@@ -1448,7 +1525,9 @@ Astronomy.SkyPos = function(gc_vector, observer) { // based on NOVAS place()
gc_vector.y - gc_observer[1],
gc_vector.z - gc_observer[2] ];
return vector2radec(topo_vector);
let sky = vector2radec(topo_vector);
sky.t = gc_vector.t; // patch in the Time object for convenience
return sky;
}
Astronomy.GeoMoon = function(date) {