'use strict';
const fs = require('fs');
const Astronomy = require('../source/js/astronomy.min.js');
const MINUTES_PER_DAY = 1440;
const SECONDS_PER_DAY = 86400;
let Verbose = false;

function Debug(text) {
    if (Verbose) console.debug(text);
}

function Fail(message) {
    console.trace();
    console.log(`FATAL(test.js): ${message}`);
    process.exit(1);
}

function Pass(name) {
    console.log(`JS ${name}: PASS`);
    return 0;
}

function v(x) {
    // Verify that 'x' is really a number.

    if (!Number.isFinite(x)) {
        console.trace();
        throw `Not a finite number: ${x}`;
    }

    return x;
}

function int(s) {
    return v(parseInt(s));
}

function float(s) {
    return v(parseFloat(s));
}

function abs(x) {
    return Math.abs(v(x));
}

function max(a, b) {
    return Math.max(v(a), v(b));
}

function min(a, b) {
    return Math.min(v(a), v(b));
}

function sin(x) {
    return Math.sin(v(x));
}

function cos(x) {
    return Math.cos(v(x));
}

function sqrt(x) {
    return v(Math.sqrt(v(x)));
}

function ReadLines(filename) {
    // Load the entire text of the given file.
    const text = fs.readFileSync(filename, {encoding:'utf8'});

    // Split it into lines, handling the different line endings
    // in Linux and Windows.
    // FIXFIXFIX: This might not work on Mac OS or other operating systems.
    const lines = text.trimEnd().split(/\r?\n/);
    return lines;
}

function SelectJupiterMoon(jm, mindex) {
    return [jm.io, jm.europa, jm.ganymede, jm.callisto][mindex] ||
        Fail(`SelectJupiterMoon: invalid mindex = ${mindex}`);
}

function AstroCheck() {
    var date = Astronomy.MakeTime(new Date('1700-01-01T00:00:00Z'));
    var stop = Astronomy.MakeTime(new Date('2200-01-01T00:00:00Z'));
    var body, pos, hor, dt, j2000, ofdate, time;
    const observer = new Astronomy.Observer(29, -81, 10);
    const bodylist = [
        Astronomy.Body.Sun, Astronomy.Body.Mercury, Astronomy.Body.Venus, Astronomy.Body.Earth, Astronomy.Body.Mars,
        Astronomy.Body.Jupiter, Astronomy.Body.Saturn, Astronomy.Body.Uranus, Astronomy.Body.Neptune, Astronomy.Body.Pluto,
        Astronomy.Body.SSB, Astronomy.Body.EMB
    ];

    console.log(`o ${observer.latitude.toFixed(6)} ${observer.longitude.toFixed(6)} ${observer.height.toFixed(6)}`);

    dt = (10 + Math.PI/100);       // 10.03141592... days; exercise different times of day
    while (date.tt < stop.tt) {
        time = Astronomy.MakeTime(date);
        for (body of bodylist) {
            pos = Astronomy.HelioVector(body, date);
            console.log(`v ${body} ${pos.t.tt.toExponential(18)} ${pos.x.toExponential(18)} ${pos.y.toExponential(18)} ${pos.z.toExponential(18)}`);

            if (body !== 'Earth' && body !== 'EMB' && body !== 'SSB') {
                j2000 = Astronomy.Equator(body, date, observer, false, false);
                ofdate = Astronomy.Equator(body, date, observer, true, true);
                hor = Astronomy.Horizon(date, observer, ofdate.ra, ofdate.dec);
                console.log(`s ${body} ${time.tt.toExponential(18)} ${time.ut.toExponential(18)} ${j2000.ra.toExponential(18)} ${j2000.dec.toExponential(18)} ${j2000.dist.toExponential(18)} ${hor.azimuth.toExponential(18)} ${hor.altitude.toExponential(18)}`);
            }
        }
        pos = Astronomy.GeoMoon(date);
        console.log(`v GM ${pos.t.tt.toExponential(18)} ${pos.x.toExponential(18)} ${pos.y.toExponential(18)} ${pos.z.toExponential(18)}`);

        j2000 = Astronomy.Equator('Moon', date, observer, false, false);
        ofdate = Astronomy.Equator('Moon', date, observer, true, true);
        hor = Astronomy.Horizon(date, observer, ofdate.ra, ofdate.dec);
        console.log(`s GM ${time.tt.toExponential(18)} ${time.ut.toExponential(18)} ${j2000.ra.toExponential(18)} ${j2000.dec.toExponential(18)} ${j2000.dist.toExponential(18)} ${hor.azimuth.toExponential(18)} ${hor.altitude.toExponential(18)}`);

        const jm = Astronomy.JupiterMoons(time);
        for (let mindex = 0; mindex < 4; ++mindex) {
            const moon = SelectJupiterMoon(jm, mindex);
            console.log(`j ${mindex} ${time.tt.toExponential(18)} ${time.ut.toExponential(18)} ${moon.x.toExponential(18)} ${moon.y.toExponential(18)} ${moon.z.toExponential(18)} ${moon.vx.toExponential(18)} ${moon.vy.toExponential(18)} ${moon.vz.toExponential(18)}`);
        }

        // Nutation calculations.
        const et = Astronomy.e_tilt(time);
        console.log(`n ${et.dpsi.toExponential(18)} ${et.deps.toExponential(18)}`);

        const sphere = Astronomy.EclipticGeoMoon(time);
        console.log(`m ${sphere.lat.toExponential(18)} ${sphere.lon.toExponential(18)} ${sphere.dist.toExponential(18)}`);

        date = date.AddDays(dt);
    }

    return 0;
}


function MoonPhase() {
    function LoadMoonPhaseData(filename) {
        // Load known moon phase times from US Naval Observatory.
        const lines = ReadLines(filename);
        let data = [];
        for (let row of lines) {
            let token = row.split(' ');
            data.push({quarter:int(token[0]), date:new Date(token[1])});
        }
        return data;
    }

    function TestLongitudes(data) {
        let max_arcmin = 0;
        for (let row of data) {
            let elong = v(Astronomy.MoonPhase(row.date));
            let expected_elong = 90 * v(row.quarter);
            let degree_error = abs(elong - expected_elong);
            if (degree_error > 180) degree_error = 360 - degree_error;
            let arcmin = 60 * degree_error;
            max_arcmin = max(max_arcmin, arcmin);
        }
        console.log(`JS TestLongitudes: nrows = ${data.length}, max_arcmin = ${max_arcmin}`);
        if (max_arcmin > 1.0) {
            console.log(`JS TestLongitudes: EXCESSIVE ANGULAR ERROR`);
            return 1;
        }
        return 0;
    }

    function SearchYear(year, data, index) {
        const millis_per_minute = 60*1000;
        const threshold_minutes = 1.5;    // max tolerable prediction error in minutes
        let date = new Date(Date.UTC(year, 0, 1));
        let maxdiff = 0;
        let count = 0;
        let mq = Astronomy.SearchMoonQuarter(date);
        while (index < data.length && data[index].date.getUTCFullYear() === year) {

            // Verify that we found anything at all.
            if (!mq) {
                console.log(`JS SearchYear: could not find next moon quarter after ${date.toISOString()}`);
                return 1;
            }

            // Verify that we found the correct quarter.
            if (data[index].quarter !== mq.quarter) {
                console.log(`JS SearchYear: predicted quarter ${mq.quarter} but correct is ${data[index].quarter} for ${date.toISOString()}`);
                return 1;
            }

            // Verify that the date and time we found is very close to the correct answer.
            // Calculate the discrepancy in minutes.
            // This is appropriate because the "correct" answers are only given to the minute.
            let diff = abs(mq.time.date - data[index].date) / millis_per_minute;
            if (diff > threshold_minutes) {
                console.log(`JS SearchYear: EXCESSIVE ERROR = ${diff.toFixed(3)} minutes, correct=${data[index].date.toISOString()}, calculated=${mq.time.toString()}`);
                return 1;
            }
            maxdiff = max(maxdiff, diff);

            ++index;
            ++count;
            mq = Astronomy.NextMoonQuarter(mq);
            date = mq.time.date;
        }
        Debug(`JS SearchYear(${year}): count=${count}, maxdiff=${maxdiff.toFixed(3)}`);
        return 0;
    }

    function TestSearch(data) {
        // Search each year finding each quarter moon phase and confirm
        // they match up with the correct answers stored in the 'data' array.
        let index = 0;  // index into 'data' for the current year we are interested in.
        for (let year=1800; year <= 2100; year += 10) {
            while (data[index].date.getUTCFullYear() < year) ++index;
            let error = SearchYear(year, data, index);
            if (error) return error;
        }
        return 0;
    }

    function Statistics(data) {
        // Figure out mean, min, and max differences between various moon quarters.
        const stats = [null, null, null, null];
        let prev = null;
        for (let curr of data) {
            if (prev && curr.date.getUTCFullYear() === prev.date.getUTCFullYear()) {
                let dt = (curr.date - prev.date) / (SECONDS_PER_DAY * 1000);
                let s = stats[prev.quarter];
                if (s) {
                    s.min = min(s.min, dt);
                    s.max = max(s.max, dt);
                    s.sum += dt;
                    ++s.count;
                } else {
                    stats[prev.quarter] = {
                        min: dt,
                        max: dt,
                        sum: dt,
                        count: 1
                    };
                }
            }
            prev = curr;
        }

        for (let q=0; q < 4; ++q) {
            let s = stats[q];
            if (s) {
                console.log(`JS MoonPhase statistics: q=${q} min=${s.min.toFixed(3)} avg=${(s.sum/s.count).toFixed(3)} max=${s.max.toFixed(3)} span=${(s.max-s.min).toFixed(3)}`);
            }
        }
    }

    const TestData = LoadMoonPhaseData('moonphase/moonphases.txt');
    if (Verbose) {
        Statistics(TestData);
    }
    if (TestLongitudes(TestData)) return 1;
    if (TestSearch(TestData)) return 1;
    console.log('JS MoonPhase: PASS');
    return 0;
}

function MoonReverse(longitude) {
    return (
        MoonReversePhase(0) ||
        MoonReversePhase(90) ||
        MoonReversePhase(180) ||
        MoonReversePhase(270)
    );
}

function MoonReversePhase(longitude) {
    // Verify that SearchMoonPhase works both forward and backward in time.

    const nphases = 5000;
    let utList = [];
    let i, result, diff;

    let dtMin = +1000;
    let dtMax = -1000;
    let time = Astronomy.MakeTime(new Date('1800-01-01T00:00:00Z'));
    for (i = 0; i < nphases; ++i) {
        result = Astronomy.SearchMoonPhase(longitude, time, +40);
        if (result === null) {
            console.error(`JS MoonReverse(i=${i}): failed to find phase ${longitude} after ${time}`);
            return 1;
        }
        utList.push(result.ut);
        if (i > 0) {
            // Verify that consecutive phase events are reasonably close to the synodic period (29.5 days) apart.
            const dt = v(utList[i] - utList[i-1]);
            if (dt < dtMin) dtMin = dt;
            if (dt > dtMax) dtMax = dt;
        }
        time = result.AddDays(+0.1);
    }

    Debug(`JS MoonReverse(${longitude}): dtMin=${dtMin.toFixed(3)} days, dtMax=${dtMax.toFixed(3)} days.`);
    if (dtMin < 29.175 || dtMax > 29.926) {
        console.error(`JS MoonReverse(${longitude}): Time between consecutive phases is suspicious.`);
        return 1;
    }

    // Do a reverse chronological search and make sure the results are consistent with the forward search.
    time = time.AddDays(20);
    let maxDiff = 0;
    for (i = nphases-1; i >= 0; --i) {
        result = Astronomy.SearchMoonPhase(longitude, time, -40);
        if (result === null) {
            console.error(`JS MoonReverse(i=${i}): failed to find phase ${longitude} before ${time}`);
            return 1;
        }
        diff = SECONDS_PER_DAY * abs(result.ut - utList[i]);
        if (diff > maxDiff) maxDiff = diff;
        time = result.AddDays(-0.1);
    }

    Debug(`JS MoonReverse(${longitude}): Maximum discrepancy in reverse search = ${maxDiff.toFixed(3)} seconds.`);
    if (maxDiff > 0.165) {
        console.error(`JS MoonReverse(${longitude}): EXCESSIVE DISCREPANCY in reverse search.`);
        return 1;
    }

    // Pick a pair of consecutive events from the middle of the list.
    // Verify forward and backward searches work correctly from many intermediate times.
    const nslots = 100;
    const k = Math.floor(nslots / 2);
    const ut1 = utList[k];
    const ut2 = utList[k+1];
    for (i = 1; i < nslots; ++i) {
        const ut = ut1 + (i/nslots)*(ut2 - ut1);
        time = Astronomy.MakeTime(ut);

        const before = Astronomy.SearchMoonPhase(longitude, time, -40);
        if (!before) {
            console.error(`JS MoonReverse(${longitude}): backward search failed from ${time}`);
            return 1;
        }
        diff = SECONDS_PER_DAY * abs(before.ut - ut1);
        if (diff > 0.07) {
            console.error(`JS MoonPhase(${longitude}): backward search from ${time} has error = ${diff.toExponential(4)} seconds.`);
            return 1;
        }

        const after = Astronomy.SearchMoonPhase(longitude, time, +40);
        if (!after) {
            console.error(`JS MoonReverse(${longitude}): forward search failed from ${time}`);
            return 1;
        }
        diff = SECONDS_PER_DAY * abs(after.ut - ut2);
        if (diff > 0.07) {
            console.error(`JS MoonPhase(${longitude}): forward search from ${time} has error = ${diff.toExponential(4)} seconds.`);
            return 1;
        }
    }

    console.log(`JS MoonReverse(${longitude}): PASS`);
    return 0;
}


function LunarApsis() {
    const filename = 'apsides/moon.txt';
    const lines = ReadLines(filename);

    const time_before = new Date();
    let evt = Astronomy.SearchLunarApsis(new Date(Date.UTC(2001, 0, 1)));

    let lnum = 0;
    let count = 0;
    let max_minute_error = 0;
    let max_dist_error = 0;
    for (let line of lines) {
        ++lnum;

        let token = line.split(/\s+/);
        if (token.length !== 3)
            continue;

        let kind = int(token[0]);
        let date = new Date(token[1]);
        let dist = int(token[2]);

        if (evt.kind !== kind) {
            console.log('line = ', line);
            throw `${filename} line ${lnum}: Expected apsis type ${kind}, found ${evt.kind}`;
        }

        let diff_minutes = abs(evt.time.date - date) / (1000 * 60);
        if (diff_minutes > 35) {
            throw `${filename} line ${lnum}: Excessive time error: ${diff_minutes} minutes`;
        }
        max_minute_error = max(max_minute_error, diff_minutes);

        let diff_dist = abs(evt.dist_km - dist);
        if (diff_dist > 25) {
            throw `${filename} line ${lnum}: Excessive distance error: ${diff_dist} km`;
        }
        max_dist_error = max(max_dist_error, diff_dist);

        ++count;
        evt = Astronomy.NextLunarApsis(evt);
    }
    const time_after = new Date();
    const elapsed = (time_after - time_before) / 1000;

    Debug(`JS LunarApsis: verified ${count} lines, max time error = ${max_minute_error.toFixed(3)} minutes, max dist error = ${max_dist_error.toFixed(3)} km.`);

    if (count !== 2651)
        throw 'FATAL: Did not process the expected number of data rows!';

    console.log(`JS LunarApsis PASS: time=${elapsed.toFixed(3)}`);
    return 0;
}


function LunarEclipseIssue78() {
    // https://github.com/cosinekitty/astronomy/issues/78

    let eclipse = Astronomy.SearchLunarEclipse(new Date(Date.UTC(2020, 11, 19)));
    const expected_peak = new Date('2021-05-26T11:18:42Z');  // https://www.timeanddate.com/eclipse/lunar/2021-may-26
    const dt = (expected_peak - eclipse.peak.date) / 1000;
    if (abs(dt) > 40.0)
        throw `LunarEclipseIssue78: Excessive prediction error = ${dt} seconds.`;
    if (eclipse.kind !== Astronomy.EclipseKind.Total)
        throw `Expected total eclipse; found: ${eclipse.kind}`;
    console.log(`JS LunarEclipseIssue78: PASS`);
    return 0;
}


function LunarEclipse() {
    Astronomy.CalcMoonCount = 0;
    const filename = 'eclipse/lunar_eclipse.txt';
    const lines = ReadLines(filename);
    const diff_limit = 2.0;
    let lnum = 0;
    let skip_count = 0;
    let max_diff_minutes = 0;
    let sum_diff_minutes = 0;
    let diff_count = 0;
    let eclipse = Astronomy.SearchLunarEclipse(new Date(Date.UTC(1701, 0)));
    for (let line of lines) {
        ++lnum;

        // Make sure numeric data are finite values.
        v(eclipse.obscuration);
        v(eclipse.sd_partial);
        v(eclipse.sd_penum);
        v(eclipse.sd_total);

        if (line.length < 17) {
            console.error(`JS LunarEclipse(${filename} line ${lnum}): line is too short.`);
            return 1;
        }
        const time_text = line.substr(0, 17);
        const peak_time = Astronomy.MakeTime(new Date(time_text));
        const token = line.substr(17).trim().split(/\s+/);
        if (token.length !== 2) {
            console.error(`JS LunarEclipse(${filename} line ${lnum}): invalid number of tokens.`);
            return 1;
        }
        const partial_minutes = float(token[0]);
        const total_minutes = float(token[1]);

        // Verify that the calculated eclipse semi-durations are consistent with the kind.
        // Verify that obscurations also make sense for the kind.
        let sd_valid = false;
        let frac_valid = false;
        switch (eclipse.kind) {
        case Astronomy.EclipseKind.Penumbral:
            sd_valid = (eclipse.sd_penum > 0.0) && (eclipse.sd_partial === 0.0) && (eclipse.sd_total === 0.0);
            frac_valid = (eclipse.obscuration === 0.0);
            break;

        case Astronomy.EclipseKind.Partial:
            sd_valid = (eclipse.sd_penum > 0.0) && (eclipse.sd_partial > 0.0) && (eclipse.sd_total === 0.0);
            frac_valid = (eclipse.obscuration > 0.0) && (eclipse.obscuration < 1.0);
            break;

        case Astronomy.EclipseKind.Total:
            sd_valid = (eclipse.sd_penum > 0.0) && (eclipse.sd_partial > 0.0) && (eclipse.sd_total > 0.0);
            frac_valid = (eclipse.obscuration === 1.0);
            break;

        default:
            console.error(`JS LunarEclipse(${filename} line ${lnum}): invalid eclipse kind ${eclipse.kind}.`);
            return 1;
        }

        if (!sd_valid) {
            console.error(`JS LunarEclipse(${filename} line ${lnum}): invalid semidurations for kind ${kind}: penum=${sd_penum}, partial=${sd_partial}, total=${sd_total}.`);
            return 1;
        }

        if (!frac_valid) {
            console.error(`JS LunarEclipse(${filename} line ${lnum}): invalid obscuration ${eclipse.obscuration} for kind ${kind}.`);
            return 1;
        }

        // Check eclipse peak.
        let diff_days = v(eclipse.peak.ut) - v(peak_time.ut);

        // Tolerate missing penumbral eclipses - skip to next input line without calculating next eclipse.
        if (partial_minutes == 0.0 && diff_days > 20.0) {
            ++skip_count;
            continue;
        }

        let diff_minutes = (24.0 * 60.0) * abs(diff_days);
        sum_diff_minutes += diff_minutes;
        ++diff_count;

        if (diff_minutes > diff_limit) {
            console.error(`JS LunarEclipse expected peak: ${peak_time}`);
            console.error(`JS LunarEclipse found    peak: ${eclipse.peak}`);
            console.error(`JS LunarEclipse(${filename} line ${lnum}): EXCESSIVE peak time error = ${diff_minutes} minutes (${diff_days} days).`);
            return 1;
        }

        if (diff_minutes > max_diff_minutes)
            max_diff_minutes = diff_minutes;

        /* check partial eclipse duration */

        diff_minutes = abs(partial_minutes - eclipse.sd_partial);
        sum_diff_minutes += diff_minutes;
        ++diff_count;

        if (diff_minutes > diff_limit) {
            console.error(`JS LunarEclipse(${filename} line ${lnum}): EXCESSIVE partial eclipse semiduration error: ${diff_minutes} minutes.`);
            return 1;
        }

        if (diff_minutes > max_diff_minutes)
            max_diff_minutes = diff_minutes;

        /* check total eclipse duration */

        diff_minutes = abs(total_minutes - eclipse.sd_total);
        sum_diff_minutes += diff_minutes;
        ++diff_count;

        if (diff_minutes > diff_limit) {
            console.error(`JS LunarEclipse(${filename} line ${lnum}): EXCESSIVE total eclipse semiduration error: ${diff_minutes} minutes.`);
            return 1;
        }

        if (diff_minutes > max_diff_minutes)
            max_diff_minutes = diff_minutes;

        /* calculate for next iteration */

        eclipse = Astronomy.NextLunarEclipse(eclipse.peak);
    }
    console.log(`JS LunarEclipse: PASS (verified ${lnum}, skipped ${skip_count}, max_diff_minutes = ${max_diff_minutes}, avg_diff_minutes = ${sum_diff_minutes / diff_count}, moon calcs = ${Astronomy.CalcMoonCount})`);
    return 0;
}

