satellite-brightness-eclipse.js

// by cmdr2, anyone is free to use and modify this, without the need to give credit

const satellite = require('satellite.js'); // https://www.npmjs.com/package/satellite.js
const SunCalc = require('suncalc'); // https://www.npmjs.com/package/suncalc

var satrec = satellite.twoline2satrec('1 71066C 20001A   20014.19313215 -.01059873  00000-0 -43936-2 0   145',
				      '2 71066  52.9993   1.7743 0000976  92.9471 326.1492 15.86148217    16');

var stdMag = 4.7; // for Starlink-2

var SUN_TO_EARTH_DIST = 147124525.068; // Km
var EARTH_RADIUS = 6378.16; // Km
var SUN_RADIUS = 695510; // Km

var observerGd = {
    latitude: satellite.degreesToRadians(12.2),
    longitude: satellite.degreesToRadians(-30.2),
    height: 0
};

(function test() {
    var date = new Date();
    var positionAndVelocity = satellite.propagate(satrec, date);
    var positionEci = positionAndVelocity.position;

    var gmst = satellite.gstime(date);

    var positionEcf = satellite.eciToEcf(positionEci, gmst);
    var lookAngles = satellite.ecfToLookAngles(observerGd, positionEcf);
    lookAngles.distance = lookAngles.rangeSat;

    var sunXYZ = getSunPosition(date);

    var eclipsed = isEclipsed(positionEci, sunXYZ);
    var brightness = getBrightness(date, lookAngles, stdMag);

    console.log(date, lookAngles, eclipsed, brightness)
})();


function isEclipsed(satXYZ, sunXYZ) {
    var a = mag(sunXYZ);
    var b = mag(satXYZ);
    var sat_to_sun = diff(sunXYZ, satXYZ);
    var c = mag(sat_to_sun);

    var phase_angle = Math.acos((Math.pow(b, 2) + Math.pow(c, 2) - Math.pow(a, 2)) / (2 * b * c));
    var R_e = EARTH_RADIUS
    var R_s = SUN_RADIUS
    var pE = b
    var pS = c

    var theta_e = Math.asin(R_e / pE)
    var theta_s = Math.asin(R_s / pS)

    var full_eclipse = ((theta_e > theta_s) && (phase_angle < (theta_e - theta_s)))
    var partial_eclipse = (Math.abs(theta_e - theta_s) < phase_angle) && (phase_angle < (theta_e + theta_s))

    return full_eclipse || partial_eclipse
}

function getBrightness(date, lookAngles, stdMag) {
    var sunPos = SunCalc.getPosition(date, observer.latitude, observer.longitude);
    sunPos.azimuth += Math.PI; // SunCalc's 0 is south, satellite-js's 0 is north
    sunPos.elevation = sunPos.altitude;
    // sunPos.altitude = satellite.degreesToRadians(sunPos.altitude);
    // sunPos.azimuth = satellite.degreesToRadians(sunPos.azimuth);

    var cosDelta = Math.sin(lookAngles.elevation) * Math.sin(sunPos.elevation)
                   + Math.cos(lookAngles.elevation) * Math.cos(sunPos.elevation) * Math.cos(lookAngles.azimuth - sunPos.azimuth);
    var angleC = Math.acos(cosDelta);

    var a = SUN_TO_EARTH_DIST - EARTH_RADIUS - SUN_RADIUS; // # Km
    var b = lookAngles.distance;
    var c = Math.sqrt(Math.pow(a, 2) + Math.pow(b, 2) - 2*a*b*Math.cos(angleC));

    var phaseAngle = Math.acos((Math.pow(b, 2) + Math.pow(c, 2) - Math.pow(a, 2)) / (2 * b * c));

    // calc magnitude
    var term1 = stdMag;
    var term2 = 5.0 * Math.log10(b / 1000);

    var arg = Math.sin(phaseAngle) + (Math.PI - phaseAngle) * Math.cos(phaseAngle);
    var term3 = -2.5 * Math.log10(arg);

    var apparentMag = term1 + term2 + term3;

    return [apparentMag, satellite.radiansToDegrees(sunPos.azimuth), satellite.radiansToDegrees(sunPos.altitude)];
}

function getSunPosition(date) {
    var rad = Math.PI / 180;
 //    var time = new Orb.Time(date)
    // var jd = time.jd();// + dt;
    var jd = satellite.jday(date);
    var t = (jd -2451545.0)/36525;
    var mean_longitude = 280.46646 + 36000.76983*t + 0.0003032*t*t;
    var mean_anomaly =  357.52911+ 35999.05029*t - 0.0001537*t*t;
    var eccentricity = 0.016708634 - 0.000042037*t - 0.0000001267*t*t;
    var equation = (1.914602 - 0.004817*t - 0.000014*t*t)*Math.sin(mean_anomaly*rad);
    equation += (0.019993 - 0.000101*t)*Math.sin(2*mean_anomaly*rad);
    equation += 0.000289 *Math.sin(3*mean_anomaly*rad);
    var true_longitude = mean_longitude + equation;
    var true_anomary = mean_anomaly + equation;
    var radius = (1.000001018*(1-eccentricity*eccentricity))/(1 + eccentricity*Math.cos(true_anomary*rad));
    var nao = new NutationAndObliquity(date)
    var nutation = nao.nutation();
    var obliquity = nao.obliquity();
    var apparent_longitude = true_longitude + nutation;
    var longitude = apparent_longitude;
    var distance=radius*149597870.691;

    var x = distance*Math.cos(longitude*rad);
    var y = distance*(Math.sin(longitude*rad)*Math.cos(obliquity*rad));
    var z = distance*(Math.sin(longitude*rad)*Math.sin(obliquity*rad));

    return {x: x, y: y, z: z};
}

function NutationAndObliquity(date) {
    var rad = Math.PI / 180;
    var jd = satellite.jday(date)
    var t = (jd -2451545.0)/36525;
    var omega = (125.04452 - 1934.136261*t+0.0020708*t*t + (t*t+t)/450000)*rad;
    var L0 = (280.4665 + 36000.7698*t)*rad
    var L1 = (218.3165 + 481267.8813*t)*rad
    return {
        nutation:function(){
            var nutation = (-17.20/3600)*Math.sin(omega)-(-1.32/3600)*Math.sin(2*L0)-(0.23/3600)*Math.sin(2*L1)+(0.21/3600)*Math.sin(2*omega)/rad;
            return nutation;
        },
        obliquity:function(){
            var obliquity_zero = 23+26.0/60+21.448/3600 -(46.8150/3600)*t -(0.00059/3600)*t*t +(0.001813/3600)*t*t*t;
            var obliquity_delta = (9.20/3600)*Math.cos(omega) + (0.57/3600)*Math.cos(2*L0) +(0.10/3600)*Math.cos(2*L1)-(0.09/3600)*Math.cos(2*omega);
            var obliquity= obliquity_zero + obliquity_delta;
            return obliquity;
        }
    }
}

function mag(v) {
    return Math.sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
}

function diff(v0, v1) {
    return {x: v0.x - v1.x, y: v0.y - v1.y, z: v0.z - v1.z};
}