mirror of
https://github.com/eddyem/bta-panels.git
synced 2025-12-06 02:35:17 +03:00
306 lines
9.6 KiB
JavaScript
306 lines
9.6 KiB
JavaScript
// http://www.esrl.noaa.gov/gmd/grad/solcalc/azel.html
|
|
(function () { 'use strict';
|
|
var NOAAcalc = {};
|
|
const
|
|
PI = Math.PI,
|
|
PI_2 = PI/2,
|
|
_2PI = 2*PI,
|
|
DegToRad = PI/180.,
|
|
RadToDeg = 180./PI,
|
|
J2000 = 2451545.0; // JD of 2000yr
|
|
var
|
|
floor = Math.floor,
|
|
abs = Math.abs,
|
|
sqrt = Math.sqrt;
|
|
// all inner angle functions are in degrees!
|
|
function sin(deg){return Math.sin(deg * DegToRad);}
|
|
function cos(deg){return Math.cos(deg * DegToRad);}
|
|
function tan(deg){return Math.tan(deg * DegToRad);}
|
|
function asin(deg){return RadToDeg * Math.asin(deg);}
|
|
function acos(deg){return RadToDeg * Math.acos(deg);}
|
|
function atan(y,x){return RadToDeg * Math.atan2(y,x);}
|
|
function d360(angle){
|
|
return(angle - floor(angle/360.0)*360);
|
|
}
|
|
function d180(angle){
|
|
var a = d360(angle);
|
|
if(a > 180) a -= 360;
|
|
return a;
|
|
}
|
|
|
|
function calcJD(d){
|
|
var dayMs = 1000 * 60 * 60 * 24,
|
|
J1970 = 2440588;
|
|
return (d.valueOf()/dayMs - 0.5 + J1970)
|
|
}
|
|
NOAAcalc.calcJD = calcJD;
|
|
function calcTimeJulianCent(jd){
|
|
return ((jd - J2000)/36525.0);
|
|
}
|
|
// argument t in all functions is in Julian cents
|
|
function calcGeomMeanAnomalySun(t){
|
|
return (357.52911 + t * (35999.05029 - 0.0001537 * t));
|
|
}
|
|
function calcSunEqOfCenter(t){
|
|
var m = calcGeomMeanAnomalySun(t);
|
|
var sinm = sin(m);
|
|
var sin2m = sin(2*m);
|
|
var sin3m = sin(3*m);
|
|
return (sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) +
|
|
sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289);
|
|
}
|
|
function calcSunTrueAnomaly(t){
|
|
var m = calcGeomMeanAnomalySun(t);
|
|
var c = calcSunEqOfCenter(t);
|
|
return d360(m + c);
|
|
}
|
|
function calcEccentricityEarthOrbit(t){
|
|
return (0.016708634 - t * (0.000042037 + 0.0000001267 * t));
|
|
}
|
|
function calcSunRadVector(t){
|
|
var v = calcSunTrueAnomaly(t);
|
|
var e = calcEccentricityEarthOrbit(t);
|
|
var R = (1.000001018 * (1 - e * e)) / (1 + e * cos(v));
|
|
return R;
|
|
}
|
|
function calcMeanObliquityOfEcliptic(t){
|
|
var seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)));
|
|
return (23.0 + (26.0 + (seconds/60.0))/60.0);
|
|
}
|
|
function calcObliquityCorrection(t){
|
|
var e0 = calcMeanObliquityOfEcliptic(t);
|
|
var omega = 125.04 - 1934.136 * t;
|
|
var e = e0 + 0.00256 * cos(omega);
|
|
return d360(e);
|
|
}
|
|
function calcGeomMeanLongSun(t){
|
|
var L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
|
|
return d360(L0);
|
|
}
|
|
function calcSunTrueLong(t){
|
|
var l0 = calcGeomMeanLongSun(t);
|
|
var c = calcSunEqOfCenter(t);
|
|
return d360(l0 + c);
|
|
}
|
|
function calcSunApparentLong(t){
|
|
var o = calcSunTrueLong(t);
|
|
var omega = 125.04 - 1934.136 * t;
|
|
return d360(o - (0.00569 - 0.00478 * sin(omega)));
|
|
}
|
|
function calcSunRtAscension(t){
|
|
var e = calcObliquityCorrection(t);
|
|
var lambda = calcSunApparentLong(t);
|
|
var tananum = cos(e) * sin(lambda);
|
|
var tanadenom = cos(lambda);
|
|
return atan(tananum, tanadenom);
|
|
}
|
|
function calcSunDeclination(t){
|
|
var e = calcObliquityCorrection(t);
|
|
var lambda = calcSunApparentLong(t);
|
|
var sint = sin(e) * sin(lambda);
|
|
return asin(sint);
|
|
}
|
|
function calcEquationOfTime(t){
|
|
var epsilon = calcObliquityCorrection(t);
|
|
var l0 = calcGeomMeanLongSun(t);
|
|
var e = calcEccentricityEarthOrbit(t);
|
|
var m = calcGeomMeanAnomalySun(t);
|
|
var y = tan(epsilon/2.0);
|
|
y *= y;
|
|
var sin2l0 = sin(2.0 * l0);
|
|
var sinm = sin(m);
|
|
var cos2l0 = cos(2.0 * l0);
|
|
var sin4l0 = sin(4.0 * l0);
|
|
var sin2m = sin(2.0 * m);
|
|
var Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
|
|
- 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;
|
|
return Etime*4.0; // in minutes of time
|
|
}
|
|
function Refract(zenith){
|
|
var refractionCorrection = 0., exoatmElevation = 90. - zenith;
|
|
if(exoatmElevation < 85){
|
|
var te = tan(exoatmElevation);
|
|
if (exoatmElevation > 5){
|
|
refractionCorrection = 58.1 / te - 0.07 / (te*te*te) +
|
|
0.000086 / (te*te*te*te*te);
|
|
} else if (exoatmElevation > -0.575) {
|
|
refractionCorrection = 1735.0 + exoatmElevation *
|
|
(-518.2 + exoatmElevation * (103.4 +
|
|
exoatmElevation * (-12.79 + exoatmElevation * 0.711) ) );
|
|
} else {
|
|
refractionCorrection = -20.774 / te;
|
|
}
|
|
refractionCorrection /= 3600.0;
|
|
}
|
|
return refractionCorrection;
|
|
}
|
|
NOAAcalc.getSunPosition = function(adate, latitude, longitude){
|
|
var JD = calcJD(adate);
|
|
var T = calcTimeJulianCent(JD);
|
|
//var R = calcSunRadVector(T);
|
|
var alpha = calcSunRtAscension(T);
|
|
var solarDec = calcSunDeclination(T);
|
|
var eqTime = calcEquationOfTime(T);
|
|
|
|
var solarTimeFix = eqTime + 4.*longitude + adate.getTimezoneOffset();
|
|
var trueSolarTime = adate.getHours() * 60. + adate.getMinutes() +
|
|
adate.getSeconds()/60. + solarTimeFix;
|
|
// in minutes
|
|
trueSolarTime -= floor(trueSolarTime/1440)*1440;
|
|
|
|
var hourAngle = d360(trueSolarTime / 4.0 - 180.0);
|
|
|
|
var csz = sin(latitude) * sin(solarDec) + cos(latitude) * cos(solarDec) * cos(hourAngle);
|
|
if(csz > 1.0){csz = 1.0;} else if(csz < -1.0){csz = -1.0;}
|
|
var azimuth, zenith = acos(csz);
|
|
var azDenom = cos(latitude) * sin(zenith);
|
|
if(abs(azDenom) > 0.001){
|
|
var azRad = (sin(latitude) * cos(zenith) - sin(solarDec)) / azDenom;
|
|
if(abs(azRad) > 1.0){
|
|
if(azRad < 0){azRad = -1.0;} else {azRad = 1.0;}
|
|
}
|
|
azimuth = acos(azRad);
|
|
}else{
|
|
if(latitude > 0.0){azimuth = 0.;}else{azimuth = 180.;}
|
|
}
|
|
|
|
var solarZen = zenith - Refract(zenith);
|
|
return{
|
|
azimuth: azimuth,
|
|
altitude: 90. - solarZen,
|
|
solarTimeFix: solarTimeFix,
|
|
RA: alpha,
|
|
Dec: solarDec
|
|
};
|
|
};
|
|
|
|
function daynumber(adate){
|
|
var Y = adate.getUTCFullYear(), M = adate.getUTCMonth() + 1;
|
|
var d = 367 * Y - floor(7*(Y + floor((M+9)/12))/4) + floor(275*M/9) + adate.getUTCDate() - 730530;
|
|
d += (adate.getUTCHours() + adate.getUTCMinutes()/60. + adate.getUTCSeconds()/3600.)/24.;
|
|
return d;
|
|
}
|
|
// http://www.stjarnhimlen.se/comp/ppcomp.html
|
|
NOAAcalc.getMoonPosition = function(adate, latitude, longitude){
|
|
/*
|
|
Primary orbital elements of the Moon:
|
|
N = 125.1228 - 0.0529538083 * d - longitude of the ascending node
|
|
i = 5.1454 - inclination to the ecliptic (plane of the Earth's orbit)
|
|
w = 318.0634 + 0.1643573223 * d - argument of perihelion
|
|
a = 60.2666 (Earth radii) - semi-major axis, or mean distance from Sun
|
|
e = 0.054900 - eccentricity (0=circle, 0-1=ellipse, 1=parabola)
|
|
M = 115.3654 + 13.0649929509 * d - mean anomaly (0 at perihelion; increases uniformly with time)
|
|
Related orbital elements are:
|
|
w1 = N + w - longitude of perihelion
|
|
L = M + w1 - mean longitude
|
|
q = a*(1-e) - perihelion distance
|
|
Q = a*(1+e) - aphelion distance
|
|
P = a ^ 1.5 - orbital period (years if a is in AU, astronomical units)
|
|
T = Epoch_of_M - (M(deg)/360_deg) / P - time of perihelion
|
|
v - true anomaly (angle between position and perihelion)
|
|
E - eccentric anomaly
|
|
|
|
*/
|
|
var d = daynumber(adate);
|
|
var N = d360(125.1228 - 0.0529538083 * d),
|
|
sinN = sin(N), cosN = cos(N),
|
|
i = 5.1454,
|
|
cosi = cos(i), sini = sin(i),
|
|
w = d360(318.0634 + 0.1643573223 * d),
|
|
a = 60.2666,
|
|
e = 0.054900,
|
|
M = d360(115.3654 + 13.0649929509 * d),
|
|
E = d360(M + RadToDeg * e * sin(M) * (1.0 + e*cos(M)));
|
|
|
|
// iterative correction of E
|
|
var err = 1, E0=E, iter;
|
|
for(iter = 0; iter < 20 && err > 1e-10; ++iter){
|
|
var E1 = d360(E - (E - RadToDeg * e * sin(E) - M) / (1 - e * cos(E)));
|
|
err = d360(abs(E - E1));
|
|
E = E1;
|
|
}
|
|
|
|
var xv = a * (cos(E) - e),
|
|
yv = a * sqrt(1.0 - e*e) * sin(E),
|
|
v = atan(yv, xv),
|
|
r = sqrt(xv*xv + yv*yv);
|
|
/*
|
|
Compute the planet's position in 3-dimensional space:
|
|
xh = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) )
|
|
yh = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) )
|
|
zh = r * ( sin(v+w) * sin(i) )
|
|
lonecl = atan2( yh, xh )
|
|
latecl = atan2( zh, sqrt(xh*xh+yh*yh) )
|
|
*/
|
|
// calculate sun mean anomaly & longitude; ALL in radians!
|
|
var T = calcTimeJulianCent(calcJD(adate));
|
|
var Ms = calcSunTrueAnomaly(T), // Mean Anomaly of the Sun
|
|
Ls = calcGeomMeanLongSun(T),// Mean Longitude of the Sun
|
|
Mm = M, // Mean Anomaly of the Moon
|
|
Lm = d360(N + w + M), // Mean longitude of the Moon
|
|
D = d360(Lm - Ls), // Mean elongation of the Moon
|
|
F = d360(Lm - N); // Argument of latitude for the Moon
|
|
var cosvw = cos(v + w), sinvw = sin(v + w),
|
|
xh = r * (cosN * cosvw - sinN * sinvw * cosi),
|
|
yh = r * (sinN * cosvw + cosN * sinvw * cosi),
|
|
zh = r * sinvw * sini,
|
|
elon = d360(atan(yh, xh)),
|
|
elat = atan(zh, sqrt(xh*xh + yh*yh));
|
|
|
|
var longAdd =
|
|
-1.274 * sin(Mm - 2*D)
|
|
+0.658 * sin(2*D)
|
|
-0.186 * sin(Ms)
|
|
-0.059 * sin(2*Mm - 2*D)
|
|
-0.057 * sin(Mm - 2*D + Ms)
|
|
+0.053 * sin(Mm + 2*D)
|
|
+0.046 * sin(2*D - Ms)
|
|
+0.041 * sin(Mm - Ms)
|
|
-0.035 * sin(D)
|
|
-0.031 * sin(Mm + Ms)
|
|
-0.015 * sin(2*F - 2*D)
|
|
+0.011 * sin(Mm - 4*D);
|
|
|
|
var latAdd =
|
|
-0.173 * sin(F - 2*D)
|
|
-0.055 * sin(Mm - F - 2*D)
|
|
-0.046 * sin(Mm + F - 2*D)
|
|
+0.033 * sin(F + 2*D)
|
|
+0.017 * sin(2*Mm + F);
|
|
elon += longAdd;
|
|
elat += latAdd;
|
|
r -= 0.58 * cos(Mm - 2*D) + 0.46 * cos(2*D);
|
|
// Geocentric (Earth-centered) coordinates
|
|
var xh = r * cos(elon) * cos(elat),
|
|
yh = r * sin(elon) * cos(elat),
|
|
zh = r * sin(elat);
|
|
var ecl = calcMeanObliquityOfEcliptic(T);
|
|
// Equatorial coordinates
|
|
var ye = yh * cos(ecl) - zh * sin(ecl),
|
|
ze = yh * sin(ecl) + zh * cos(ecl),
|
|
RA = d360(atan(ye, xh)),
|
|
Dec = atan(ze, sqrt(xh*xh+ye*ye));
|
|
// Horizontal coords
|
|
var GMST0=(Ls+180),
|
|
UT=d-Math.floor(d),
|
|
LST = GMST0+UT*360+longitude;
|
|
var HA = LST - RA,
|
|
x = cos(HA) * cos(Dec),
|
|
y = sin(HA) * cos(Dec),
|
|
z = sin(Dec),
|
|
xhor = x * sin(latitude) - z * cos(latitude),
|
|
zhor = x * cos(latitude) + z * sin(latitude),
|
|
az = d180(atan(y, xhor)),
|
|
zd = 90. - asin(zhor);
|
|
var MoonZen = zd - Refract(zd);
|
|
return{
|
|
azimuth: az,
|
|
zenith: MoonZen,
|
|
altitude: 90. - MoonZen,
|
|
RA: RA,
|
|
Dec: Dec
|
|
};
|
|
}
|
|
window.NOAAcalc = NOAAcalc;
|
|
}());
|