function VectorFromAngles(lat, lon) {
    const latrad = Astronomy.DEG2RAD * lat;
    const lonrad = Astronomy.DEG2RAD * lon;
    const coslat = cos(latrad);
    return [
        cos(lonrad) * coslat,
        sin(lonrad) * coslat,
        sin(latrad)
    ];
}


function AngleDiff(alat, alon, blat, blon) {
    const a = VectorFromAngles(alat, alon);
    const b = VectorFromAngles(blat, blon);
    const dot = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
    if (dot <= -1.0) {
        return 180.0;
    }
    if (dot >= +1.0) {
        return 0.0;
    }
    return v(Astronomy.RAD2DEG * Math.acos(dot));
}


function GlobalSolarEclipse() {
    const expected_count = 1180;
    const filename = 'eclipse/solar_eclipse.txt';
    const lines = ReadLines(filename);
    let max_minutes = 0.0;
    let max_angle = 0.0;
    let skip_count = 0;
    let eclipse = Astronomy.SearchGlobalSolarEclipse(new Date(Date.UTC(1701, 0)));
    let lnum = 0;
    for (let line of lines) {
        ++lnum;
        // 1889-12-22T12:54:15Z   -6 T   -12.7   -12.8
        let token = line.trim().split(/\s+/);
        if (token.length !== 5) {
            console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): invalid token count = ${token.length}`);
            return 1;
        }
        const peak = Astronomy.MakeTime(new Date(token[0]));
        const typeChar = token[2];
        const lat = float(token[3]);
        const lon = float(token[4]);
        const expected_kind = {
            'P': Astronomy.EclipseKind.Partial,
            'A': Astronomy.EclipseKind.Annular,
            'T': Astronomy.EclipseKind.Total,
            'H': Astronomy.EclipseKind.Total
        }[typeChar];

        let diff_days = eclipse.peak.tt - peak.tt;
        // Sometimes we find marginal eclipses that aren't listed in the test data.
        // Ignore them if the distance between the Sun/Moon shadow axis and the Earth's center is large.
        while (diff_days < -25.0 && eclipse.distance > 9000.0) {
            ++skip_count;
            eclipse = Astronomy.NextGlobalSolarEclipse(eclipse.peak);
            diff_days = eclipse.peak.ut - peak.ut;
        }

        // Validate the eclipse prediction.
        const diff_minutes = (24 * 60) * abs(diff_days);
        if (diff_minutes > 7.56) {
            console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): EXCESSIVE TIME ERROR = ${diff_minutes} minutes`);
            return 1;
        }

        if (diff_minutes > max_minutes) {
            max_minutes = diff_minutes;
        }

        // Validate the eclipse kind, but only when it is not a "glancing" eclipse.
        if ((eclipse.distance < 6360) && (eclipse.kind != expected_kind)) {
            console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): WRONG ECLIPSE KIND: expected ${expected_kind}, found ${eclipse.kind}`);
            return 1;
        }

        if (eclipse.kind === Astronomy.EclipseKind.Total || eclipse.kind === Astronomy.EclipseKind.Annular) {
            // When the distance between the Moon's shadow ray and the Earth's center is beyond 6100 km,
            // it creates a glancing blow whose geographic coordinates are excessively sensitive to
            // slight changes in the ray. Therefore, it is unreasonable to count large errors there.
            if (eclipse.distance < 6100.0) {
                const diff_angle = AngleDiff(lat, lon, eclipse.latitude, eclipse.longitude);
                if (diff_angle > 0.247) {
                    console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): EXCESSIVE GEOGRAPHIC LOCATION ERROR = ${diff_angle} degrees`);
                    console.log(`   calculated lat=${eclipse.latitude}, lon=${eclipse.longitude}; expected ${lat}, ${lon}`);
                    return 1;
                }
                if (diff_angle > max_angle) {
                    max_angle = diff_angle;
                }
            }
        }

        // Verify the obscuration value is consistent with the eclipse kind.
        switch (eclipse.kind) {
            case Astronomy.EclipseKind.Partial:
                if (eclipse.obscuration !== undefined) {
                    console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): Expected NAN obscuration for partial eclipse, but found ${eclipse.obscuration}`);
                    return 1;
                }
                break;

            case Astronomy.EclipseKind.Annular:
                if (!Number.isFinite(eclipse.obscuration)) {
                    console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): Expected finite obscuration for annular eclipse.`);
                    return 1;
                }
                if (eclipse.obscuration < 0.8 || eclipse.obscuration >= 1.0) {
                    console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): Invalid obscuration = ${eclipse.obscuration} for annular eclipse.`);
                    return 1;
                }
                break;

            case Astronomy.EclipseKind.Total:
                if (!Number.isFinite(eclipse.obscuration)) {
                    console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): Expected finite obscuration for total eclipse.`);
                    return 1;
                }
                if (eclipse.obscuration !== 1.0) {
                    console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): Invalid obscuration = ${eclipse.obscuration} for total eclipse.`);
                    return 1;
                }
                break;

            default:
                console.error(`JS GlobalSolarEclipse(${filename} line ${lnum}): Unhandled eclipse kind ${eclipse.kind} for obscuration check.`);
                return 1;
        }

        eclipse = Astronomy.NextGlobalSolarEclipse(eclipse.peak);
    }

    if (lnum != expected_count) {
        console.error(`JS GlobalSolarEclipse: WRONG LINE COUNT = ${lnum}, expected ${expected_count}`);
        return 1;
    }

    if (skip_count > 2) {
        console.error(`JS GlobalSolarEclipse: EXCESSIVE SKIP COUNT = ${skip_count}`);
        return 1;
    }

    console.log(`JS GlobalSolarEclipse: PASS (${lnum} verified, ${skip_count} skipped, max minutes = ${max_minutes}, max angle = ${max_angle})`);
    return 0;
}


function LocalSolarEclipse1() {
    // Re-use the test data for global solar eclipses, only feed the given coordinates
    // into the local solar eclipse predictor as the observer's location.
    // In each case, start the search 20 days before the expected eclipse.
    // Then verify that the peak time and eclipse type is correct in each case.

    const expected_count = 1180;
    const filename = 'eclipse/solar_eclipse.txt';
    const lines = ReadLines(filename);
    let max_minutes = 0.0;
    let skip_count = 0;
    let lnum = 0;
    for (let line of lines) {
        ++lnum;
        // 1889-12-22T12:54:15Z   -6 T   -12.7   -12.8
        let token = line.trim().split(/\s+/);
        if (token.length !== 5) {
            console.error(`JS LocalSolarEclipse1(${filename} line ${lnum}): invalid token count = ${token.length}`);
            return 1;
        }
        const peak = Astronomy.MakeTime(new Date(token[0]));
        //const typeChar = token[2];
        const lat = float(token[3]);
        const lon = float(token[4]);
        const observer = new Astronomy.Observer(lat, lon, 0);

        // Start the search 20 days before we know the eclipse should peak.
        const search_start = peak.AddDays(-20);
        const eclipse = Astronomy.SearchLocalSolarEclipse(search_start, observer);

        // Validate the predicted eclipse peak time.
        const diff_days = v(eclipse.peak.time.tt) - v(peak.tt);
        if (diff_days > 20) {
            ++skip_count;
            continue;
        }

        const diff_minutes = (24 * 60) * abs(diff_days);
        if (diff_minutes > 7.737) {
            console.error(`JS LocalSolarEclipse1(${filename} line ${lnum}): EXCESSIVE TIME ERROR = ${diff_minutes} minutes`);
            return 1;
        }

        if (diff_minutes > max_minutes) {
            max_minutes = diff_minutes;
        }

        // Verify obscuration makes sense for this kind of eclipse.
        if (!Number.isFinite(eclipse.obscuration)) {
            console.error(`JS LocalSolarEclipse1(${filename} line ${lnum}) FAIL: obscuration is not finite.`);
            return 1;
        }

        switch (eclipse.kind) {
            case Astronomy.EclipseKind.Annular:
            case Astronomy.EclipseKind.Partial:
                if (eclipse.obscuration <= 0.0 || eclipse.obscuration >= 1.0) {
                    console.error(`JS LocalSolarEclipse1(${filename} line ${lnum}) FAIL: ${eclipse.kind} eclipse has invalid obscuration ${obscuration}`);
                    return 1;
                }
                break;

            case Astronomy.EclipseKind.Total:
                if (eclipse.obscuration !== 1.0) {
                    console.error(`JS LocalSolarEclipse1(${filename} line ${lnum}) FAIL: total eclipse has invalid obscuration ${obscuration}`);
                    return 1;
                }
                break;

            default:
                console.error(`JS LocalSolarEclipse1(${filename} line ${lnum}) FAIL: invalid eclipse kind ${eclipse.kind}`);
                return 1;
        }

        // Confirm that we can call NextLocalSolarEclipse and it doesn't throw an exception.
        // Could add more stringent checking here later.
        const next = Astronomy.NextLocalSolarEclipse(eclipse.peak.time, observer);
        v(next.peak.time.tt);
    }

    if (lnum != expected_count) {
        console.error(`JS LocalSolarEclipse1: WRONG LINE COUNT = ${lnum}, expected ${expected_count}`);
        return 1;
    }

    if (skip_count > 6) {
        console.error(`JS LocalSolarEclipse1: EXCESSIVE SKIP COUNT = ${skip_count}`);
        return 1;
    }

    console.log(`JS LocalSolarEclipse1: PASS (${lnum} verified, ${skip_count} skipped, max minutes = ${max_minutes})`);
    return 0;
}


function TrimLine(line) {
    // Treat '#' as a comment character.
    const poundIndex = line.indexOf('#');
    if (poundIndex >= 0) {
        line = line.substr(0, poundIndex);
    }
    line = line.trim();
    return line;
}


function ParseEvent(time_str, alt_str, required) {
    if (required) {
        return {
            time: new Astronomy.MakeTime(new Date(time_str)),
            altitude: float(alt_str)
        };
    }
    if (time_str !== '-') {
        throw `Expected event time to be "-" but found "${time_str}"`;
    }
    return null;
}


function LocalSolarEclipse2() {
    // Test ability to calculate local solar eclipse conditions away from
    // the peak position on the Earth.

    const filename = 'eclipse/local_solar_eclipse.txt';
    const lines = ReadLines(filename);
    let lnum = 0;
    let verify_count = 0;
    let max_minutes = 0.0;
    let max_degrees = 0.0;

    function CheckEvent(calc, expect) {
        const diff_minutes = (24 * 60) * abs(expect.time.ut - calc.time.ut);
        if (diff_minutes > max_minutes) {
            max_minutes = diff_minutes;
        }
        if (diff_minutes > 1.0) {
            throw `CheckEvent(${filename} line ${lnum}): EXCESSIVE TIME ERROR: ${diff_minutes} minutes.`;
        }

        // Ignore discrepancies for negative altitudes, because of quirky and irrelevant differences in refraction models.
        if (expect.altitude > 0.0) {
            const diff_alt = abs(expect.altitude - calc.altitude);
            if (diff_alt > max_degrees) {
                max_degrees = diff_alt;
            }
            if (diff_alt > 0.5) {
                throw `CheckEvent(${filename} line ${lnum}): EXCESSIVE ALTITUDE ERROR: ${diff_alt} degrees.`;
            }
        }
    }

    for (let line of lines) {
        ++lnum;
        line = TrimLine(line);
        if (line === '') continue;
        const token = line.split(/\s+/);
        if (token.length !== 13) {
            console.error(`JS LocalSolarEclipse2(${filename} line ${lnum}): Incorrect token count = ${token.length}`);
            return 1;
        }
        const latitude = float(token[0]);
        const longitude = float(token[1]);
        const observer = new Astronomy.Observer(latitude, longitude, 0);
        const typeChar = token[2];
        const expected_kind = {
            'P': Astronomy.EclipseKind.Partial,
            'A': Astronomy.EclipseKind.Annular,
            'T': Astronomy.EclipseKind.Total,
            'H': Astronomy.EclipseKind.Total
        }[typeChar];
        const p1    = ParseEvent(token[3],  token[4],   true);
        const t1    = ParseEvent(token[5],  token[6],   (typeChar !== 'P'));
        const peak  = ParseEvent(token[7],  token[8],   true);
        const t2    = ParseEvent(token[9],  token[10],  (typeChar !== 'P'));
        const p2    = ParseEvent(token[11], token[12],  true);

        const search_time = p1.time.AddDays(-20);
        const eclipse = Astronomy.SearchLocalSolarEclipse(search_time, observer);
        if (eclipse.kind !== expected_kind) {
            console.error(`JS LocalSolarEclipse2(${filename} line ${lnum}): expected eclipse kind "${expected_kind}" but found "${eclipse.kind}".`);
            return 1;
        }
        CheckEvent(eclipse.peak, peak);
        CheckEvent(eclipse.partial_begin, p1);
        CheckEvent(eclipse.partial_end, p2);
        if (typeChar != 'P') {
            CheckEvent(eclipse.total_begin, t1);
            CheckEvent(eclipse.total_end, t2);
        }
        ++verify_count;
    }

    console.log(`JS LocalSolarEclipse2: PASS (${verify_count} verified, max_minutes = ${max_minutes}, max_degrees = ${max_degrees})`);
    return 0;
}


function LocalSolarEclipse() {
    return (
        LocalSolarEclipse1() ||
        LocalSolarEclipse2()
    );
}


function GlobalAnnularCase(dateText, obscuration) {
    // Search for the first solar eclipse that occurs after the given date.
    const time = Astronomy.MakeTime(new Date(dateText));
    const eclipse = Astronomy.SearchGlobalSolarEclipse(time);

    // Verify the eclipse is within 1 day after the search basis time.
    const dt = v(eclipse.peak.ut - time.ut);
    if (dt < 0.0 || dt > 1.0) {
        console.error(`JS GlobalAnnularCase(${dateText}) FAIL: found eclipse ${dt.toFixed(4)} days after search time.`);
        return 1;
    }

    // Verify we found an annular solar eclipse.
    if (eclipse.kind !== Astronomy.EclipseKind.Annular) {
        console.error(`JS GlobalAnnularCase(${dateText}) FAIL: expected annular eclipse but found ${eclipse.kind}.`);
        return 1;
    }

    // Check how accurately we calculated obscuration.
    const diff = v(eclipse.obscuration - obscuration);
    if (abs(diff) > 0.0000904) {
        console.error(`JS GlobalAnnularCase(${dateText}) FAIL: excessive obscuration error = ${diff.toFixed(8)}, expected = ${obscuration.toFixed(8)}, calculated = ${eclipse.obscuration.toFixed(8)}.`);
        return 1;
    }
    Debug(`JS GlobalAnnularCase(${dateText}) obscuration error = ${diff.toFixed(8)}`);
    return 0;
}


function LocalSolarCase(dateText, latitude, longitude, kind, obscuration, tolerance) {
    const time = Astronomy.MakeTime(new Date(dateText));
    const observer = new Astronomy.Observer(latitude, longitude, 0.0);
    const eclipse = Astronomy.SearchLocalSolarEclipse(time, observer);
    const dt = v(eclipse.peak.time.ut - time.ut);
    if (dt < 0.0 || dt > 1.0) {
        console.error(`JS LocalSolarCase(${dateText}) FAIL: elapsed time = ${dt.toFixed(4)} days.`);
        return 1;
    }

    if (eclipse.kind !== kind) {
        console.error(`JS LocalSolarCase(${dateText}) FAIL: expected ${kind} but found ${eclipse.kind}.`);
        return 1;
    }

    const diff = v(eclipse.obscuration - obscuration);
    if (abs(diff) > tolerance) {
        console.error(`JS LocalSolarCase(${dateText}) FAIL: obscuration diff = ${diff}, expected = ${obscuration}, actual = ${eclipse.obscuration}.`);
        return 1;
    }
    Debug(`JS LocalSolarCase(${dateText}): obscuration diff = ${diff.toFixed(8)}`);
    return 0;
}


function SolarFraction() {
    return (
        // Verify global solar eclipse obscurations for annular eclipses only.
        // This is because they are the only nontrivial values for global solar eclipses.
        // The trivial values are all validated exactly by GlobalSolarEclipseTest().
        GlobalAnnularCase('2023-10-14', 0.90638) ||    // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2023Oct14Aprime.html
        GlobalAnnularCase('2024-10-02', 0.86975) ||    // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2024Oct02Aprime.html
        GlobalAnnularCase('2027-02-06', 0.86139) ||    // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2027Feb06Aprime.html
        GlobalAnnularCase('2028-01-26', 0.84787) ||    // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2028Jan26Aprime.html
        GlobalAnnularCase('2030-06-01', 0.89163) ||    // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2030Jun01Aprime.html

        // Verify obscuration values for specific locations on the Earth.
        // Local solar eclipse calculations include obscuration for all types of eclipse, not just annular and total.
        LocalSolarCase('2023-10-14',  11.3683,  -83.1017, Astronomy.EclipseKind.Annular, 0.90638, 0.000080) ||   // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2023Oct14Aprime.html
        LocalSolarCase('2023-10-14',  25.78,    -80.22,   Astronomy.EclipseKind.Partial, 0.578,   0.000023) ||   // https://aa.usno.navy.mil/calculated/eclipse/solar?eclipse=22023&lat=25.78&lon=-80.22&label=Miami%2C+FL&height=0&submit=Get+Data
        LocalSolarCase('2023-10-14',  30.2666,  -97.7000, Astronomy.EclipseKind.Partial, 0.8867,  0.001016) ||   // http://astro.ukho.gov.uk/eclipse/0332023/Austin_TX_United_States_2023Oct14.png
        LocalSolarCase('2024-04-08',  25.2900, -104.1383, Astronomy.EclipseKind.Total,   1.0,     0.0     ) ||   // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2024Apr08Tprime.html
        LocalSolarCase('2024-04-08',  37.76,   -122.44,   Astronomy.EclipseKind.Partial, 0.340,   0.000604) ||   // https://aa.usno.navy.mil/calculated/eclipse/solar?eclipse=12024&lat=37.76&lon=-122.44&label=San+Francisco%2C+CA&height=0&submit=Get+Data
        LocalSolarCase('2024-10-02', -21.9533, -114.5083, Astronomy.EclipseKind.Annular, 0.86975, 0.000061) ||   // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2024Oct02Aprime.html
        LocalSolarCase('2024-10-02', -33.468,   -70.636,  Astronomy.EclipseKind.Partial, 0.436,   0.000980) ||   // https://aa.usno.navy.mil/calculated/eclipse/solar?eclipse=22024&lat=-33.468&lon=-70.636&label=Santiago%2C+Chile&height=0&submit=Get+Data
        LocalSolarCase('2030-06-01',  56.525,    80.0617, Astronomy.EclipseKind.Annular, 0.89163, 0.000067) ||   // https://www.eclipsewise.com/solar/SEprime/2001-2100/SE2030Jun01Aprime.html
        LocalSolarCase('2030-06-01',  40.388,    49.914,  Astronomy.EclipseKind.Partial, 0.67240, 0.000599) ||   // http://xjubier.free.fr/en/site_pages/SolarEclipseCalc_Diagram.html
        LocalSolarCase('2030-06-01',  40.3667,   49.8333, Astronomy.EclipseKind.Partial, 0.6736,  0.001464) ||   // http://astro.ukho.gov.uk/eclipse/0132030/Baku_Azerbaijan_2030Jun01.png

        Pass('SolarFraction')
    );
}


function PlanetApsis() {
    const start_time = Astronomy.MakeTime(new Date('1700-01-01T00:00:00Z'));
    let pindex = 0;
    for (let body of ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto']) {
        let count = 1;
        const period = Astronomy.PlanetOrbitalPeriod(body);
        let min_interval = -1.0;
        let max_interval = -1.0;
        let max_diff_days = 0.0;
        let max_dist_ratio = 0.0;
        let apsis = Astronomy.SearchPlanetApsis(body, start_time);
        const filename = `apsides/apsis_${pindex}.txt`;
        const lines = ReadLines(filename);
        for (const line of lines) {
            if (line.trim() == '') {
                continue;
            }
            const token = line.split(/\s+/);
            if (token.length != 3) {
                throw `${filename} line ${count}: Invalid data format: ${token.length} tokens.`;
            }
            const expected_kind = int(token[0]);
            const expected_time = Astronomy.MakeTime(new Date(token[1]));
            const expected_distance = float(token[2]);
            if (expected_kind !== apsis.kind) {
                throw `${filename} line ${count}: WRONG APSIS KIND: expected ${expected_kind}, found ${apsis.kind}`;
            }
            const diff_days = abs(expected_time.tt - apsis.time.tt);
            max_diff_days = max(max_diff_days, diff_days);
            const diff_degrees = (diff_days / period) * 360;
            const degree_threshold = 0.1;
            if (diff_degrees > 0.1) {
                throw `APSIS FAIL: ${body} exceeded angular threshold (${diff_degrees} vs ${degree_threshold} degrees).`;
            }
            const diff_dist_ratio = abs(expected_distance - apsis.dist_au) / expected_distance;
            max_dist_ratio = max(max_dist_ratio, diff_dist_ratio);
            if (diff_dist_ratio > 1.05e-4) {
                throw `${filename} line ${count}: distance ratio ${diff_dist_ratio} is too large.`;
            }

            // Calculate the next apsis
            const prev_time = apsis.time;
            apsis = Astronomy.NextPlanetApsis(body, apsis);
            ++count;
            const interval = apsis.time.tt - prev_time.tt;
            if (min_interval < 0.0) {
                min_interval = max_interval = interval;
            } else {
                min_interval = min(min_interval, interval);
                max_interval = max(max_interval, interval);
            }
        }
        if (count < 2) {
            throw `Failed to find apsides for ${body}`;
        }
        Debug(`JS PlanetApsis: ${count} apsides for ${body} -- intervals: min=${min_interval}, max=${max_interval}, ratio=${max_interval/min_interval}; max day=${max_diff_days}, degrees=${(max_diff_days / period) * 360}, dist ratio=${max_dist_ratio}`);
        ++pindex;
    }
    return 0;
}


function Elongation() {
    function SearchElongTest() {
        // Max elongation data obtained from:
        // http://www.skycaramba.com/greatest_elongations.shtml
        TestMaxElong('Mercury', '2010-01-17T05:22Z', '2010-01-27T05:22Z', 24.80, 'morning');
        TestMaxElong('Mercury', '2010-05-16T02:15Z', '2010-05-26T02:15Z', 25.10, 'morning');
        TestMaxElong('Mercury', '2010-09-09T17:24Z', '2010-09-19T17:24Z', 17.90, 'morning');
        TestMaxElong('Mercury', '2010-12-30T14:33Z', '2011-01-09T14:33Z', 23.30, 'morning');
        TestMaxElong('Mercury', '2011-04-27T19:03Z', '2011-05-07T19:03Z', 26.60, 'morning');
        TestMaxElong('Mercury', '2011-08-24T05:52Z', '2011-09-03T05:52Z', 18.10, 'morning');
        TestMaxElong('Mercury', '2011-12-13T02:56Z', '2011-12-23T02:56Z', 21.80, 'morning');
        TestMaxElong('Mercury', '2012-04-08T17:22Z', '2012-04-18T17:22Z', 27.50, 'morning');
        TestMaxElong('Mercury', '2012-08-06T12:04Z', '2012-08-16T12:04Z', 18.70, 'morning');
        TestMaxElong('Mercury', '2012-11-24T22:55Z', '2012-12-04T22:55Z', 20.60, 'morning');
        TestMaxElong('Mercury', '2013-03-21T22:02Z', '2013-03-31T22:02Z', 27.80, 'morning');
        TestMaxElong('Mercury', '2013-07-20T08:51Z', '2013-07-30T08:51Z', 19.60, 'morning');
        TestMaxElong('Mercury', '2013-11-08T02:28Z', '2013-11-18T02:28Z', 19.50, 'morning');
        TestMaxElong('Mercury', '2014-03-04T06:38Z', '2014-03-14T06:38Z', 27.60, 'morning');
        TestMaxElong('Mercury', '2014-07-02T18:22Z', '2014-07-12T18:22Z', 20.90, 'morning');
        TestMaxElong('Mercury', '2014-10-22T12:36Z', '2014-11-01T12:36Z', 18.70, 'morning');
        TestMaxElong('Mercury', '2015-02-14T16:20Z', '2015-02-24T16:20Z', 26.70, 'morning');
        TestMaxElong('Mercury', '2015-06-14T17:10Z', '2015-06-24T17:10Z', 22.50, 'morning');
        TestMaxElong('Mercury', '2015-10-06T03:20Z', '2015-10-16T03:20Z', 18.10, 'morning');
        TestMaxElong('Mercury', '2016-01-28T01:22Z', '2016-02-07T01:22Z', 25.60, 'morning');
        TestMaxElong('Mercury', '2016-05-26T08:45Z', '2016-06-05T08:45Z', 24.20, 'morning');
        TestMaxElong('Mercury', '2016-09-18T19:27Z', '2016-09-28T19:27Z', 17.90, 'morning');
        TestMaxElong('Mercury', '2017-01-09T09:42Z', '2017-01-19T09:42Z', 24.10, 'morning');
        TestMaxElong('Mercury', '2017-05-07T23:19Z', '2017-05-17T23:19Z', 25.80, 'morning');
        TestMaxElong('Mercury', '2017-09-02T10:14Z', '2017-09-12T10:14Z', 17.90, 'morning');
        TestMaxElong('Mercury', '2017-12-22T19:48Z', '2018-01-01T19:48Z', 22.70, 'morning');
        TestMaxElong('Mercury', '2018-04-19T18:17Z', '2018-04-29T18:17Z', 27.00, 'morning');
        TestMaxElong('Mercury', '2018-08-16T20:35Z', '2018-08-26T20:35Z', 18.30, 'morning');
        TestMaxElong('Mercury', '2018-12-05T11:34Z', '2018-12-15T11:34Z', 21.30, 'morning');
        TestMaxElong('Mercury', '2019-04-01T19:40Z', '2019-04-11T19:40Z', 27.70, 'morning');
        TestMaxElong('Mercury', '2019-07-30T23:08Z', '2019-08-09T23:08Z', 19.00, 'morning');
        TestMaxElong('Mercury', '2019-11-18T10:31Z', '2019-11-28T10:31Z', 20.10, 'morning');
        TestMaxElong('Mercury', '2010-03-29T23:32Z', '2010-04-08T23:32Z', 19.40, 'evening');
        TestMaxElong('Mercury', '2010-07-28T01:03Z', '2010-08-07T01:03Z', 27.40, 'evening');
        TestMaxElong('Mercury', '2010-11-21T15:42Z', '2010-12-01T15:42Z', 21.50, 'evening');
        TestMaxElong('Mercury', '2011-03-13T01:07Z', '2011-03-23T01:07Z', 18.60, 'evening');
        TestMaxElong('Mercury', '2011-07-10T04:56Z', '2011-07-20T04:56Z', 26.80, 'evening');
        TestMaxElong('Mercury', '2011-11-04T08:40Z', '2011-11-14T08:40Z', 22.70, 'evening');
        TestMaxElong('Mercury', '2012-02-24T09:39Z', '2012-03-05T09:39Z', 18.20, 'evening');
        TestMaxElong('Mercury', '2012-06-21T02:00Z', '2012-07-01T02:00Z', 25.70, 'evening');
        TestMaxElong('Mercury', '2012-10-16T21:59Z', '2012-10-26T21:59Z', 24.10, 'evening');
        TestMaxElong('Mercury', '2013-02-06T21:24Z', '2013-02-16T21:24Z', 18.10, 'evening');
        TestMaxElong('Mercury', '2013-06-02T16:45Z', '2013-06-12T16:45Z', 24.30, 'evening');
        TestMaxElong('Mercury', '2013-09-29T09:59Z', '2013-10-09T09:59Z', 25.30, 'evening');
        TestMaxElong('Mercury', '2014-01-21T10:00Z', '2014-01-31T10:00Z', 18.40, 'evening');
        TestMaxElong('Mercury', '2014-05-15T07:06Z', '2014-05-25T07:06Z', 22.70, 'evening');
        TestMaxElong('Mercury', '2014-09-11T22:20Z', '2014-09-21T22:20Z', 26.40, 'evening');
        TestMaxElong('Mercury', '2015-01-04T20:26Z', '2015-01-14T20:26Z', 18.90, 'evening');
        TestMaxElong('Mercury', '2015-04-27T04:46Z', '2015-05-07T04:46Z', 21.20, 'evening');
        TestMaxElong('Mercury', '2015-08-25T10:20Z', '2015-09-04T10:20Z', 27.10, 'evening');
        TestMaxElong('Mercury', '2015-12-19T03:11Z', '2015-12-29T03:11Z', 19.70, 'evening');
        TestMaxElong('Mercury', '2016-04-08T14:00Z', '2016-04-18T14:00Z', 19.90, 'evening');
        TestMaxElong('Mercury', '2016-08-06T21:24Z', '2016-08-16T21:24Z', 27.40, 'evening');
        TestMaxElong('Mercury', '2016-12-01T04:36Z', '2016-12-11T04:36Z', 20.80, 'evening');
        TestMaxElong('Mercury', '2017-03-22T10:24Z', '2017-04-01T10:24Z', 19.00, 'evening');
        TestMaxElong('Mercury', '2017-07-20T04:34Z', '2017-07-30T04:34Z', 27.20, 'evening');
        TestMaxElong('Mercury', '2017-11-14T00:32Z', '2017-11-24T00:32Z', 22.00, 'evening');
        TestMaxElong('Mercury', '2018-03-05T15:07Z', '2018-03-15T15:07Z', 18.40, 'evening');
        TestMaxElong('Mercury', '2018-07-02T05:24Z', '2018-07-12T05:24Z', 26.40, 'evening');
        TestMaxElong('Mercury', '2018-10-27T15:25Z', '2018-11-06T15:25Z', 23.30, 'evening');
        TestMaxElong('Mercury', '2019-02-17T01:23Z', '2019-02-27T01:23Z', 18.10, 'evening');
        TestMaxElong('Mercury', '2019-06-13T23:14Z', '2019-06-23T23:14Z', 25.20, 'evening');
        TestMaxElong('Mercury', '2019-10-10T04:00Z', '2019-10-20T04:00Z', 24.60, 'evening');
        TestMaxElong('Venus', '2010-12-29T15:57Z', '2011-01-08T15:57Z', 47.00, 'morning');
        TestMaxElong('Venus', '2012-08-05T08:59Z', '2012-08-15T08:59Z', 45.80, 'morning');
        TestMaxElong('Venus', '2014-03-12T19:25Z', '2014-03-22T19:25Z', 46.60, 'morning');
        TestMaxElong('Venus', '2015-10-16T06:57Z', '2015-10-26T06:57Z', 46.40, 'morning');
        TestMaxElong('Venus', '2017-05-24T13:09Z', '2017-06-03T13:09Z', 45.90, 'morning');
        TestMaxElong('Venus', '2018-12-27T04:24Z', '2019-01-06T04:24Z', 47.00, 'morning');
        TestMaxElong('Venus', '2010-08-10T03:19Z', '2010-08-20T03:19Z', 46.00, 'evening');
        TestMaxElong('Venus', '2012-03-17T08:03Z', '2012-03-27T08:03Z', 46.00, 'evening');
        TestMaxElong('Venus', '2013-10-22T08:00Z', '2013-11-01T08:00Z', 47.10, 'evening');
        TestMaxElong('Venus', '2015-05-27T18:46Z', '2015-06-06T18:46Z', 45.40, 'evening');
        TestMaxElong('Venus', '2017-01-02T13:19Z', '2017-01-12T13:19Z', 47.10, 'evening');
        TestMaxElong('Venus', '2018-08-07T17:02Z', '2018-08-17T17:02Z', 45.90, 'evening');
    }

    function LoadData(filename) {
        const lines = ReadLines(filename);
        let data = [];
        for (let row of lines) {
            let token = row.split(/\s+/);
            data.push({date:new Date(token[0]), body:token[1]});
        }
        return data;
    }

    function TestFile(filename, startSearchYear, targetRelLon) {
        const data = LoadData(filename);
        const startDate = new Date(Date.UTC(startSearchYear, 0, 1));
        for (let item of data) {
            let time = Astronomy.SearchRelativeLongitude(item.body, targetRelLon, startDate);
            let diff_minutes = (time.date - item.date) / 60000;
            Debug(`JS ${item.body}: error = ${diff_minutes.toFixed(3)} minutes`);
            if (abs(diff_minutes) > 6.8)
                throw `!!! Excessive error for body ${item.body}`;
        }
    }

    function TestPlanet(outFileName, body, startYear, stopYear, zeroLonEventName) {
        let rlon = 0;
        let date = new Date(Date.UTC(startYear, 0, 1));
        let stopDate = new Date(Date.UTC(stopYear, 0, 1));
        let text = '';
        let count = 0;
        let prev_time, min_diff, max_diff, sum_diff=0;

        while (date < stopDate) {
            let event = (rlon === 0) ? zeroLonEventName : 'sup';
            let evt_time = Astronomy.SearchRelativeLongitude(body, rlon, date);
            if (prev_time) {
                // Check for consistent intervals.
                // Mainly I don't want to accidentally skip over an event!
                let day_diff = evt_time.tt - prev_time.tt;
                if (min_diff === undefined) {
                    min_diff = max_diff = day_diff;
                } else {
                    min_diff = min(min_diff, day_diff);
                    max_diff = max(max_diff, day_diff);
                }
                sum_diff += day_diff;
            }
            let geo = Astronomy.GeoVector(body, evt_time, false);
            let dist = sqrt(geo.x*geo.x + geo.y*geo.y + geo.z*geo.z);
            text += `e ${body} ${event} ${evt_time.tt} ${dist}\n`;
            rlon = 180 - rlon;
            date = evt_time.date;
            ++count;
            prev_time = evt_time;
        }

        fs.writeFileSync(outFileName, text);

        const ratio = max_diff / min_diff;
        Debug(`JS TestPlanet(${body}): ${count} events, interval min=${min_diff.toFixed(1)}, max=${max_diff.toFixed(1)}, avg=${(sum_diff/count).toFixed(1)}, ratio=${ratio.toFixed(3)}`);

        let thresh = {Mercury:1.65, Mars:1.30}[body] || 1.07;
        if (ratio > thresh)
            throw `TestPlanet: Excessive event interval ratio for ${body} = ${ratio}`;
    }

    function TestMaxElong(body, startText, verifyText, verifyAngle, verifyVisibility) {
        let startDate = new Date(startText);
        let verifyDate = new Date(verifyText);
        let evt = Astronomy.SearchMaxElongation(body, startDate);

        let hour_diff = (verifyDate - evt.time.date) / (1000 * 3600);
        let arcmin_diff = 60.0 * abs(evt.elongation - verifyAngle);
        Debug(`JS TestMaxElong: ${body.padStart(8)} ${evt.visibility.padStart(8)} elong=${evt.elongation.toFixed(2).padStart(5)} (${arcmin_diff.toFixed(2).padStart(4)} arcmin)  ${evt.time.toString()} (err ${hour_diff.toFixed(2).padStart(5)} hours)`);

        if (evt.visibility !== verifyVisibility)
            throw `TestMaxElong: expected visibility ${verifyVisibility}, but found ${evt.visibility}`;

        if (arcmin_diff > 3.4)
            throw `TestMaxElong: excessive angular error = ${angle_diff} arcmin`;

        if (abs(hour_diff) > 0.6)
            throw `TestMaxElong: excessive hour error = ${hour_diff}`;
    }

    TestFile('longitude/opposition_2018.txt', 2018, 0);

    for (let body of ['Mercury', 'Venus'])
        TestPlanet(`temp/js_longitude_${body}.txt`, body, 1700, 2200, 'inf');

    for (let body of ['Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto'])
        TestPlanet(`temp/js_longitude_${body}.txt`, body, 1700, 2200, 'opp');

    SearchElongTest();

    console.log('JS Elongation: PASS');
    return 0;
}


function Seasons() {
    function LoadTestData(filename) {
        // Moon  150 -45 2050-03-07T19:13Z s
        const lines = ReadLines(filename);
        let data = [];
        let lnum = 0;
        let minByMonth = [];
        let maxByMonth = [];
        for (let row of lines) {
            let token = row.split(/\s+/g);
            let item = {
                lnum: ++lnum,
                date: new Date(token[0]),
                name: token[1]  // Perihelion, Equinox, Solstice, Aphelion
            };
            data.push(item);

            if (item.name === 'Equinox' || item.name == 'Solstice') {
                let month = 1 + item.date.getUTCMonth();
                let format = item.date.toISOString().substring(8);
                if (!minByMonth[month] || format < minByMonth[month])
                    minByMonth[month] = format;
                if (!maxByMonth[month] || format > maxByMonth[month])
                    maxByMonth[month] = format;
            }
        }

        if (Verbose) {
            console.log(`JS Seasons LoadTestData: count = ${data.length}`);
            for (let month of [3, 6, 9, 12]) {
                console.log(`Month ${month}: earliest ${minByMonth[month]}, latest ${maxByMonth[month]}`);
            }
        }
        return data;
    }
    const data = LoadTestData('seasons/seasons.txt');
    let current_year;
    let seasons;
    let calc_date;
    let min_diff, max_diff, sum_diff = 0, count = 0;
    let month_max_diff = [];
    for (let item of data) {
        let year = item.date.getUTCFullYear();
        if (current_year !== year) {
            current_year = year;
            seasons = Astronomy.Seasons(year);
        }
        calc_date = null;
        let month = 1 + item.date.getUTCMonth();
        switch (item.name) {
        case 'Equinox':
            switch (month) {
            case 3:
                calc_date = seasons.mar_equinox.date;
                break;
            case 9:
                calc_date = seasons.sep_equinox.date;
                break;
            default:
                throw `ERROR: Invalid equinox date in test data: ${item.date.toISOString()}`;
            }
            break;
        case 'Solstice':
            switch (month) {
            case 6:
                calc_date = seasons.jun_solstice.date;
                break;
            case 12:
                calc_date = seasons.dec_solstice.date;
                break;
            default:
                throw `ERROR: Invalid solstice date in test data: ${item.date.toISOString()}`;
            }
            break;
        default:
            continue;   // ignore the other kinds of events for now
        }
        if (!calc_date)
            throw `ERROR: Missing calc_date for test date ${item.date.toISOString()}`;
        let diff_minutes = (calc_date - item.date) / 60000;
        if (abs(diff_minutes) > 2.37) {
            throw `ERROR: Excessive error in season calculation: ${diff_minutes.toFixed(3)} minutes`;
        }

        if (min_diff === undefined) {
            min_diff = max_diff = diff_minutes;
        } else {
            min_diff = min(min_diff, diff_minutes);
            max_diff = max(max_diff, diff_minutes);
        }
        if (month_max_diff[month] === undefined) {
            month_max_diff[month] = abs(diff_minutes);
        } else {
            month_max_diff[month] = max(month_max_diff[month], abs(diff_minutes));
        }
        sum_diff += diff_minutes;
        ++count;
    }

    if (Verbose) {
        console.log(`JS Seasons: n=${count}, minute errors: min=${min_diff.toFixed(3)}, avg=${(sum_diff/count).toFixed(3)}, max=${max_diff.toFixed(3)}`);
        for (let month of [3, 6, 9, 12]) {
            console.log(`JS Seasons: max diff by month ${month} = ${month_max_diff[month].toFixed(3)}`);
        }
    }
    console.log('JS Seasons: PASS');
    return 0;
}


function SeasonsIssue187() {
    // This is a regression test for:
    // https://github.com/cosinekitty/astronomy/issues/187
    // For years far from the present, the seasons search was sometimes failing.

    for (let year = -2000; year <= +9999; ++year) {
        try {
            Astronomy.Seasons(year);
        } catch (e) {
            console.error(`JS SeasonsIssue187: FAIL (year = ${year}): `, e);
            return 1;
        }
    }

    console.log('JS SeasonsIssue187: PASS');
    return 0;
}


function RiseSet() {
    function LoadTestData(filename) {
        // Moon  150 -45 2050-03-07T19:13Z s
        const lines = ReadLines(filename);
        let data = [];
        let lnum = 0;
        for (let row of lines) {
            let token = row.split(/\s+/g);
            data.push({
                lnum: ++lnum,
                body: token[0],
                lon: float(token[1]),
                lat: float(token[2]),
                date: new Date(token[3]),
                direction: { r:+1, s:-1 }[token[4]] || Fail(`Invalid event code ${token[4]}`)
            });
        }
        return data;
    }
    const data = LoadTestData('riseset/riseset.txt');

    // The test data is sorted by body, then geographic location, then date/time.

    const before_date = new Date();

    let body;
    let observer;
    let r_search_date, r_date;
    let s_search_date, s_date;
    let a_date, b_date, a_dir, b_dir;
    let sum_minutes = 0;
    let max_minutes = 0;
    for (let evt of data) {
        if (!observer || observer.latitude !== evt.lat || observer.longitude !== evt.lon || body !== evt.body) {
            // Every time we see a new geographic location, start a new iteration
            // of finding all rise/set times for that UTC calendar year.
            body = evt.body;
            observer = new Astronomy.Observer(evt.lat, evt.lon, 0);
            r_search_date = s_search_date = new Date(Date.UTC(evt.date.getUTCFullYear(), 0, 1));
            b_date = null;
            Debug(`JS RiseSet: ${body} ${evt.lat} ${evt.lon}`);
        }

        if (b_date) {
            // recycle the second event from the previous iteration as the first event
            a_date = b_date;
            a_dir = b_dir;
            b_date = null;
        } else {
            r_date = Astronomy.SearchRiseSet(body, observer, +1, r_search_date, 366) ||
                Fail(`JS RiseSet: Did not find ${body} rise after ${r_search_date.toISOString()}`);

            s_date = Astronomy.SearchRiseSet(body, observer, -1, s_search_date, 366) ||
                Fail(`JS RiseSet: Did not find ${body} set after ${s_search_date.toISOString()}`);

            // Expect the current event to match the earlier of the found dates.
            if (r_date.tt < s_date.tt) {
                a_date = r_date;
                b_date = s_date;
                a_dir = +1;
                b_dir = -1;
            } else {
                a_date = s_date;
                b_date = r_date;
                a_dir = -1;
                b_dir = +1;
            }

            r_search_date = r_date.AddDays(0.01).date;
            s_search_date = s_date.AddDays(0.01).date;
        }

        if (a_dir !== evt.direction) {
            Fail(`[line ${evt.lnum}] Expected ${body} dir=${evt.direction} at ${evt.date.toISOString()} but found ${a_dir} ${a_date.toString()}`);
        }

        let error_minutes = abs(a_date.date - evt.date) / 60000;
        sum_minutes += error_minutes * error_minutes;
        if (error_minutes > max_minutes) {
            max_minutes = error_minutes;
            Debug(`Line ${evt.lnum} : error = ${error_minutes.toFixed(4)}`);
        }
        if (error_minutes > 1.18) {
            console.log(`Expected ${evt.date.toISOString()}`);
            console.log(`Found    ${a_date.toString()}`);
            Fail("Excessive prediction time error.");
        }
    }

    const after_date = new Date();
    const elapsed_seconds = (after_date - before_date) / 1000;

    console.log(`JS RiseSet PASS: elapsed=${elapsed_seconds.toFixed(3)}, error in minutes: rms=${sqrt(sum_minutes/data.length).toFixed(4)}, max=${max_minutes.toFixed(4)}`);
    return 0;
}


function RiseSetSlot(ut1, ut2, dir, observer) {
    const nslots = 100;
    let maxDiff = 0.0;
    for (let i = 1; i < nslots; ++i) {
        const ut = ut1 + (i/nslots)*(ut2 - ut1);
        const time = Astronomy.MakeTime(ut);
        let result = Astronomy.SearchRiseSet(Astronomy.Body.Sun, observer, dir, time, -1.0);
        if (!result) {
            console.error(`JS RiseSetSlot(${dir}): backward slot search failed for sunrise before ${time}`);
            return 1;
        }
        let diff = SECONDS_PER_DAY * abs(result.ut - ut1);
        if (diff > maxDiff) maxDiff = diff;

        result = Astronomy.SearchRiseSet(Astronomy.Body.Sun, observer, dir, time, +1.0);
        if (!result) {
            console.error(`JS RiseSetSlot(${dir}): forward slot search failed for sunrise before ${time}`);
            return 1;
        }
        diff = SECONDS_PER_DAY * abs(result.ut - ut2);
        if (diff > maxDiff) maxDiff = diff;
    }

    if (maxDiff > 0.13) {
        console.error(`JS RiseSetSlot(${dir}): EXCESSIVE slot-test discrepancy = ${maxDiff.toFixed(6)} seconds.`);
        return 1;
    }
    Debug(`JS RiseSetSlot(${dir}): slot-test discrepancy = ${maxDiff.toFixed(6)} seconds.`);
    return 0;
}


function RiseSetReverse() {
    // Verify that the rise/set search works equally well forwards and backwards in time.
    const nsamples = 5000;
    const nudge = 0.1;
    const utList = [];
    const observer = new Astronomy.Observer(30.5, -90.7, 0.0);
    let dtMin = +1000;
    let dtMax = -1000;
    let maxDiff = 0;

    // Find alternating sunrise/sunset events in forward chronological order.
    let dir = +1;
    let time = Astronomy.MakeTime(new Date('2022-01-01T00:00:00Z'));
    let i;
    for (i = 0; i < nsamples; ++i) {
        const result = Astronomy.SearchRiseSet(Astronomy.Body.Sun, observer, dir, time, +1.0);
        if (!result) {
            console.error(`JS RiseSetReverse: cannot find dir=${dir} event after ${time}.`);
            return 1;
        }
        utList.push(result.ut);
        if (i > 0) {
            const dt = v(utList[i] - utList[i-1]);
            if (dt < dtMin) dtMin = dt;
            if (dt > dtMax) dtMax = dt;
        }
        dir *= -1;
        time = result.AddDays(+nudge);
    }

    Debug(`JS RiseSetReverse: dtMin=${dtMin.toFixed(6)} days, dtMax=${dtMax.toFixed(6)} days.`);
    if (dtMin < 0.411 || dtMax > 0.589) {
        console.error(`JS RiseSetReverse: Invalid intervals between sunrise/sunset.`);
        return 1;
    }

    // Perform the same search in reverse. Verify we get consistent rise/set times.
    for (i = nsamples-1; i >= 0; --i) {
        dir *= -1;
        const result = Astronomy.SearchRiseSet(Astronomy.Body.Sun, observer, dir, time, -1.0);
        if (!result) {
            console.error(`JS RiseSetReverse: cannot find dir=${dir} event before ${time}.`);
            return 1;
        }
        const diff = SECONDS_PER_DAY * abs(utList[i] - result.ut);
        if (diff > maxDiff) maxDiff = diff;
        time = result.AddDays(-nudge);
    }

    if (maxDiff > 0.1) {
        console.error(`C# RiseSetReverse: EXCESSIVE forward/backward discrepancy = ${maxDiff.toFixed(6)} seconds.`);
        return 1;
    }
    Debug(`JS RiseSetReverse: forward/backward discrepancy = ${maxDiff.toFixed(6)} seconds.`);

    // All even indexes in utList hold sunrise times.
    // All odd indexes in utList hold sunset times.
    // Verify that forward/backward searches for consecutive sunrises/sunsets
    // resolve correctly for 100 time slots between them.
    const k = (nsamples / 2) & ~1;
    if (RiseSetSlot(utList[k+0], utList[k+2], +1, observer)) return 1;
    if (RiseSetSlot(utList[k+1], utList[k+3], -1, observer)) return 1;
    console.log('JS RiseSetReverse: PASS');
    return 0;
}

