From 005edb555bbd5c268ca4d34cba89fe7909774ac7 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Thu, 11 Apr 2019 14:15:31 -0400 Subject: [PATCH] Unconfirmed calculation of horizontal coordinates. 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. --- generate/astro_check.js | 8 ++- generate/debug.html | 22 +++++++ generate/generate.c | 42 +++++++++++-- generate/template/astronomy.js | 105 +++++++++++++++++++++++++++++---- source/js/astronomy.js | 105 +++++++++++++++++++++++++++++---- 5 files changed, 249 insertions(+), 33 deletions(-) create mode 100644 generate/debug.html diff --git a/generate/astro_check.js b/generate/astro_check.js index e16d75d4..9a35b344 100644 --- a/generate/astro_check.js +++ b/generate/astro_check.js @@ -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)); } diff --git a/generate/debug.html b/generate/debug.html new file mode 100644 index 00000000..0f340d42 --- /dev/null +++ b/generate/debug.html @@ -0,0 +1,22 @@ + + + + Astronomy debug + + + + + + + diff --git a/generate/generate.c b/generate/generate.c index 811330b1..bbf6797b 100644 --- a/generate/generate.c +++ b/generate/generate.c @@ -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; } diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index 5184ecb3..b657e6fb 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -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) { diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 7b9a0e41..2fc88cfd 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -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) {