function VectorDiff(a, b) {
    const dx = a.x - b.x;
    const dy = a.y - b.y;
    const dz = a.z - b.z;
    return sqrt(dx*dx + dy*dy + dz*dz);
}

function Rotation() {
    function CompareMatrices(caller, a, b, tolerance) {
        for (let i=0; i<3; ++i) {
            for (let j=0; j<3; ++j) {
                const diff = abs(a.rot[i][j] - b.rot[i][j]);
                if (diff > tolerance) {
                    throw `ERROR(${caller}): matrix[${i}][${j}] = ${a.rot[i][j]}, expected ${b.rot[i][j]}, diff ${diff}`;
                }
            }
        }
    }

    function CompareVectors(caller, a, b, tolerance) {
        let diff;

        diff = abs(a.x - b.x);
        if (diff > tolerance) {
            throw `ERROR(${caller}): vector x = ${a.x}, expected ${b.x}, diff ${diff}`;
        }

        diff = abs(a.y - b.y);
        if (diff > tolerance) {
            throw `ERROR(${caller}): vector y = ${a.y}, expected ${b.y}, diff ${diff}`;
        }

        diff = abs(a.z - b.z);
        if (diff > tolerance) {
            throw `ERROR(${caller}): vector z = ${a.z}, expected ${b.z}, diff ${diff}`;
        }
    }

    function Rotation_MatrixInverse() {
        const a = Astronomy.MakeRotation([
            [1, 4, 7],
            [2, 5, 8],
            [3, 6, 9]
        ]);

        const v = Astronomy.MakeRotation([
            [1, 2, 3],
            [4, 5, 6],
            [7, 8, 9]
        ]);

        const b = Astronomy.InverseRotation(a);
        CompareMatrices('Rotation_MatrixInverse', b, v, 0);
    }

    function Rotation_MatrixMultiply() {
        const a = Astronomy.MakeRotation([
            [1, 4, 7],
            [2, 5, 8],
            [3, 6, 9]
        ]);

        const b = Astronomy.MakeRotation([
            [10, 13, 16],
            [11, 14, 17],
            [12, 15, 18]
        ]);

        const v = Astronomy.MakeRotation([
            [84, 201, 318],
            [90, 216, 342],
            [96, 231, 366]
        ]);

        const c = Astronomy.CombineRotation(b, a);
        CompareMatrices('Rotation_MatrixMultiply', c, v, 0);
    }

    function Test_GAL_EQJ_NOVAS(filename) {
        const THRESHOLD_SECONDS = 8.8;
        const rot = Astronomy.Rotation_EQJ_GAL();
        const lines = ReadLines(filename);
        const time = new Astronomy.AstroTime(0);    // placeholder time - value does not matter
        let lnum = 0;
        let max_diff = 0;
        for (let row of lines) {
            ++lnum;
            let token = row.trim().split(/\s+/);
            if (token.length !== 4)
                throw `Test_GAL_EQJ_NOVAS(${filename} line ${lnum}): expected 4 tokens, found ${token.length}`;
            const ra = float(token[0]);
            const dec = float(token[1]);
            const glon = float(token[2]);
            const glat = float(token[3]);

            const eqj_sphere = new Astronomy.Spherical(dec, 15*ra, 1);
            const eqj_vec = Astronomy.VectorFromSphere(eqj_sphere, time);
            const gal_vec = Astronomy.RotateVector(rot, eqj_vec);
            const gal_sphere = Astronomy.SphereFromVector(gal_vec);
            const dlat = v(gal_sphere.lat - glat);
            const dlon = cos(Astronomy.DEG2RAD * glat) * v(gal_sphere.lon - glon);
            const diff = 3600 * sqrt(dlon*dlon + dlat*dlat);
            if (diff > THRESHOLD_SECONDS)
                throw `Test_GAL_EQJ_NOVAS(${filename} line ${lnum}): EXCESSIVE ERROR = ${diff.toFixed(3)} arcseconds.`;
            if (diff > max_diff)
                max_diff = diff;
        }
        Debug(`JS Test_GAL_EQJ_NOVAS: PASS. max_diff = ${max_diff.toFixed(3)} arcseconds.`);
        return 0;
    }

    function Test_EQJ_EQD(body) {
        /* Verify conversion of equatorial J2000 to equatorial of-date, and back. */
        /* Use established functions to calculate spherical coordinates for the body, in both EQJ and EQD. */
        const time = Astronomy.MakeTime(new Date('2019-12-08T20:50:00Z'));
        const observer = new Astronomy.Observer(+35, -85, 0);
        const eq2000 = Astronomy.Equator(body, time, observer, false, true);
        const eqdate = Astronomy.Equator(body, time, observer, true, true);

        /* Convert EQJ spherical coordinates to vector. */
        const v2000 = eq2000.vec;

        /* Find rotation matrix. */
        const r = Astronomy.Rotation_EQJ_EQD(time);

        /* Rotate EQJ vector to EQD vector. */
        const vdate = Astronomy.RotateVector(r, v2000);

        /* Convert vector back to angular equatorial coordinates. */
        let equcheck = Astronomy.EquatorFromVector(vdate);

        /* Compare the result with the eqdate. */
        const ra_diff = abs(equcheck.ra - eqdate.ra);
        const dec_diff = abs(equcheck.dec - eqdate.dec);
        const dist_diff = abs(equcheck.dist - eqdate.dist);
        Debug(`JS Test_EQJ_EQD: ${body} ra=${eqdate.ra}, dec=${eqdate.dec}, dist=${eqdate.dist}, ra_diff=${ra_diff}, dec_diff=${dec_diff}, dist_diff=${dist_diff}`);
        if (ra_diff > 1.0e-14 || dec_diff > 1.0e-14 || dist_diff > 4.0e-15)
            throw 'Test_EQJ_EQD: EXCESSIVE ERROR';

        /* Perform the inverse conversion back to equatorial J2000 coordinates. */
        const ir = Astronomy.Rotation_EQD_EQJ(time);
        const t2000 = Astronomy.RotateVector(ir, vdate);
        const diff = VectorDiff(t2000, v2000);
        Debug(`JS Test_EQJ_EQD: ${body} inverse diff = ${diff}`);
        if (diff > 5.0e-15)
            throw 'Test_EQJ_EQD: EXCESSIVE INVERSE ERROR';
    }

    function Test_EQD_HOR(body) {
        /* Use existing functions to calculate horizontal coordinates of the body for the time+observer. */
        const time = Astronomy.MakeTime(new Date('1970-12-13T05:15:00Z'));
        const observer = new Astronomy.Observer(-37, +45, 0);
        const eqd = Astronomy.Equator(body, time, observer, true, true);
        Debug(`JS Test_EQD_HOR ${body}: OFDATE ra=${eqd.ra}, dec=${eqd.dec}`);
        const hor = Astronomy.Horizon(time, observer, eqd.ra, eqd.dec, 'normal');

        /* Calculate the position of the body as an equatorial vector of date. */
        const vec_eqd = eqd.vec;

        /* Calculate rotation matrix to convert equatorial J2000 vector to horizontal vector. */
        const rot = Astronomy.Rotation_EQD_HOR(time, observer);

        /* Rotate the equator of date vector to a horizontal vector. */
        const vec_hor = Astronomy.RotateVector(rot, vec_eqd);

        /* Convert the horizontal vector to horizontal angular coordinates. */
        const xsphere = Astronomy.HorizonFromVector(vec_hor, 'normal');
        const diff_alt = abs(xsphere.lat - hor.altitude);
        const diff_az = abs(xsphere.lon - hor.azimuth);

        Debug(`JS Test_EQD_HOR ${body}: trusted alt=${hor.altitude}, az=${hor.azimuth}; test alt=${xsphere.lat}, az=${xsphere.lon}; diff_alt=${diff_alt}, diff_az=${diff_az}`);
        if (diff_alt > 4.3e-14 || diff_az > 1.2e-13)
            throw 'Test_EQD_HOR: EXCESSIVE HORIZONTAL ERROR.';

        /* Confirm that we can convert back to horizontal vector. */
        const check_hor = Astronomy.VectorFromHorizon(xsphere, time, 'normal');
        let diff = VectorDiff(check_hor, vec_hor);
        Debug(`JS Test_EQD_HOR ${body}: horizontal recovery: diff = ${diff}`);
        if (diff > 2.0e-15)
            throw 'Test_EQD_HOR: EXCESSIVE ERROR IN HORIZONTAL RECOVERY.';

        /* Verify the inverse translation from horizontal vector to equatorial of-date vector. */
        const irot = Astronomy.Rotation_HOR_EQD(time, observer);
        const check_eqd = Astronomy.RotateVector(irot, vec_hor);
        diff = VectorDiff(check_eqd, vec_eqd);
        Debug(`JS Test_EQD_HOR ${body}: OFDATE inverse rotation diff = ${diff}`);
        if (diff > 2.3e-15)
            throw 'Test_EQD_HOR: EXCESSIVE OFDATE INVERSE HORIZONTAL ERROR.';

        /* Exercise HOR to EQJ translation. */
        const eqj = Astronomy.Equator(body, time, observer, false, true);
        const vec_eqj = eqj.vec;
        const yrot = Astronomy.Rotation_HOR_EQJ(time, observer);
        const check_eqj = Astronomy.RotateVector(yrot, vec_hor);
        diff = VectorDiff(check_eqj, vec_eqj);
        Debug(`JS Test_EQD_HOR ${body}: J2000 inverse rotation diff = ${diff}`);
        if (diff > 6.0e-15)
            throw 'Test_EQD_HOR: EXCESSIVE J2000 INVERSE HORIZONTAL ERROR.';

        /* Verify the inverse translation: EQJ to HOR. */
        const zrot = Astronomy.Rotation_EQJ_HOR(time, observer);
        const another_hor = Astronomy.RotateVector(zrot, vec_eqj);
        diff = VectorDiff(another_hor, vec_hor);
        Debug(`JS Test_EQD_HOR ${body}: EQJ inverse rotation diff = ${diff}`);
        if (diff > 3.0e-15)
            throw 'Test_EQD_HOR: EXCESSIVE EQJ INVERSE HORIZONTAL ERROR.';
    }

    function Test_EQD_ECT() {
        let time = Astronomy.MakeTime(new Date('1900-01-01T00:00:00Z'));
        const stopTime = Astronomy.MakeTime(new Date('2100-01-01T00:00:00Z'));
        let max_diff = 0.0;
        let count = 0;
        while (time.ut <= stopTime.ut) {
            // Get Moon's geocentric position in EQJ.
            const eqj = Astronomy.GeoMoon(time);

            // Convert EQJ to EQD.
            const eqj_eqd = Astronomy.Rotation_EQJ_EQD(time);
            const eqd = Astronomy.RotateVector(eqj_eqd, eqj);

            // Convert EQD to ECT.
            const eqd_ect = Astronomy.Rotation_EQD_ECT(time);
            const ect = Astronomy.RotateVector(eqd_ect, eqd);

            // Independently get the Moon's spherical coordinates in ECT.
            const sphere = Astronomy.EclipticGeoMoon(time);

            // Convert spherical coordinates to ECT vector.
            const check_ect = Astronomy.VectorFromSphere(sphere, time);

            // Verify the two ECT vectors are identical, within tolerance.
            max_diff = Math.max(max_diff, VectorDiff(ect, check_ect));

            time = time.AddDays(10.0);
            ++count;
        }
        if (max_diff > 3.743e-18)
            throw `Test_EQD_ECT: excessive vector diff = ${max_diff.toExponential(6)} AU.`;

        Debug(`JS Test_EQD_ECT: PASS: count = ${count}, max_diff = ${max_diff.toExponential(6)} AU.`);
    }

    const IdentityMatrix = Astronomy.IdentityMatrix();

    function CheckInverse(aname, bname, arot, brot) {
        const crot = Astronomy.CombineRotation(arot, brot);
        CompareMatrices(`CheckInverse(${aname},${bname})`, crot, IdentityMatrix, 2.0e-15);
    }

    function CheckCycle(cyclename, arot, brot, crot) {
        const xrot = Astronomy.CombineRotation(arot, brot);
        const irot = Astronomy.InverseRotation(xrot);
        CompareMatrices(cyclename, crot, irot, 2.0e-15);
    }

    function Test_RotRoundTrip() {
        const time = Astronomy.MakeTime(new Date('2067-05-30T14:45:00Z'));
        const observer = new Astronomy.Observer(+28, -82, 0);

        /*
            In each round trip, calculate a forward rotation and a backward rotation.
            Verify the two are inverse matrices.
        */

        /* Round trip #1: EQJ <==> EQD. */
        const eqj_eqd = Astronomy.Rotation_EQJ_EQD(time);
        const eqd_eqj = Astronomy.Rotation_EQD_EQJ(time);
        CheckInverse('eqj_eqd', 'eqd_eqj', eqj_eqd, eqd_eqj);

        /* Round trip #2: EQJ <==> ECL. */
        const eqj_ecl = Astronomy.Rotation_EQJ_ECL();
        const ecl_eqj = Astronomy.Rotation_ECL_EQJ();
        CheckInverse('eqj_ecl', 'ecl_eqj', eqj_ecl, ecl_eqj);

        /* Round trip #3: EQJ <==> HOR. */
        const eqj_hor = Astronomy.Rotation_EQJ_HOR(time, observer);
        const hor_eqj = Astronomy.Rotation_HOR_EQJ(time, observer);
        CheckInverse('eqj_hor', 'hor_eqj', eqj_hor, hor_eqj);

        /* Round trip #4: EQD <==> HOR. */
        const eqd_hor = Astronomy.Rotation_EQD_HOR(time, observer);
        const hor_eqd = Astronomy.Rotation_HOR_EQD(time, observer);
        CheckInverse('eqd_hor', 'hor_eqd', eqd_hor, hor_eqd);

        /* Round trip #5: EQD <==> ECL. */
        const eqd_ecl = Astronomy.Rotation_EQD_ECL(time);
        const ecl_eqd = Astronomy.Rotation_ECL_EQD(time);
        CheckInverse('eqd_ecl', 'ecl_eqd', eqd_ecl, ecl_eqd);

        /* Round trip #6: HOR <==> ECL. */
        const hor_ecl = Astronomy.Rotation_HOR_ECL(time, observer);
        const ecl_hor = Astronomy.Rotation_ECL_HOR(time, observer);
        CheckInverse('hor_ecl', 'ecl_hor', hor_ecl, ecl_hor);

        /* Round trip #7: EQD <==> ECT. */
        const eqd_ect = Astronomy.Rotation_EQD_ECT(time);
        const ect_eqd = Astronomy.Rotation_ECT_EQD(time);
        CheckInverse('eqd_ect', 'ect_eqd', eqd_ect, ect_eqd);

        /* Round trip #8: EQJ <==> ECT. */
        const eqj_ect = Astronomy.Rotation_EQJ_ECT(time);
        const ect_eqj = Astronomy.Rotation_ECT_EQJ(time);
        CheckInverse('eqj_ect', 'ect_eqj', eqj_ect, ect_eqj);

        // Verify that combining different sequences of rotations result
        // in the expected combination.
        // For example, (EQJ ==> HOR ==> ECL) must be the same matrix as (EQJ ==> ECL).
        CheckCycle('eqj_ecl, ecl_eqd, eqd_eqj', eqj_ecl, ecl_eqd, eqd_eqj);
        CheckCycle('eqj_hor, hor_ecl, ecl_eqj', eqj_hor, hor_ecl, ecl_eqj);
        CheckCycle('eqj_hor, hor_eqd, eqd_eqj', eqj_hor, hor_eqd, eqd_eqj);
        CheckCycle('ecl_eqd, eqd_hor, hor_ecl', ecl_eqd, eqd_hor, hor_ecl);
        CheckCycle('eqj_eqd, eqd_ect, ect_eqj', eqj_eqd, eqd_ect, ect_eqj);

        Debug('JS Test_RotRoundTrip: PASS');
    }

    function Rotation_Pivot() {
        const tolerance = 1.0e-15;

        /* Test #1 */

        /* Start with an identity matrix. */
        const ident = Astronomy.IdentityMatrix();

        /* Pivot 90 degrees counterclockwise around the z-axis. */
        let r = Astronomy.Pivot(ident, 2, +90.0);

        /* Put the expected answer in 'a'. */
        const a = Astronomy.MakeRotation([
            [ 0, +1,  0],
            [-1,  0,  0],
            [ 0,  0, +1],
        ]);

        /* Compare actual 'r' with expected 'a'. */
        CompareMatrices('Rotation_Pivot #1', r, a, tolerance);

        /* Test #2. */

        /* Pivot again, -30 degrees around the x-axis. */
        r = Astronomy.Pivot(r, 0, -30.0);

        /* Pivot a third time, 180 degrees around the y-axis. */
        r = Astronomy.Pivot(r, 1, +180.0);

        /* Use the 'r' matrix to rotate a vector. */
        const v1 = new Astronomy.Vector(1, 2, 3, Astronomy.MakeTime(0));

        const v2 = Astronomy.RotateVector(r, v1);

        /* Initialize the expected vector 've'. */
        const ve = new Astronomy.Vector(+2.0, +2.3660254037844390, -2.0980762113533156, v1.t);

        CompareVectors('Rotation_Pivot #2', v2, ve, tolerance);

        Debug('JS Rotation_Pivot: PASS');
    }


    Rotation_MatrixInverse();
    Rotation_MatrixMultiply();
    Rotation_Pivot();
    Test_GAL_EQJ_NOVAS('temp/galeqj.txt');

    Test_EQJ_EQD('Mercury');
    Test_EQJ_EQD('Venus');
    Test_EQJ_EQD('Mars');
    Test_EQJ_EQD('Jupiter');
    Test_EQJ_EQD('Saturn');

    Test_EQD_HOR('Mercury');
    Test_EQD_HOR('Venus');
    Test_EQD_HOR('Mars');
    Test_EQD_HOR('Jupiter');
    Test_EQD_HOR('Saturn');

    Test_EQD_ECT();

    Test_RotRoundTrip();

    console.log('JS Rotation: PASS');
    return 0;
}


function Refraction() {
    for (let alt = -90.1; alt <= +90.1; alt += 0.001) {
        const refr = Astronomy.Refraction('normal', alt);
        const corrected = alt + refr;
        const inv_refr = Astronomy.InverseRefraction('normal', corrected);
        const check_alt = corrected + inv_refr;
        const diff = abs(check_alt - alt);
        if (diff > 2.0e-14)
            throw `JS Refraction: alt=${alt}, refr=${refr}, diff=${diff}`;
    }

    console.log('JS Refraction: PASS');
    return 0;
}


function MonthNumber(mtext) {
    const index = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'].indexOf(mtext);
    if (index < 0)
        throw `Invalid month text "${mtext}"`;
    return 1 + index;
}


function Magnitude() {
    function LoadMagnitudeData(filename) {
        const lines = ReadLines(filename);
        let lnum = 0;
        let rows = [];
        for (let line of lines) {
            ++lnum;

            // [ Date__(UT)__HR:MN      APmag  S-brt               r        rdot            delta      deldot    S-T-O]
            // [ 2023-Mar-30 00:00      -4.01   1.17  0.719092953368  -0.1186373 1.20453495004726 -11.0204917  55.9004]
            let m = line.match(/^\s(\d{4})-(Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)-(\d{2})\s(\d{2}):(\d{2})\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s*$/);
            if (m) {
                const year = int(m[1]);
                const month = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'].indexOf(m[2]);
                const day = int(m[3]);
                const hour = int(m[4]);
                const minute = int(m[5]);

                const item = {
                    lnum: lnum,
                    date: new Date(Date.UTC(year, month, day, hour, minute)),
                    mag: float(m[6]),
                    sbrt: float(m[7]),
                    helio_dist: float(m[8]),
                    helio_radvel: float(m[9]),
                    geo_dist: float(m[10]),
                    geo_radvel: float(m[11]),
                    phase_angle: float(m[12])
                };

                rows.push(item);
            }
        }

        return {
            filename: filename,
            rows: rows
        };
    }

    function CheckMagnitudeData(body, data) {
        let diff_lo, diff_hi;
        let sum_squared_diff = 0;
        for (let item of data.rows) {
            let illum = Astronomy.Illumination(body, item.date);
            let diff = illum.mag - item.mag;
            sum_squared_diff += diff*diff;
            if (diff_lo === undefined) {
                diff_lo = diff_hi = diff;
            } else {
                diff_lo = min(diff_lo, diff);
                diff_hi = max(diff_hi, diff);
            }
        }
        let rms = sqrt(sum_squared_diff / data.rows.length);
        const limit = 0.012;
        const pass = (abs(diff_lo) < limit && abs(diff_hi) < limit);
        if (!pass || Verbose) console.log(`JS ${body.padEnd(8)} ${pass?"    ":"FAIL"}  diff_lo=${diff_lo.toFixed(4).padStart(8)}, diff_hi=${diff_hi.toFixed(4).padStart(8)}, rms=${rms.toFixed(4).padStart(8)}`);
        return pass;
    }

    function CheckSaturn() {
        // JPL Horizons does not include Saturn's rings in its magnitude models.
        // I still don't have authoritative test data for Saturn's magnitude.
        // For now, I just test for consistency with Paul Schlyter's formulas at:
        // http://www.stjarnhimlen.se/comp/ppcomp.html#15

        let success = true;

        const data = [
            { date: "1972-01-01T00:00Z", mag: -0.31725492,  tilt: +24.43386475 },
            { date: "1980-01-01T00:00Z", mag: +0.85796177,  tilt:  -1.72627324 },
            { date: "2009-09-04T00:00Z", mag: +1.01932560,  tilt:  +0.01834451 },
            { date: "2017-06-15T00:00Z", mag: -0.12303373,  tilt: -26.60068380 },
            { date: "2019-05-01T00:00Z", mag: +0.33124502,  tilt: -23.47173574 },
            { date: "2025-09-25T00:00Z", mag: +0.50543708,  tilt:  +1.69118986 },
            { date: "2032-05-15T00:00Z", mag: -0.04649573,  tilt: +26.95238680 }
        ];

        for (let item of data) {
            let illum = Astronomy.Illumination('Saturn', new Date(item.date));
            Debug(`JS Saturn: date=${illum.time.date.toISOString()}  mag=${illum.mag.toFixed(8).padStart(12)}  ring_tilt=${illum.ring_tilt.toFixed(8).padStart(12)}`);
            const mag_diff = abs(illum.mag - item.mag);
            if (mag_diff > 1.0e-8) {
                console.log(`ERROR: Excessive magnitude error ${mag_diff}`);
                success = false;
            }
            const tilt_diff = abs(illum.ring_tilt - item.tilt);
            if (tilt_diff > 1.0e-8) {
                console.log(`ERROR: Excessive ring tilt error ${tilt_diff}`);
                success = false;
            }
        }

        return success;
    }

    function TestMaxMag(filename, body) {
        // Test that we can find maximum magnitude events for Venus within
        // ranges found using JPL Horizons ephemeris data that has been
        // pre-processed by magnitude/findmax.py.

        const lines = ReadLines(filename);
        let date = new Date(Date.UTC(2000, 0, 1));
        let max_diff = 0;
        for (let line of lines) {
            let token = line.split(/\s+/);
            let date1 = new Date(token[0]);
            let date2 = new Date(token[1]);
            let evt = Astronomy.SearchPeakMagnitude(body, date);
            if (evt.time.date < date1 || evt.time.date > date2)
                throw `Event time ${evt.time.toString()} is outside the range ${date1.toISOString()} .. ${date2.toISOString()}`;

            // How close are we to the center date?
            let date_center = new Date((date1.getTime() + date2.getTime())/2);
            let diff_hours = abs(evt.time.date - date_center) / (1000 * 3600);
            if (diff_hours > 7.1)
                throw `Excessive diff_hours = ${diff_hours} from center date ${date_center.toISOString()}`;

            max_diff = max(max_diff, diff_hours);
            date = date2;
        }
        Debug(`JS TestMaxMag: ${lines.length} events, max error = ${max_diff.toFixed(3)} hours.`);
        return true;
    }

    let all_passed = true;
    const bodies = ['Sun', 'Moon', 'Mercury', 'Venus', 'Mars', 'Jupiter', 'Uranus', 'Neptune', 'Pluto'];
    for (let body of bodies) {
        const data = LoadMagnitudeData(`magnitude/${body}.txt`);
        if (!CheckMagnitudeData(body, data))
            all_passed = false;
    }

    if (!CheckSaturn())
        all_passed = false;

    if (!TestMaxMag('magnitude/maxmag_Venus.txt', 'Venus'))
        all_passed = false;

    all_passed || Fail('Found excessive error in at least one test.');
    console.log('JS Magnitude: PASS');
    return 0;
}


function Constellation() {
    const filename = 'constellation/test_input.txt';
    const lines = ReadLines(filename);
    let lnum = 0;
    let failcount = 0;
    for (let row of lines) {
        ++lnum;
        let token = row.trim().split(/\s+/);
        if (token.length !== 4) {
            Fail(`Bad data in ${filename} line ${lnum}: found ${token.length} tokens.`);
        }
        const id = int(token[0]);
        const ra = float(token[1]);
        const dec = float(token[2]);
        const symbol = token[3];
        if (!/^[A-Z][A-Za-z]{2}$/.test(symbol)) {
            Fail(`Invalid symbol "${symbol}" in ${filename} line ${lnum}`);
        }
        const constel = Astronomy.Constellation(ra, dec);
        if (constel.symbol !== symbol) {
            ++failcount;
            console.error(`Star ${id}: expected ${symbol}, found ${constel.symbol} at B1875 RA=${constel.ra1875}, DEC=${constel.dec1875}`);
        }
    }
    if (failcount > 0) {
        Fail(`JS Constellation: ${failcount} failures.`);
    }
    console.log(`JS Constellation: PASS (verified ${lines.length})`);
    return 0;
}


function TransitFile(body, filename, limit_minutes, limit_sep) {
    const lines = ReadLines(filename);
    let lnum = 0;
    let max_minutes = 0;
    let max_sep = 0;
    let transit = Astronomy.SearchTransit(body, Astronomy.MakeTime(new Date(Date.UTC(1600, 0))));
    for (let row of lines) {
        ++lnum;
        let token = row.trim().split(/\s+/);

        /* 22:17 1881-11-08T00:57Z 03:38  3.8633 */
        if (token.length !== 4) {
            Fail(`JS TransitFile(${filename} line ${lnum}): bad data format.`);
        }
        const textp = token[1];
        const text1 = textp.substr(0, 11) + token[0] + 'Z';
        const text2 = textp.substr(0, 11) + token[2] + 'Z';
        let timep = Astronomy.MakeTime(new Date(textp));
        let time1 = Astronomy.MakeTime(new Date(text1));
        let time2 = Astronomy.MakeTime(new Date(text2));
        const separation = float(token[3]);

        // If the start time is after the peak time, it really starts on the previous day.
        if (time1.ut > timep.ut)
            time1 = time1.AddDays(-1.0);

        // If the finish time is before the peak time, it really starts on the following day.
        if (time2.ut < timep.ut)
            time2 = time2.AddDays(+1.0);

        const diff_start  = (24.0 * 60.0) * abs(time1.ut - transit.start.ut );
        const diff_peak   = (24.0 * 60.0) * abs(timep.ut - transit.peak.ut  );
        const diff_finish = (24.0 * 60.0) * abs(time2.ut - transit.finish.ut);
        const diff_sep = abs(separation - transit.separation);

        max_minutes = max(max_minutes, diff_start);
        max_minutes = max(max_minutes, diff_peak);
        max_minutes = max(max_minutes, diff_finish);
        if (max_minutes > limit_minutes) {
            Fail(`JS TransitFile(${filename} line ${lnum}): EXCESSIVE TIME ERROR = ${max_minutes} minutes.`);
        }

        max_sep = max(max_sep, diff_sep);
        if (max_sep > limit_sep) {
            Fail(`JS TransitFile(${filename} line ${lnum}): EXCESSIVE SEPARATION ERROR = ${max_sep} arcminutes.`);
        }

        transit = Astronomy.NextTransit(body, transit.finish);
    }

    console.log(`JS TransitFile(${filename}): PASS - verified ${lnum}, max minutes = ${max_minutes}, max sep arcmin = ${max_sep}`);
    return 0;
}


function Transit() {
    if (0 !== TransitFile('Mercury', 'eclipse/mercury.txt', 10.710, 0.2121)) {
        return 1;
    }

    if (0 !== TransitFile('Venus', 'eclipse/venus.txt', 9.109, 0.6772)) {
        return 1;
    }

    return 0;
}


function PlutoCheckDate(ut, arcmin_tolerance, x, y, z) {
    const time = Astronomy.MakeTime(ut);
    const timeText = time.toString();
    Debug(`JS PlutoCheck: ${timeText} = ${time.ut} UT = ${time.tt} TT`);
    const vector = Astronomy.HelioVector('Pluto', time);
    const dx = v(vector.x - x);
    const dy = v(vector.y - y);
    const dz = v(vector.z - z);
    const diff = sqrt(dx*dx + dy*dy + dz*dz);
    const dist = (sqrt(x*x + y*y + z*z) - 1.0);       /* worst-case distance between Pluto and Earth */
    const arcmin = (diff / dist) * (180.0 * 60.0 / Math.PI);
    Debug(`JS PlutoCheck: calc pos = [${vector.x}, ${vector.y}, ${vector.z}]`);
    Debug(`JS PlutoCheck: ref  pos = [${x}, ${y}, ${z}]`);
    Debug(`JS PlutoCheck: del  pos = [${vector.x - x}, ${vector.y - y}, ${vector.z - z}]`);
    Debug(`JS PlutoCheck: diff = ${diff} AU, ${arcmin} arcmin`);
    Debug("");
    if (v(arcmin) > arcmin_tolerance) {
        console.error("JS PlutoCheck: EXCESSIVE ERROR");
        return 1;
    }
    return 0;
}


function PlutoCheck() {
    if (0 !== PlutoCheckDate(  +18250.0,  0.089, +37.4377303523676090, -10.2466292454075898, -14.4773101310875809)) return 1;
    if (0 !== PlutoCheckDate( -856493.0,  4.067, +23.4292113199166252, +42.1452685817740829,  +6.0580908436642940)) return 1;
    if (0 !== PlutoCheckDate( +435633.0,  0.016, -27.3178902095231813, +18.5887022581070305, +14.0493896259306936)) return 1;
    if (0 !== PlutoCheckDate(       0.0,  8.e-9,  -9.8753673425269000, -27.9789270580402771,  -5.7537127596369588)) return 1;
    if (0 !== PlutoCheckDate( +800916.0,  2.286, -29.5266052645301365, +12.0554287322176474, +12.6878484911631091)) return 1;
    console.log("JS PlutoCheck: PASS");
    return 0;
}


function GeoidTestCase(time, observer, ofdate) {
    let topo_moon = Astronomy.Equator(Astronomy.Body.Moon, time, observer, ofdate, false);
    let surface = Astronomy.ObserverVector(time, observer, ofdate);
    let geo_moon = Astronomy.GeoVector(Astronomy.Body.Moon, time, false);

    if (ofdate) {
        // GeoVector() returns J2000 coordinates. Convert to equator-of-date coordinates.
        const rot = Astronomy.Rotation_EQJ_EQD(time);
        geo_moon = Astronomy.RotateVector(rot, geo_moon);
    }

    const dx = Astronomy.KM_PER_AU * v((geo_moon.x - surface.x) - topo_moon.vec.x);
    const dy = Astronomy.KM_PER_AU * v((geo_moon.y - surface.y) - topo_moon.vec.y);
    const dz = Astronomy.KM_PER_AU * v((geo_moon.z - surface.z) - topo_moon.vec.z);
    const diff = sqrt(dx*dx + dy*dy + dz*dz);
    Debug(`JS GeoidTestCase: ofdate=${ofdate}, time=${time.toISOString()}, lat=${observer.latitude}, lon=${observer.longitude}, ht=${observer.height}, surface=(${Astronomy.KM_PER_AU * surface.x}, ${Astronomy.KM_PER_AU * surface.y}, ${Astronomy.KM_PER_AU * surface.z}), diff = ${diff} km`);

    // Require 1 millimeter accuracy! (one millionth of a kilometer).
    if (diff > 1.0e-6) {
        console.error('JS GeoidTestCase: EXCESSIVE POSITION ERROR.');
        return 1;
    }

    // Verify that we can convert the surface vector back to an observer.
    const vobs = Astronomy.VectorObserver(surface, ofdate);
    const lat_diff = abs(vobs.latitude - observer.latitude);
    let lon_diff;

    /* Longitude is meaningless at the poles, so don't bother checking it there. */
    if (-89.99 <= observer.latitude && observer.latitude <= +89.99) {
        lon_diff = abs(vobs.longitude - observer.longitude);
        if (lon_diff > 180.0)
            lon_diff = 360.0 - lon_diff;
        lon_diff = abs(lon_diff * Math.cos(Astronomy.DEG2RAD * observer.latitude));
        if (lon_diff > 1.0e-6) {
            console.error(`JS GeoidTestCase: EXCESSIVE longitude check error = ${lon_diff}`);
            return 1;
        }
    } else {
        lon_diff = 0.0;
    }

    const h_diff = abs(vobs.height - observer.height);
    Debug(`JS GeoidTestCase: vobs=(lat=${vobs.latitude}, lon=${vobs.longitude}, height=${vobs.height}), lat_diff=${lat_diff}, lon_diff=${lon_diff}, h_diff=${h_diff}`);

    if (lat_diff > 1.0e-6) {
        console.error(`JS GeoidTestCase: EXCESSIVE latitude check error = ${lat_diff}`);
        return 1;
    }

    if (h_diff > 0.001) {
        console.error(`JS GeoidTestCase: EXCESSIVE height check error = ${h_diff}`);
        return 1;
    }

    return 0;
}


function Geoid() {
    const time_list = [
        new Date('1066-09-27T18:00:00Z'),
        new Date('1970-12-13T15:42:00Z'),
        new Date('1970-12-13T15:43:00Z'),
        new Date('2015-03-05T02:15:45Z')
    ];

    const observer_list = [
        new Astronomy.Observer( +1.5,   +2.7,    7.4),
        new Astronomy.Observer(-53.7, +141.7, +100.0),
        new Astronomy.Observer(+30.0,  -85.2,  -50.0),
        new Astronomy.Observer(+90.0,  +45.0,  -50.0),
        new Astronomy.Observer(-90.0, -180.0,    0.0)
    ];

    // Test a variety of times and locations, in both supported orientation systems.

    let observer, time;
    for (observer of observer_list) {
        for (time of time_list) {
            if (0 != GeoidTestCase(time, observer, false))
                return 1;
            if (0 != GeoidTestCase(time, observer, true))
                return 1;
        }
    }

    // More exhaustive tests for a single time value across many different geographic coordinates.

    time = new Date('2021-06-20T15:08:00Z');
    for (let lat = -90; lat <= +90; lat += 1) {
        for (let lon = -175; lon <= +180; lon += 5) {
            observer = new Astronomy.Observer(lat, lon, 0.0);
            if (0 != GeoidTestCase(time, observer, true))
                return 1;
        }
    }

    console.log('JS GeoidTest: PASS');
    return 0;
}


function JupiterMoons_CheckJpl(mindex, tt, pos, vel) {
    const pos_tolerance = 9.0e-4;
    const vel_tolerance = 9.0e-4;
    const time = Astronomy.AstroTime.FromTerrestrialTime(tt);
    const jm = Astronomy.JupiterMoons(time);
    const moon = SelectJupiterMoon(jm, mindex);

    let dx = v(pos[0] - moon.x);
    let dy = v(pos[1] - moon.y);
    let dz = v(pos[2] - moon.z);
    let mag = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
    const pos_diff = sqrt(dx*dx + dy*dy + dz*dz) / mag;
    if (pos_diff > pos_tolerance) {
        console.error(`JS JupiterMoons_CheckJpl(mindex=${mindex}, tt=${tt}): excessive position error ${pos_diff}`);
        return 1;
    }

    dx = v(vel[0] - moon.vx);
    dy = v(vel[1] - moon.vy);
    dz = v(vel[2] - moon.vz);
    mag = sqrt(vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
    const vel_diff = sqrt(dx*dx + dy*dy + dz*dz) / mag;
    if (vel_diff > vel_tolerance) {
        console.error(`JS JupiterMoons_CheckJpl(mindex=${mindex}, tt=${tt}): excessive velocity error ${vel_diff}`);
        return 1;
    }

    Debug(`JS JupiterMoons_CheckJpl: mindex=${mindex}, tt=${tt}, pos_diff=${pos_diff}, vel_diff=${vel_diff}`);
    return 0;
}


class JplStateRecord {
    constructor(lnum, state) {
        this.lnum = lnum;
        this.state = state;
    }
}

function JplHorizonsStateVectors(filename) {
    const lines = ReadLines(filename);
    let lnum = 0;
    let found = false;
    let part = -1;
    let tt, match, pos, vel, time;
    let reclist = [];
    for (let line of lines) {
        ++lnum;
        if (!found) {
            if (line == '$$SOE') {
                found = true;
                part = 0;
            }
        } else if (line == '$$EOE') {
            break;
        } else {
            switch (part) {
            case 0:
                // 2446545.000000000 = A.D. 1986-Apr-24 12:00:00.0000 TDB
                tt = float(line.split()[0]) - 2451545.0;    // convert JD to J2000 TT
                time = Astronomy.AstroTime.FromTerrestrialTime(tt);
                break;

            case 1:
                // X = 1.134408131605554E-03 Y =-2.590904586750408E-03 Z =-7.490427225904720E-05
                match = /\s*X =\s*(\S+) Y =\s*(\S+) Z =\s*(\S+)/.exec(line);
                if (!match) {
                    console.error(`JS JplHorizonsStateVectors(${filename} line ${lnum}): cannot parse position vector.`);
                    return 1;
                }
                pos = [ float(match[1]), float(match[2]), float(match[3]) ];
                break;

            case 2:
                // VX= 9.148038778472862E-03 VY= 3.973823407182510E-03 VZ= 2.765660368640458E-04
                match = /\s*VX=\s*(\S+) VY=\s*(\S+) VZ=\s*(\S+)/.exec(line);
                if (!match) {
                    console.error(`JS JplHorizonsStateVectors(${filename} line ${lnum}): cannot parse velocity vector.`);
                    return 1;
                }
                vel = [ float(match[1]), float(match[2]), float(match[3]) ];
                let state = new Astronomy.StateVector(pos[0], pos[1], pos[2], vel[0], vel[1], vel[2], time);
                reclist.push(new JplStateRecord(lnum, state));
                break;

            default:
                console.error(`JS JplHorizonsStateVectors(${filename} line ${lnum}): unexpected part = ${part}`);
                return 1;
            }
            part = (part + 1) % 3;
        }
    }
    return reclist;
}



function JupiterMoons() {
    for (let mindex = 0; mindex < 4; ++mindex) {
        const filename = `jupiter_moons/horizons/jm${mindex}.txt`;
        const lines = ReadLines(filename);
        let lnum = 0;
        let found = false;
        let part = -1;
        const expected_count = 5001;
        let count = 0;
        let tt, match, pos, vel;
        for (let line of lines) {
            ++lnum;
            if (!found) {
                if (line == '$$SOE') {
                    found = true;
                    part = 0;
                } else if (line.startsWith('Revised:')) {
                    if (line.length !== 79) {
                        console.error(`JS JupiterMoons(${filename} line ${lnum}): unexpected line length.`);
                        return 1;
                    }
                    console.log(line.substr(76));
                    const check_mindex = int(line.substr(76)) - 501;
                    if (mindex !== check_mindex) {
                        console.error(`JS JupiterMoons(${filename} line ${lnum}): moon index does not match: check=${check_mindex}, mindex=${mindex}.`);
                        return 1;
                    }
                }
            } else if (line == '$$EOE') {
                break;
            } else {
                switch (part) {
                case 0:
                    // 2446545.000000000 = A.D. 1986-Apr-24 12:00:00.0000 TDB
                    tt = float(line.split()[0]) - 2451545.0;    // convert JD to J2000 TT
                    break;

                case 1:
                    // X = 1.134408131605554E-03 Y =-2.590904586750408E-03 Z =-7.490427225904720E-05
                    match = /\s*X =\s*(\S+) Y =\s*(\S+) Z =\s*(\S+)/.exec(line);
                    if (!match) {
                        console.error(`JS JupiterMoons(${filename} line ${lnum}): cannot parse position vector.`);
                        return 1;
                    }
                    pos = [ float(match[1]), float(match[2]), float(match[3]) ];
                    break;

                case 2:
                    // VX= 9.148038778472862E-03 VY= 3.973823407182510E-03 VZ= 2.765660368640458E-04
                    match = /\s*VX=\s*(\S+) VY=\s*(\S+) VZ=\s*(\S+)/.exec(line);
                    if (!match) {
                        console.error(`JS JupiterMoons(${filename} line ${lnum}): cannot parse velocity vector.`);
                        return 1;
                    }
                    vel = [ float(match[1]), float(match[2]), float(match[3]) ];
                    if (JupiterMoons_CheckJpl(mindex, tt, pos, vel)) {
                        console.error(`JS JupiterMoons(${filename} line ${lnum}): FAILED VERIFICATION.`);
                        return 1;
                    }
                    ++count;
                    break;

                default:
                    console.error(`JS JupiterMoons(${filename} line ${lnum}): unexpected part = ${part}`);
                    return 1;
                }
                part = (part + 1) % 3;
            }
        }
        if (count !== expected_count) {
            console.error(`JS JupiterMoons: Expected ${expected_count} test cases, but found ${count}`);
            return 1;
        }
    }
    console.log(`JS JupiterMoons: PASS`);
    return 0;
}


function Issue103() {
    // https://github.com/cosinekitty/astronomy/issues/103

    const observer = new Astronomy.Observer(29, -81, 10);
    const ut = -51279.9420508868643083;
    const time = Astronomy.MakeTime(ut);
    const ofdate = Astronomy.Equator(Astronomy.Body.Moon, time, observer, true, true);

    console.log(`tt  = ${time.tt.toFixed(16)}`);
    console.log(`ra  = ${ofdate.ra.toFixed(16)}`);
    console.log(`dec = ${ofdate.dec.toFixed(16)}`);
    const hor = Astronomy.Horizon(time, observer, ofdate.ra, ofdate.dec, false);
    console.log(`az  = ${hor.azimuth.toFixed(16)}`);
    console.log(`alt = ${hor.altitude.toFixed(16)}`);

    return 0;
}


function AberrationTest() {
    const THRESHOLD_SECONDS = 0.453;
    const filename = 'equatorial/Mars_j2000_ofdate_aberration.txt';
    const lines = ReadLines(filename);
    let lnum = 0;
    let found_begin = false;
    let max_diff_seconds = 0;
    let count = 0;
    for (let line of lines) {
        ++lnum;
        if (!found_begin) {
            if (line === '$$SOE')
                found_begin = true;
        } else if (line === '$$EOE') {
            break;
        } else {
            // 2459371.500000000 *   118.566080210  22.210647456 118.874086738  22.155784122
            const ut = float(line.trim().split(/\s+/)[0]) - 2451545.0;    // convert JD to J2000 UT day value
            const time = Astronomy.MakeTime(ut);
            const tokens = line.substr(22).trim().split(/\s+/);
            const jra = float(tokens[0]);
            const jdec = float(tokens[1]);
            const dra = float(tokens[2]);
            const ddec = float(tokens[3]);

            // Create spherical coordinates with magnitude = speed of light.
            // This helps us perform the aberration correction later.
            const eqj_sphere = new Astronomy.Spherical(jdec, jra, Astronomy.C_AUDAY);

            // Convert EQJ angular coordinates (jra, jdec) to an EQJ vector.
            // This vector represents a ray of light, only it is travelling from the observer toward the star.
            // The backwards direction makes the math simpler.
            const eqj_vec = Astronomy.VectorFromSphere(eqj_sphere, time);

            // Calculate the Earth's barycentric velocity vector in EQJ coordinates.
            const eqj_earth = Astronomy.BaryState(Astronomy.Body.Earth, time);

            // Use non-relativistic approximation: add light vector to Earth velocity vector.
            // This gives aberration-corrected apparent position of the star in EQJ.
            eqj_vec.x += eqj_earth.vx;
            eqj_vec.y += eqj_earth.vy;
            eqj_vec.z += eqj_earth.vz;

            // Calculate the rotation matrix that converts J2000 coordinates to of-date coordinates.
            const rot = Astronomy.Rotation_EQJ_EQD(time);

            // Use the rotation matrix to re-orient the EQJ vector to a EQD vector.
            const eqd_vec = Astronomy.RotateVector(rot, eqj_vec);

            // Convert the EQD vector back to spherical angular coordinates.
            const eqd_sphere = Astronomy.SphereFromVector(eqd_vec);

            // Calculate the differences in RA and DEC between expected and calculated values.
            const factor = cos(eqd_sphere.lat * Astronomy.DEG2RAD);     // RA errors are less important toward the poles.
            const xra = factor * abs(eqd_sphere.lon - dra);
            const xdec = abs(eqd_sphere.lat - ddec);
            const diff_seconds = 3600 * sqrt(xra*xra + xdec*xdec);
            Debug(`JS AberrationTest(${filename} line ${lnum}): xra=${xra.toFixed(6)} deg, xdec=${xdec.toFixed(6)} deg, diff_seconds=${diff_seconds.toFixed(3)}`);
            if (diff_seconds > THRESHOLD_SECONDS) {
                console.error(`JS AberrationTest(${filename} line ${lnum}): EXCESSIVE ANGULAR ERROR = ${diff_seconds.toFixed(3)} seconds.`);
                return 1;
            }
            if (diff_seconds > max_diff_seconds)
                max_diff_seconds = diff_seconds;
            ++count;
        }
    }
    console.log(`JS AberrationTest(${filename}): PASS - Tested ${count} cases. max_diff_seconds = ${max_diff_seconds.toFixed(3)}`);
    return 0;
}


function StateVectorDiff(relative, vec, x, y, z) {
    const dx = v(vec[0] - x);
    const dy = v(vec[1] - y);
    const dz = v(vec[2] - z);
    let diff_squared = dx*dx + dy*dy + dz*dz;
    if (relative)
        diff_squared /= (vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
    return sqrt(diff_squared);
}


function VerifyState(func, score, filename, lnum, time, pos, vel, r_thresh, v_thresh) {
    const state = func.Eval(time);

    const rdiff = StateVectorDiff((r_thresh > 0.0), pos, state.x, state.y, state.z);
    if (rdiff > score.max_rdiff)
        score.max_rdiff = rdiff;

    const vdiff = StateVectorDiff((v_thresh > 0.0), vel, state.vx, state.vy, state.vz);
    if (vdiff > score.max_vdiff)
        score.max_vdiff = vdiff;

    if (rdiff > Math.abs(r_thresh)) {
        console.error(`JS VerifyState(${filename} line ${lnum}): EXCESSIVE POSITION ERROR = ${rdiff.toExponential(3)}`);
        return 1;
    }

    if (vdiff > Math.abs(v_thresh)) {
        console.error(`JS VerifyState(${filename} line ${lnum}): EXCESSIVE VELOCITY ERROR = ${vdiff.toExponential(3)}`);
        return 1;
    }

    return 0;
}


function VerifyStateBody(func, filename, r_thresh, v_thresh) {
    let score = { max_rdiff:0, max_vdiff:0 };
    let count = 0;

    for (let rec of JplHorizonsStateVectors(filename)) {
        let pos = [rec.state.x, rec.state.y, rec.state.z];
        let vel = [rec.state.vx, rec.state.vy, rec.state.vz];
        if (VerifyState(func, score, filename, rec.lnum, rec.state.t, pos, vel, r_thresh, v_thresh))
            return 1;
        ++count;
    }

    Debug(`JS VerifyStateBody(${filename}): PASS - Tested ${count} cases. max rdiff=${score.max_rdiff.toExponential(3)}, vdiff=${score.max_vdiff.toExponential(3)}`);
    return 0;
}


// Constants for use inside unit tests only; they doesn't make sense for public consumption.
const Body_GeoMoon = -100;
const Body_Geo_EMB = -101;

class BaryStateFunc {
    constructor(body) {
        this.body = body;
    }

    Eval(time) {
        if (this.body === Body_GeoMoon)
            return Astronomy.GeoMoonState(time);

        if (this.body === Body_Geo_EMB)
            return Astronomy.GeoEmbState(time);

        return Astronomy.BaryState(this.body, time);
    }
}


function BaryStateTest() {
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Sun),     'barystate/Sun.txt',     -1.224e-05, -1.134e-07)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Mercury), 'barystate/Mercury.txt',  1.672e-04,  2.698e-04)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Venus),   'barystate/Venus.txt',    4.123e-05,  4.308e-05)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Earth),   'barystate/Earth.txt',    2.296e-05,  6.359e-05)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Mars),    'barystate/Mars.txt',     3.107e-05,  5.550e-05)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Jupiter), 'barystate/Jupiter.txt',  7.389e-05,  2.471e-04)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Saturn),  'barystate/Saturn.txt',   1.067e-04,  3.220e-04)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Uranus),  'barystate/Uranus.txt',   9.035e-05,  2.519e-04)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Neptune), 'barystate/Neptune.txt',  9.838e-05,  4.446e-04)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Pluto),   'barystate/Pluto.txt',    4.259e-05,  7.827e-05)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.Moon),    "barystate/Moon.txt",     2.354e-05,  6.604e-05)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Astronomy.Body.EMB),     "barystate/EMB.txt",      2.353e-05,  6.511e-05)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Body_GeoMoon),           "barystate/GeoMoon.txt",  4.086e-05,  5.347e-05)) return 1;
    if (VerifyStateBody(new BaryStateFunc(Body_Geo_EMB),           "barystate/GeoEMB.txt",   4.076e-05,  5.335e-05)) return 1;
    console.log('JS BaryStateTest: PASS');
    return 0;
}


class HelioStateFunc {
    constructor(body) {
        this.body = body;
    }

    Eval(time) {
        return Astronomy.HelioState(this.body, time);
    }
}


function HelioStateTest() {
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.SSB),     'heliostate/SSB.txt',     -1.209e-05, -1.125e-07)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Mercury), 'heliostate/Mercury.txt',  1.481e-04,  2.756e-04)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Venus),   'heliostate/Venus.txt',    3.528e-05,  4.485e-05)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Earth),   'heliostate/Earth.txt',    1.476e-05,  6.105e-05)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Mars),    'heliostate/Mars.txt',     3.154e-05,  5.603e-05)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Jupiter), 'heliostate/Jupiter.txt',  7.455e-05,  2.562e-04)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Saturn),  'heliostate/Saturn.txt',   1.066e-04,  3.150e-04)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Uranus),  'heliostate/Uranus.txt',   9.034e-05,  2.712e-04)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Neptune), 'heliostate/Neptune.txt',  9.834e-05,  4.534e-04)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Pluto),   'heliostate/Pluto.txt',    4.271e-05,  1.198e-04)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.Moon),    'heliostate/Moon.txt',     1.477e-05,  6.195e-05)) return 1;
    if (VerifyStateBody(new HelioStateFunc(Astronomy.Body.EMB),     'heliostate/EMB.txt',      1.476e-05,  6.106e-05)) return 1;
    console.log('JS HelioStateTest: PASS');
    return 0;
}


class TopoStateFunc {
    constructor(body) {
        this.body = body;
    }

    Eval(time) {
        const observer = new Astronomy.Observer(30.0, -80.0, 1000.0);

        let observer_state = Astronomy.ObserverState(time, observer, false);
        let state;
        if (this.body == Body_Geo_EMB) {
            state = Astronomy.GeoEmbState(time);
        } else if (this.body == Astronomy.Body.Earth) {
            state = new Astronomy.StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time);
        } else {
            throw `JS TopoStateFunction: unsupported body ${this.body}`;
        }

        state.x  -= observer_state.x;
        state.y  -= observer_state.y;
        state.z  -= observer_state.z;
        state.vx -= observer_state.vx;
        state.vy -= observer_state.vy;
        state.vz -= observer_state.vz;

        return state;
    }
}

function TopoStateTest() {
    if (VerifyStateBody(new TopoStateFunc(Astronomy.Body.Earth),  "topostate/Earth_N30_W80_1000m.txt",  2.108e-04, 2.430e-04)) return 1;
    if (VerifyStateBody(new TopoStateFunc(Body_Geo_EMB),          "topostate/EMB_N30_W80_1000m.txt",    7.197e-04, 2.497e-04)) return 1;
    console.log("JS TopoStateTest: PASS");
    return 0;
}



function TwilightTest() {
    const tolerance_seconds = 60.0;
    const filename = 'riseset/twilight.txt';
    const lines = ReadLines(filename);
    let lnum = 0;
    let max_diff = 0.0;

    const name = [
        "astronomical dawn",
        "nautical dawn",
        "civil dawn",
        "civil dusk",
        "nautical dusk",
        "astronomical dusk"
    ];

    for (let line of lines) {
        ++lnum;
        const tokens = line.split(/\s+/);
        if (tokens.length !== 9) {
            console.error(`JS TwilightTest: FAIL(${filename} line ${lnum}): incorrect number of tokens = ${tokens.length}.`);
            return 1;
        }
        const lat = float(tokens[0]);
        const lon = float(tokens[1]);
        const observer = new Astronomy.Observer(lat, lon, 0);

        const searchDate = Astronomy.MakeTime(new Date(tokens[2]));
        const correctTimes = [
            Astronomy.MakeTime(new Date(tokens[3])),      // astronomical dawn
            Astronomy.MakeTime(new Date(tokens[4])),      // nautical dawn
            Astronomy.MakeTime(new Date(tokens[5])),      // civil dawn
            Astronomy.MakeTime(new Date(tokens[6])),      // civil dusk
            Astronomy.MakeTime(new Date(tokens[7])),      // nautical dusk
            Astronomy.MakeTime(new Date(tokens[8]))       // astronomical dusk
        ];

        const calcTimes = [
            Astronomy.SearchAltitude(Astronomy.Body.Sun, observer, +1, searchDate, 1, -18),   // astronomical dawn
            Astronomy.SearchAltitude(Astronomy.Body.Sun, observer, +1, searchDate, 1, -12),   // nautical dawn
            Astronomy.SearchAltitude(Astronomy.Body.Sun, observer, +1, searchDate, 1, -6),    // civil dawn
            Astronomy.SearchAltitude(Astronomy.Body.Sun, observer, -1, searchDate, 1, -6),    // civil dusk
            Astronomy.SearchAltitude(Astronomy.Body.Sun, observer, -1, searchDate, 1, -12),   // nautical dusk
            Astronomy.SearchAltitude(Astronomy.Body.Sun, observer, -1, searchDate, 1, -18),   // astronomical dusk
        ];

        for (let i = 0; i < correctTimes.length; ++i) {
            const correct = correctTimes[i];
            const calc = calcTimes[i];
            const diff = SECONDS_PER_DAY * abs(calc.ut - correct.ut);
            if (diff > tolerance_seconds) {
                console.error(`JS TwilightTest(${filename} line ${lnum}): EXCESSIVE ERROR = ${diff} seconds for ${name[i]}`);
                console.error(`Expected ${correct} but calculated ${calc}`);
                return 1;
            }
            if (diff > max_diff)
                max_diff = diff;
        }
    }
    console.log(`JS TwilightTest: PASS (${lnum} test cases, max error = ${max_diff.toFixed(3)} seconds)`);
    return 0;
}


function VectorObserverCase(body) {
    for (let ts = 1717780096000; ts <= 1717780096200; ts++) {
        var date = new Date(ts);
        var vect = Astronomy.GeoVector(body, date, true)
        //console.log('ts, date, vector', ts, date.valueOf(), vect)
        var point = Astronomy.VectorObserver(vect)
        //console.log('point', point)
        //console.log('*')
     }
     return Pass(`VectorObserverCase(${body})`);
 }


function VectorObserverIssue347() {
    // https://github.com/cosinekitty/astronomy/issues/347
    return (
        VectorObserverCase(Astronomy.Body.Sun) ||
        VectorObserverCase(Astronomy.Body.Jupiter) ||
        VectorObserverCase(Astronomy.Body.Saturn) ||
        VectorObserverCase(Astronomy.Body.Uranus) ||
        VectorObserverCase(Astronomy.Body.Neptune)
    );
}


function Libration(filename) {
    const lines = ReadLines(filename);
    let max_diff_elon = 0.0;
    let max_diff_elat = 0.0;
    let max_diff_distance = 0.0;
    let max_diff_diam = 0.0;
    let max_eclip_lon = -900.0;
    let count = 0;
    let lnum = 0;
    for (let line of lines) {
        ++lnum;
        if (lnum === 1) {
            if (line !== "   Date       Time    Phase    Age    Diam    Dist     RA        Dec      Slon      Slat     Elon     Elat   AxisA") {
                console.error(`JS Libration(${filename} line ${lnum}): unexpected header line.`);
                return 1;
            }
        } else {
            const token = line.split(/\s+/);
            if (token.length !== 16) {
                console.error(`JS Libration: FAIL(${filename} line ${lnum}): incorrect number of tokens = ${token.length}.`);
                return 1;
            }

            const day = int(token[0]);
            const month = MonthNumber(token[1]);
            const year = int(token[2]);
            const hmtoken = token[3].split(':');
            if (hmtoken.length !== 2) {
                console.error(`JS Libration(${filename} line ${lnum}): expected hh:mm but found '${token[3]}'`);
                return 1;
            }
            const hour = int(hmtoken[0]);
            const minute = int(hmtoken[1]);
            const time = Astronomy.MakeTime(new Date(Date.UTC(year, month-1, day, hour, minute)));

            const diam = float(token[7]) / 3600.0;
            const dist = float(token[8]);
            const elon = float(token[13]);
            const elat = float(token[14]);

            const lib = Astronomy.Libration(time);

            const diff_elon = 60.0 * abs(lib.elon - elon);
            if (diff_elon > max_diff_elon)
                max_diff_elon = diff_elon;

            const diff_elat = 60.0 * abs(lib.elat - elat);
            if (diff_elat > max_diff_elat)
                max_diff_elat = diff_elat;

            const diff_distance = abs(lib.dist_km - dist);
            if (diff_distance > max_diff_distance)
                max_diff_distance = diff_distance;

            const diff_diam = abs(lib.diam_deg - diam);
            if (diff_diam > max_diff_diam)
                max_diff_diam = diff_diam;

            if (lib.mlon > max_eclip_lon)
                max_eclip_lon = lib.mlon;

            if (diff_elon > 0.1304) {
                console.error(`JS Libration(${filename} line ${lnum}): EXCESSIVE diff_elon = ${diff_elon} arcmin`);
                return 1;
            }

            if (diff_elat > 1.6476) {
                console.error(`JS Libration(${filename} line ${lnum}): EXCESSIVE diff_elat = ${diff_elat} arcmin`);
                return 1;
            }

            if (diff_distance > 54.377) {
                console.error(`JS Libration(${filename} line ${lnum}): EXCESSIVE diff_distance = ${diff_distance} km`);
                return 1;
            }

            if (diff_diam > 0.00009) {
                console.error(`JS Libration(${filename}): EXCESSIVE diff_diam = ${diff_diam} degrees.`);
                return 1;
            }
            ++count;
        }
    }
    if (max_eclip_lon < 359.0 || max_eclip_lon > 360.0) {
        console.error(`JS Libration(${filename}): INVALID max ecliptic longitude = ${max_eclip_lon.toFixed(3)} degrees.`);
        return 1;
    }
    console.log(`JS Libration(${filename}): PASS (${count} test cases, max_diff_elon = ${max_diff_elon} arcmin, max_diff_elat = ${max_diff_elat} arcmin, max_diff_distance = ${max_diff_distance} km, max_diff_diam = ${max_diff_diam} deg)`);
    return 0;
}


function LibrationTest() {
    return (
        Libration("libration/mooninfo_2020.txt") ||
        Libration("libration/mooninfo_2021.txt") ||
        Libration("libration/mooninfo_2022.txt")
    );
}


function AxisTest() {
    return (
        AxisTestBody(Astronomy.Body.Sun,      "axis/Sun.txt",       0.0)        ||
        AxisTestBody(Astronomy.Body.Mercury,  "axis/Mercury.txt",   0.074340)   ||
        AxisTestBody(Astronomy.Body.Venus,    "axis/Venus.txt",     0.0)        ||
        AxisTestBody(Astronomy.Body.Earth,    "axis/Earth.txt",     0.002033)   ||
        AxisTestBody(Astronomy.Body.Moon,     "axis/Moon.txt",      0.264845)   ||
        AxisTestBody(Astronomy.Body.Mars,     "axis/Mars.txt",      0.075323)   ||
        AxisTestBody(Astronomy.Body.Jupiter,  "axis/Jupiter.txt",   0.000324)   ||
        AxisTestBody(Astronomy.Body.Saturn,   "axis/Saturn.txt",    0.000304)   ||
        AxisTestBody(Astronomy.Body.Uranus,   "axis/Uranus.txt",    0.0)        ||
        AxisTestBody(Astronomy.Body.Neptune,  "axis/Neptune.txt",   0.000464)   ||
        AxisTestBody(Astronomy.Body.Pluto,    "axis/Pluto.txt",     0.0)        ||
        Pass("AxisTest")
    );
}


function AxisTestBody(body, filename, arcmin_tolerance) {
    const lines = ReadLines(filename);
    let max_arcmin = 0;
    let lnum = 0;
    let count = 0;
    let found_data = false;
    for (let line of lines) {
        ++lnum;
        if (!found_data) {
            if (line === '$$SOE')
                found_data = true;
        } else {
            if (line === '$$EOE')
                break;

            const token = line.trim().split(/\s+/);
            // [ '1970-Jan-01', '00:00', '2440587.500000000', '281.01954', '61.41577' ]
            if (token.length !== 5) {
                console.error(`JS AxisBodyTest(${filename} line ${lnum}): expected 5 tokens, found ${tokens.length}`);
                return 1;
            }
            const jd = float(token[2]);
            const ra = float(token[3]);
            const dec = float(token[4]);
            const time = Astronomy.MakeTime(jd - 2451545.0);
            const axis = Astronomy.RotationAxis(body, time);
            const sphere = new Astronomy.Spherical(dec, ra, 1);
            const north = Astronomy.VectorFromSphere(sphere, time);
            const arcmin = 60 * Astronomy.AngleBetween(north, axis.north);
            if (arcmin > max_arcmin)
                max_arcmin = arcmin;

            ++count;
        }
    }
    Debug(`JS AxisTestBody(${body}): ${count} test cases, max arcmin error = ${max_arcmin}.`);
    if (max_arcmin > arcmin_tolerance) {
        console.log(`JS AxisTestBody(${body}): EXCESSIVE ERROR = ${max_arcmin} arcmin.`);
        return 1;
    }
    return 0;
}


function MoonNodes() {
    const filename = 'moon_nodes/moon_nodes.txt';
    const lines = ReadLines(filename);
    let lnum = 0;
    let prev_kind = '?';
    let max_angle = 0;
    let max_minutes = 0;
    let node;
    for (let line of lines) {
        ++lnum;
        const token = line.trim().split(/\s+/);
        if (token.length !== 4) {
            console.error(`JS MoonNodes(${filename} line ${lnum}): Unexpected number of tokens = ${tokens.length}`);
            return 1;
        }
        const kind = token[0];
        if (kind !== 'A' && kind !== 'D') {
            console.error(`JS MoonNodes(${filename} line ${lnum}): Invalid node kind '${kind}'`);
            return 1;
        }
        if (kind === prev_kind) {
            console.error(`JS MoonNodes(${filename} line ${lnum}): Duplicate node kind.`);
            return 1;
        }

        // Convert file's EQD Moon position angles to vector form.
        const time = Astronomy.MakeTime(new Date(token[1]));
        const ra = float(token[2]);
        const dec = float(token[3]);
        const sphere = new Astronomy.Spherical(dec, 15*ra, 1);
        const vec_test = Astronomy.VectorFromSphere(sphere, time);

        // Calculate EQD coordinates of the Moon. Verify against input file.
        const vec_eqj = Astronomy.GeoMoon(time);
        const rot = Astronomy.Rotation_EQJ_EQD(time);
        const vec_eqd = Astronomy.RotateVector(rot, vec_eqj);
        const angle = Astronomy.AngleBetween(vec_test, vec_eqd);
        const diff_angle = 60 * abs(angle);
        if (diff_angle > max_angle)
            max_angle = diff_angle;
        if (diff_angle > 1.54) {
            console.error(`JS MoonNodes(${filename} line ${lnum}): EXCESSIVE equatorial error = ${diff_angle.toFixed(3)} arcmin.`);
            return 1;
        }

        // Test the Astronomy Engine moon node searcher.
        if (lnum === 1) {
            // The very first time, so search for the first node in the series.
            // Back up a few days to make sure we really are finding it ourselves.
            const earlier = time.AddDays(-6.5472);      // back up by a weird amount of time
            node = Astronomy.SearchMoonNode(earlier);
        } else {
            node = Astronomy.NextMoonNode(node);
        }

        // Verify the ecliptic latitude is very close to zero at the alleged node.
        const ecl = Astronomy.EclipticGeoMoon(node.time);
        const diff_lat = 60 * abs(ecl.lat);
        if (diff_lat > 8.1e-4) {
            console.error(`JS MoonNodes(${filename} line ${lnum}): found node has excessive latitude = ${diff_lat} arcmin.`);
            return 1;
        }

        // Verify the time agrees with Espenak's time to within a few minutes.
        const diff_minutes = (24 * 60) * abs(node.time.tt - time.tt);
        if (diff_minutes > max_minutes)
            max_minutes = diff_minutes;

        // Verify the kind of node matches what Espenak says (ascending or descending).
        if (kind === 'A' && node.kind !== Astronomy.NodeEventKind.Ascending) {
            console.error(`JS MoonNodes(${filename} line ${lnum}): did not find ascending node as expected.`);
            return 1;
        }
        if (kind === 'D' && node.kind !== Astronomy.NodeEventKind.Descending) {
            console.error(`JS MoonNodes(${filename} line ${lnum}): did not find descending node as expected.`);
            return 1;
        }

        // Prepare for the next iteration.
        prev_kind = kind;
    }
    if (max_minutes > 3.681) {
        console.error(`JS MoonNodes: EXCESSIVE time prediction error = ${max_minutes.toFixed(3)} minutes.`);
        return 1;
    }
    console.log(`JS MoonNodes: PASS (${lnum} nodes, max equ error = ${max_angle.toFixed(3)} arcmin, max time error = ${max_minutes.toFixed(3)} minutes.)`);
    return 0;
}


class LagrangeFunc {
    constructor(point, major_body, minor_body) {
        this.point = point;
        this.major_body = major_body;
        this.minor_body = minor_body;
    }

    Eval(time) {
        return Astronomy.LagrangePoint(this.point, time, this.major_body, this.minor_body);
    }
}


function VerifyStateLagrange(major_body, minor_body, point, filename, r_thresh, v_thresh) {
    const func = new LagrangeFunc(point, major_body, minor_body);
    return VerifyStateBody(func, filename, r_thresh, v_thresh);
}


function LagrangeTest() {
    // Test Sun/EMB Lagrange points.
    if (0 !== VerifyStateLagrange(Astronomy.Body.Sun, Astronomy.Body.EMB, 1, "lagrange/semb_L1.txt",   1.33e-5, 6.13e-5)) return 1;
    if (0 !== VerifyStateLagrange(Astronomy.Body.Sun, Astronomy.Body.EMB, 2, "lagrange/semb_L2.txt",   1.33e-5, 6.13e-5)) return 1;
    if (0 !== VerifyStateLagrange(Astronomy.Body.Sun, Astronomy.Body.EMB, 4, "lagrange/semb_L4.txt",   3.75e-5, 5.28e-5)) return 1;
    if (0 !== VerifyStateLagrange(Astronomy.Body.Sun, Astronomy.Body.EMB, 5, "lagrange/semb_L5.txt",   3.75e-5, 5.28e-5)) return 1;

    // Test Earth/Moon Lagrange points.
    if (0 !== VerifyStateLagrange(Astronomy.Body.Earth, Astronomy.Body.Moon, 1, "lagrange/em_L1.txt",  3.79e-5, 5.06e-5)) return 1;
    if (0 !== VerifyStateLagrange(Astronomy.Body.Earth, Astronomy.Body.Moon, 2, "lagrange/em_L2.txt",  3.79e-5, 5.06e-5)) return 1;
    if (0 !== VerifyStateLagrange(Astronomy.Body.Earth, Astronomy.Body.Moon, 4, "lagrange/em_L4.txt",  3.79e-5, 1.59e-3)) return 1;
    if (0 !== VerifyStateLagrange(Astronomy.Body.Earth, Astronomy.Body.Moon, 5, "lagrange/em_L5.txt",  3.79e-5, 1.59e-3)) return 1;

    console.log("JS LagrangeTest: PASS");
    return 0;
}


function SiderealTimeTest() {
    const date = new Date('2022-03-15T21:50:00Z');
    const gast = Astronomy.SiderealTime(date);
    const correct = 9.3983699280076483;
    const diff = abs(gast - correct);
    console.log(`JS SiderealTimeTest: gast=${gast.toFixed(15)}, correct=${correct.toFixed(15)}, diff=${diff.toExponential(3)}`);
    if (diff > 1.0e-15) {
        console.error('JS SiderealTimeTest: FAIL - excessive error.');
        return 1;
    }
    console.log('JS SiderealTimeTest: PASS');
    return 0;
}

//-------------------------------------------------------------------------------------------------


function GravSimTest() {
    Debug("");

    if (0 !== GravSimEmpty("barystate/Sun.txt",      Astronomy.Body.SSB, Astronomy.Body.Sun,      0.0269, 1.9635)) return 1;
    if (0 !== GravSimEmpty("barystate/Mercury.txt",  Astronomy.Body.SSB, Astronomy.Body.Mercury,  0.5725, 0.9332)) return 1;
    if (0 !== GravSimEmpty("barystate/Venus.txt",    Astronomy.Body.SSB, Astronomy.Body.Venus,    0.1433, 0.1458)) return 1;
    if (0 !== GravSimEmpty("barystate/Earth.txt",    Astronomy.Body.SSB, Astronomy.Body.Earth,    0.0651, 0.2098)) return 1;
    if (0 !== GravSimEmpty("barystate/Mars.txt",     Astronomy.Body.SSB, Astronomy.Body.Mars,     0.1150, 0.1896)) return 1;
    if (0 !== GravSimEmpty("barystate/Jupiter.txt",  Astronomy.Body.SSB, Astronomy.Body.Jupiter,  0.2546, 0.8831)) return 1;
    if (0 !== GravSimEmpty("barystate/Saturn.txt",   Astronomy.Body.SSB, Astronomy.Body.Saturn,   0.3660, 1.0818)) return 1;
    if (0 !== GravSimEmpty("barystate/Uranus.txt",   Astronomy.Body.SSB, Astronomy.Body.Uranus,   0.3107, 0.9321)) return 1;
    if (0 !== GravSimEmpty("barystate/Neptune.txt",  Astronomy.Body.SSB, Astronomy.Body.Neptune,  0.3382, 1.5586)) return 1;

    if (0 !== GravSimEmpty("heliostate/Mercury.txt", Astronomy.Body.Sun, Astronomy.Body.Mercury,  0.5087, 0.9473)) return 1;
    if (0 !== GravSimEmpty("heliostate/Venus.txt",   Astronomy.Body.Sun, Astronomy.Body.Venus,    0.1214, 0.1543)) return 1;
    if (0 !== GravSimEmpty("heliostate/Earth.txt",   Astronomy.Body.Sun, Astronomy.Body.Earth,    0.0508, 0.2099)) return 1;
    if (0 !== GravSimEmpty("heliostate/Mars.txt",    Astronomy.Body.Sun, Astronomy.Body.Mars,     0.1085, 0.1927)) return 1;
    if (0 !== GravSimEmpty("heliostate/Jupiter.txt", Astronomy.Body.Sun, Astronomy.Body.Jupiter,  0.2564, 0.8805)) return 1;
    if (0 !== GravSimEmpty("heliostate/Saturn.txt",  Astronomy.Body.Sun, Astronomy.Body.Saturn,   0.3664, 1.0826)) return 1;
    if (0 !== GravSimEmpty("heliostate/Uranus.txt",  Astronomy.Body.Sun, Astronomy.Body.Uranus,   0.3106, 0.9322)) return 1;
    if (0 !== GravSimEmpty("heliostate/Neptune.txt", Astronomy.Body.Sun, Astronomy.Body.Neptune,  0.3381, 1.5584)) return 1;

    Debug("");
    const nsteps = 20;

    if (0 !== GravSimFile("barystate/Ceres.txt",    Astronomy.Body.SSB,   nsteps, 0.6640, 0.6226)) return 1;
    if (0 !== GravSimFile("barystate/Pallas.txt",   Astronomy.Body.SSB,   nsteps, 0.4687, 0.3474)) return 1;
    if (0 !== GravSimFile("barystate/Vesta.txt",    Astronomy.Body.SSB,   nsteps, 0.5806, 0.5462)) return 1;
    if (0 !== GravSimFile("barystate/Juno.txt",     Astronomy.Body.SSB,   nsteps, 0.6760, 0.5750)) return 1;
    if (0 !== GravSimFile("barystate/Bennu.txt",    Astronomy.Body.SSB,   nsteps, 3.7444, 2.6581)) return 1;
    if (0 !== GravSimFile("barystate/Halley.txt",   Astronomy.Body.SSB,   nsteps, 0.0539, 0.0825)) return 1;

    if (0 !== GravSimFile("heliostate/Ceres.txt",   Astronomy.Body.Sun,   nsteps, 0.0445, 0.0355)) return 1;
    if (0 !== GravSimFile("heliostate/Pallas.txt",  Astronomy.Body.Sun,   nsteps, 0.1062, 0.0854)) return 1;
    if (0 !== GravSimFile("heliostate/Vesta.txt",   Astronomy.Body.Sun,   nsteps, 0.1432, 0.1308)) return 1;
    if (0 !== GravSimFile("heliostate/Juno.txt",    Astronomy.Body.Sun,   nsteps, 0.1554, 0.1328)) return 1;

    if (0 !== GravSimFile("geostate/Ceres.txt",     Astronomy.Body.Earth, nsteps, 6.5689, 6.4797)) return 1;
    if (0 !== GravSimFile("geostate/Pallas.txt",    Astronomy.Body.Earth, nsteps, 9.3288, 7.3533)) return 1;
    if (0 !== GravSimFile("geostate/Vesta.txt",     Astronomy.Body.Earth, nsteps, 3.2980, 3.8863)) return 1;
    if (0 !== GravSimFile("geostate/Juno.txt",      Astronomy.Body.Earth, nsteps, 6.0962, 7.7147)) return 1;

    Debug("");
    console.log("JS GravSimTest: PASS");

    return 0;
}


function GravSimEmpty(filename, origin, body, rthresh, vthresh) {
    let max_rdiff = 0;
    let max_vdiff = 0;
    let sim = null;
    for (let rec of JplHorizonsStateVectors(filename)) {
        if (sim === null)
            sim = new Astronomy.GravitySimulator(origin, rec.state.t, []);

        sim.Update(rec.state.t);
        let calc = sim.SolarSystemBodyState(body);

        let rdiff = (
            (origin === Astronomy.Body.SSB && body === Astronomy.Body.Sun)
            ? SsbArcminPosError(rec.state, calc)
            : ArcminPosError(rec.state, calc)
        );

        if (rdiff > rthresh) {
            console.log(`JS GravSimEmpty(${filename} line ${rec.lnum}): excessive position error = ${rdiff} arcmin.`);
            return 1;
        }
        if (rdiff > max_rdiff)
            max_rdiff = rdiff;

        let vdiff = ArcminVelError(rec.state, calc);
        if (vdiff > vthresh) {
            console.log(`JS GravSimEmpty(${filename} line ${rec.lnum}): excessive velocity error = ${vdiff} arcmin.`);
            return 1;
        }
        if (vdiff > max_vdiff)
            max_vdiff = vdiff;
    }
    Debug(`JS GravSimFile(${filename}): PASS - max pos error = ${max_rdiff.toFixed(4)} arcmin, max vel error = ${max_vdiff.toFixed(4)} arcmin.`);
    return 0;
}

function GravSimFile(filename, originBody, nsteps, rthresh, vthresh) {
    let sim = null;
    let prev = null;
    let time = null;
    let max_rdiff = 0;
    let max_vdiff = 0;
    let smallBodyArray;
    for (let rec of JplHorizonsStateVectors(filename)) {
        if (sim === null) {
            sim = new Astronomy.GravitySimulator(originBody, rec.state.t, [rec.state]);
            time = rec.state.t;
            smallBodyArray = sim.Update(time);
        } else {
            let tt1 = prev.state.t.tt;
            let tt2 = rec.state.t.tt;
            let dt = (tt2 - tt1) / nsteps;
            for (let k = 1; k <= nsteps; ++k) {
                time = Astronomy.AstroTime.FromTerrestrialTime(tt1 + k*dt);
                smallBodyArray = sim.Update(time);
                if (smallBodyArray.length !== 1) {
                    console.log(`JS GravSimFile(${filename} line ${rec.lnum}): unexpected smallBodyArray.length = ${smallBodyArray.length}.`);
                    return 1;
                }
                if (time.tt !== sim.Time.tt) {
                    console.log(`JS GravSimFile(${filename} line ${rec.lnum}): expected ${time} but simulator reports ${sim.Time}.`);
                    return 1;
                }
            }
        }

        const rdiff = ArcminPosError(rec.state, smallBodyArray[0]);
        if (rdiff > rthresh) {
            console.log(`JS GravSimFile(${filename} line ${rec.lnum}): excessive position error = ${rdiff} arcmin.`);
            return 1;
        }
        if (rdiff > max_rdiff)
            max_rdiff = rdiff;

        const vdiff = ArcminVelError(rec.state, smallBodyArray[0]);
        if (vdiff > vthresh) {
            console.log(`JS GravSimFile(${filename} line ${rec.lnum}): excessive velocity error = ${vdiff} arcmin.`);
            return 1;
        }
        if (vdiff > max_vdiff)
            max_vdiff = vdiff;
        prev = rec;
    }
    Debug(`JS GravSimFile(${filename}): PASS - max pos error = ${max_rdiff.toFixed(4)} arcmin, max vel error = ${max_vdiff.toFixed(4)} arcmin.`);
    return 0;
}


function SsbArcminPosError(correct, calc) {
    // Scale the SSB based on 1 AU, not on its absolute magnitude, which can become very close to zero.
    const dx = calc.x - correct.x;
    const dy = calc.y - correct.y;
    const dz = calc.z - correct.z;
    const diffSquared = dx*dx + dy*dy + dz*dz;
    const radians = sqrt(diffSquared);
    return (60.0 * Astronomy.RAD2DEG) * radians;
}

function ArcminPosError(correct, calc) {
    const dx = calc.x - correct.x;
    const dy = calc.y - correct.y;
    const dz = calc.z - correct.z;
    const diffSquared = dx*dx + dy*dy + dz*dz;
    const magSquared = correct.x*correct.x + correct.y*correct.y + correct.z*correct.z;
    const radians = sqrt(diffSquared / magSquared);
    return (60.0 * Astronomy.RAD2DEG) * radians;
}

function ArcminVelError(correct, calc) {
    const dx = calc.vx - correct.vx;
    const dy = calc.vy - correct.vy;
    const dz = calc.vz - correct.vz;
    const diffSquared = dx*dx + dy*dy + dz*dz;
    const magSquared = correct.vx*correct.vx + correct.vy*correct.vy + correct.vz*correct.vz;
    const radians = sqrt(diffSquared / magSquared);
    return (60.0 * Astronomy.RAD2DEG) * radians;
}


//-------------------------------------------------------------------------------------------------

function CheckDecemberSolstice(year, expected) {
    const si = Astronomy.Seasons(year);
    const actual = si.dec_solstice.toString();
    if (actual != expected) {
        console.error(`JS DatesIssue250 FAIL: year ${year}, expected [${expected}], actual [${actual}]`);
        return 1;
    }
    return 0;
}

function DatesIssue250() {
    return (
        CheckDecemberSolstice( 2022, "2022-12-21T21:47:54.455Z") ||
        CheckDecemberSolstice(-2300, "-002300-12-19T16:22:27.930Z") ||
        CheckDecemberSolstice(12345, "+012345-12-11T13:30:10.275Z") ||
        Pass('DatesIssue250')
    );
}

//-------------------------------------------------------------------------------------------------

function LunarFractionCase(dateText, obscuration) {
    // Search for the first lunar eclipse to occur after the given date.
    const time = Astronomy.MakeTime(new Date(dateText));
    const eclipse = Astronomy.SearchLunarEclipse(time);

    // This should be a partial lunar eclipse.
    if (eclipse.kind != Astronomy.EclipseKind.Partial) {
        console.error(`JS LunarFractionCase(${dateText}) FAIL: expected partial eclipse, found ${eclipse.kind}.`);
        return 1;
    }

    // The partial eclipse should always happen within 24 hours of the given date.
    const dt = v(eclipse.peak.ut - time.ut);
    if (dt < 0.0 || dt > 1.0) {
        console.error(`JS LunarFractionCase(${dateText}) FAIL: eclipse occurs ${dt.toFixed(4)} days after search date.`);
        return 1;
    }

    // Check accuracy of the calculated obscuration.
    const diff = v(eclipse.obscuration - obscuration);
    if (abs(diff) > 0.00763) {
        console.error(`JS LunarFractionCase(${dateText}) FAIL: obscuration error = ${diff}, expected = ${obscuration}, actual = ${eclipse.obscuration}`);
        return 1;
    }

    Debug(`JS LunarFractionCase(${dateText}) obscuration error = ${diff.toFixed(8)}`);
    return 0;
}

function LunarFraction() {
    // Verify calculation of the fraction of the Moon's disc covered by the Earth's umbra during a partial eclipse.
    // Data for this is more tedious to gather, because Espenak data does not contain it.
    // We already verify fraction=0.0 for penumbral eclipses and fraction=1.0 for total eclipses in LunarEclipseTest.
    return (
        LunarFractionCase('2010-06-26', 0.506) ||  // https://www.timeanddate.com/eclipse/lunar/2010-june-26
        LunarFractionCase('2012-06-04', 0.304) ||  // https://www.timeanddate.com/eclipse/lunar/2012-june-4
        LunarFractionCase('2013-04-25', 0.003) ||  // https://www.timeanddate.com/eclipse/lunar/2013-april-25
        LunarFractionCase('2017-08-07', 0.169) ||  // https://www.timeanddate.com/eclipse/lunar/2017-august-7
        LunarFractionCase('2019-07-16', 0.654) ||  // https://www.timeanddate.com/eclipse/lunar/2019-july-16
        LunarFractionCase('2021-11-19', 0.991) ||  // https://www.timeanddate.com/eclipse/lunar/2021-november-19
        LunarFractionCase('2023-10-28', 0.060) ||  // https://www.timeanddate.com/eclipse/lunar/2023-october-28
        LunarFractionCase('2024-09-18', 0.035) ||  // https://www.timeanddate.com/eclipse/lunar/2024-september-18
        LunarFractionCase('2026-08-28', 0.962) ||  // https://www.timeanddate.com/eclipse/lunar/2026-august-28
        LunarFractionCase('2028-01-12', 0.024) ||  // https://www.timeanddate.com/eclipse/lunar/2028-january-12
        LunarFractionCase('2028-07-06', 0.325) ||  // https://www.timeanddate.com/eclipse/lunar/2028-july-6
        LunarFractionCase('2030-06-15', 0.464) ||  // https://www.timeanddate.com/eclipse/lunar/2030-june-15
        Pass("LunarFraction")
    );
}

//-------------------------------------------------------------------------------------------------

function StarRiseSetCulmCase(starName, ra, dec, distLy, observer, searchDateText, riseText, culmText, setText) {
    Astronomy.DefineStar(Astronomy.Body.Star1, ra, dec, distLy);
    const searchTime = new Astronomy.AstroTime(new Date(searchDateText));
    const riseTime = Astronomy.SearchRiseSet(Astronomy.Body.Star1, observer, +1, searchTime, 1.0) || Fail('Star rise search failed.');
    const setTime  = Astronomy.SearchRiseSet(Astronomy.Body.Star1, observer, -1, searchTime, 1.0) || Fail('Star set search failed.');
    const culm = Astronomy.SearchHourAngle(Astronomy.Body.Star1, observer, 0.0, searchTime) || Fail('Star culmination search failed.');
    const expectedRiseTime = new Astronomy.AstroTime(new Date(searchDateText + 'T' + riseText + 'Z'));
    const expectedCulmTime = new Astronomy.AstroTime(new Date(searchDateText + 'T' + culmText + 'Z'));
    const expectedSetTime = new Astronomy.AstroTime(new Date(searchDateText + 'T' + setText + 'Z'));
    const rdiff = MINUTES_PER_DAY * abs(expectedRiseTime.ut - riseTime.ut);
    const cdiff = MINUTES_PER_DAY * abs(expectedCulmTime.ut - culm.time.ut);
    const sdiff = MINUTES_PER_DAY * abs(expectedSetTime.ut - setTime.ut);
    Debug(`StarRiseSetCulmCase(${starName}): rise=${riseTime} (err=${rdiff.toFixed(4)} min), culm=${culm.time} (err=${cdiff.toFixed(4)} min), set=${setTime} (err=${sdiff.toFixed(4)} min).`);
    if (rdiff > 0.5) Fail(`StarRiseSetCulmCase(${starName}): excessive rise time error = ${rdiff} minutes.`);
    if (cdiff > 0.5) Fail(`StarRiseSetCulmCase(${starName}): excessive culm time error = ${cdiff} minutes.`);
    if (sdiff > 0.5) Fail(`StarRiseSetCulmCase(${starName}): excessive set time error = ${sdiff} minutes.`);
    return 0;
}

function StarRiseSetCulm() {
    // https://aa.usno.navy.mil/calculated/mrst?body=-20&date=2022-11-21&reps=5&lat=25.77&lon=-80.19&label=&tz=0.00&tz_sign=-1&height=0&submit=Get+Data
    //       Date               Rise  Az.       Transit Alt.       Set  Az.
    // 2022 Nov 21 (Mon)        02:37 108        08:06 47S        13:34 252

    const observer = new Astronomy.Observer(+25.77, -80.19, 0.0);
    return (
        StarRiseSetCulmCase('Sirius', 6.7525, -16.7183, 8.6, observer, '2022-11-21', '02:37', '08:06', '13:34') ||
        StarRiseSetCulmCase('Sirius', 6.7525, -16.7183, 8.6, observer, '2022-11-25', '02:22', '07:50', '13:18') ||
        StarRiseSetCulmCase('Canopus', 6.3992, -52.6956, 310, observer, '2022-11-21', '04:17', '07:44', '11:11') ||
        StarRiseSetCulmCase('Canopus', 6.3992, -52.6956, 310, observer, '2022-11-25', '04:01', '07:28', '10:56') ||
        Pass('StarRiseSetCulm')
    );
}

//-------------------------------------------------------------------------------------------------

function EclipticTest() {
    let time = Astronomy.MakeTime(new Date('1900-01-01T00:00:00Z'));
    const stopTime = Astronomy.MakeTime(new Date('2100-01-01T00:00:00Z'));
    let max_vec_diff = 0.0;
    let max_angle_diff = 0.0;
    let count = 0;
    while (time.ut <= stopTime.ut) {
        // Get Moon's geocentric position in EQJ.
        const eqj = Astronomy.GeoMoon(time);

        // Convert EQJ to ECT.
        const eclip = Astronomy.Ecliptic(eqj);

        // Confirm that the ecliptic angles and ecliptic vectors are consistent.
        const check_sphere = new Astronomy.Spherical(eclip.elat, eclip.elon, eclip.vec.Length());
        const check_vec = Astronomy.VectorFromSphere(check_sphere, time);
        max_angle_diff = Math.max(max_angle_diff, VectorDiff(eclip.vec, check_vec));

        // Independently get the Moon's spherical coordinates in ECT.
        const sphere = Astronomy.EclipticGeoMoon(time);

        // Convert spherical coordinates to ECT vector.
        const check_ect = Astronomy.VectorFromSphere(sphere, time);

        // Verify the two ECT vectors are identical, within tolerance.
        max_vec_diff = Math.max(max_vec_diff, VectorDiff(eclip.vec, check_ect));

        time = time.AddDays(10.0);
        ++count;
    }

    if (max_angle_diff > 3.040e-18)
        throw `EclipticTest: EXCESSIVE ANGLE DIFF = ${max_angle_diff.toExponential(6)} AU.`;

    if (max_vec_diff > 3.743e-18)
        throw `EclipticTest: EXCESSIVE VECTOR DIFF = ${max_vec_diff.toExponential(6)} AU.`;

    console.log(`JS EclipticTest: PASS: count = ${count}, max_vec_diff = ${max_vec_diff.toExponential(6)} AU., max_angle_diff = ${max_angle_diff.toExponential(6)} AU.`);
    return 0;
}

//-------------------------------------------------------------------------------------------------

function HourAngleCase(latitude, longitude, hourAngle, context) {
    const threshold = 0.1 / 3600;   // SearchHourAngle() accuracy: 0.1 seconds converted to hours
    const observer = new Astronomy.Observer(latitude, longitude, 0);
    const startTime = new Astronomy.AstroTime(new Date('2023-02-11T00:00:00Z'));
    const search = Astronomy.SearchHourAngle(Astronomy.Body.Sun, observer, hourAngle, startTime, +1);
    const calc = Astronomy.HourAngle(Astronomy.Body.Sun, search.time, observer);
    let diff = abs(calc - hourAngle);
    if (diff > 12.0)
        diff = 24.0 - diff;
    if (diff > context.maxdiff)
        context.maxdiff = diff;
    ++context.cases;
    if (diff > threshold) {
        console.error(`JS HourAngleCase: EXCESSIVE ERROR = ${diff.toExponential(6)}, calc HA = ${calc.toFixed(16)}, for hourAngle=${hourAngle.toFixed(1)}`);
        return false;
    }
    Debug(`JS HourAngleCase: Hour angle = ${hourAngle.toFixed(1)}, longitude = ${longitude.toFixed(1)}, diff = ${diff.toExponential(4)}`);
    return true;
}

function HourAngleTest() {
    const context = {
        cases: 0,
        maxdiff: 0.0
    };

    const latitude = 35.0;
    for (let longitude = -170; longitude <= 180; longitude += 5) {
        for (let hour = 0; hour < 24; ++hour) {
            if (!HourAngleCase(latitude, longitude, hour, context)) {
                return 1;
            }
        }
    }

    return Pass(`HourAngleTest (${context.cases} cases, maxdiff = ${context.maxdiff.toExponential(4)})`);
}

//-------------------------------------------------------------------------------------------------

function Atmosphere() {
    const filename = 'riseset/atmosphere.csv';
    const lines = ReadLines(filename);
    const tolerance = 8.8e-11;
    let lnum = 0;
    let ncases = 0;
    let maxdiff = 0.0;
    for (let line of lines) {
        ++lnum;
        if (lnum === 1) {
            if (line != "elevation,temperature,pressure,density,relative_density")
                Fail(`Atmosphere: Expected header line, but found: [${line}]`);
        } else {
            const tokens = line.split(',');
            if (tokens.length !== 5)
                Fail(`Atmosphere(${filename} line ${line}): expected 5 numeric tokens but found ${tokens.length}`);

            const elevation = float(tokens[0]);
            const temperature = float(tokens[1]);
            const pressure = float(tokens[2]);
            // ignore tokens[3] = absolute_density
            const relative_density = float(tokens[4]);

            const atmos = Astronomy.Atmosphere(elevation);

            let diff = abs(atmos.temperature - temperature);
            if (diff > maxdiff) maxdiff = diff;
            if (diff > tolerance) Fail(`Atmosphere(${filename} line ${lnum}): EXCESSIVE temperature difference = ${diff}`);

            diff = abs(atmos.pressure - pressure);
            if (diff > maxdiff) maxdiff = diff;
            if (diff > tolerance) Fail(`Atmosphere(${filename} line ${lnum}): EXCESSIVE pressure difference = ${diff}`);

            diff = abs(atmos.density - relative_density);
            if (diff > maxdiff) maxdiff = diff;
            if (diff > tolerance) Fail(`Atmosphere(${filename} line ${lnum}): EXCESSIVE density difference = ${diff}`);

            ++ncases;
        }
    }
    if (ncases !== 34)
        Fail(`Atmosphere: expected 34 test cases but found ${ncases}.`);
    return Pass('Atmosphere');
}

//-------------------------------------------------------------------------------------------------

function RiseSetElevationBodyCase(body, observer, direction, metersAboveGround, startTime, eventOffsetDays) {
    const time = Astronomy.SearchRiseSet(body, observer, direction, startTime, 2.0, metersAboveGround);
    if (!time) Fail(`Failed to find event for ${body} direction=${direction}`);
    let diff = v(time.ut - (startTime.ut + eventOffsetDays));
    if (diff > 0.5)
        diff -= 1.0;    // assume event actually takes place on the next day

    diff = abs(MINUTES_PER_DAY * diff);
    if (diff > 0.5)
        Fail(`RiseSetElevationBodyCase: body=${body}, direction=${direction}, EXCESSIVE diff = ${diff} minutes.`);
}


function RiseSetElevation() {
    const re = /^(\d+)-(\d+)-(\d+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\d+):(\d+)\s+(\d+):(\d+)\s+(\d+):(\d+)\s+(\d+):(\d+)\s+(\S+)\s*$/;
    const filename = 'riseset/elevation.txt';
    const lines = ReadLines(filename);
    let lnum = 0;
    for (let line of lines) {
        ++lnum;
        if (line.startsWith('#')) continue;
        const m = line.match(re);
        if (!m) Fail(`RiseSetElevation(${filename} line ${lnum}): invalid data format.`);
        const year = int(m[1]);
        const month = int(m[2]);
        const day = int(m[3]);
        const latitude = float(m[4]);
        const longitude = float(m[5]);
        const height = float(m[6]);
        const metersAboveGround = float(m[7]);
        const srh = int(m[ 8]);
        const srm = int(m[ 9]);
        const ssh = int(m[10]);
        const ssm = int(m[11]);
        const mrh = int(m[12]);
        const mrm = int(m[13]);
        const msh = int(m[14]);
        const msm = int(m[15]);

        // Get search time origin.
        const time = Astronomy.MakeTime(new Date(Date.UTC(year, month-1, day, 0, 0)));

        // Convert scanned values into sunrise, sunset, moonrise, moonset day offsets.
        const sr = (srh + srm/60)/24;
        const ss = (ssh + ssm/60)/24;
        const mr = (mrh + mrm/60)/24;
        const ms = (msh + msm/60)/24;

        const observer = new Astronomy.Observer(latitude, longitude, height);

        RiseSetElevationBodyCase(Astronomy.Body.Sun,  observer, +1, metersAboveGround, time, sr);
        RiseSetElevationBodyCase(Astronomy.Body.Sun,  observer, -1, metersAboveGround, time, ss);
        RiseSetElevationBodyCase(Astronomy.Body.Moon, observer, +1, metersAboveGround, time, mr);
        RiseSetElevationBodyCase(Astronomy.Body.Moon, observer, -1, metersAboveGround, time, ms);
    }
    return Pass('RiseSetElevation');
}

//-------------------------------------------------------------------------------------------------

const UnitTests = {
    aberration:             AberrationTest,
    atmosphere:             Atmosphere,
    axis:                   AxisTest,
    barystate:              BaryStateTest,
    constellation:          Constellation,
    dates250:               DatesIssue250,
    ecliptic:               EclipticTest,
    elongation:             Elongation,
    geoid:                  Geoid,
    global_solar_eclipse:   GlobalSolarEclipse,
    gravsim:                GravSimTest,
    heliostate:             HelioStateTest,
    hour_angle:             HourAngleTest,
    issue_103:              Issue103,
    jupiter_moons:          JupiterMoons,
    lagrange:               LagrangeTest,
    libration:              LibrationTest,
    local_solar_eclipse:    LocalSolarEclipse,
    lunar_apsis:            LunarApsis,
    lunar_eclipse:          LunarEclipse,
    lunar_eclipse_78:       LunarEclipseIssue78,
    lunar_fraction:         LunarFraction,
    magnitude:              Magnitude,
    moon_nodes:             MoonNodes,
    moon_phase:             MoonPhase,
    moon_reverse:           MoonReverse,
    planet_apsis:           PlanetApsis,
    pluto:                  PlutoCheck,
    refraction:             Refraction,
    riseset:                RiseSet,
    riseset_reverse:        RiseSetReverse,
    riseset_elevation:      RiseSetElevation,
    rotation:               Rotation,
    seasons:                Seasons,
    seasons187:             SeasonsIssue187,
    sidereal:               SiderealTimeTest,
    solar_fraction:         SolarFraction,
    star_risesetculm:       StarRiseSetCulm,
    topostate:              TopoStateTest,
    transit:                Transit,
    twilight:               TwilightTest,
    vecobs:                 VectorObserverIssue347,
};


function TestAll() {
    for (let name in UnitTests) {
        if (UnitTests[name]()) {
            return 1;
        }
    }
    console.log('JS ALL PASS');
}


function main() {
    let args = process.argv.slice(2);
    if (args.length > 0 && args[0] === '-v') {
        Verbose = true;
        args = args.slice(1);
    }

    if (args.length === 1) {
        const name = args[0];
        if (name === 'all') {
            return TestAll();
        }
        if (name === 'astro_check') {
            // This is a special case because the output is redirected and parsed.
            // It needs to be invoked separately, and cannot emit any extraneous output.
            return AstroCheck();
        }
        const func = UnitTests[name];
        if (typeof func !== 'function') {
            console.log(`test.js: Unknown unit test "${name}"`);
            return 1;
        }
        return func();
    }

    console.log('test.js: Invalid command line arguments.');
    return 1;
}


process.exit(main());