The current solar terminator is shown in blue, assuming a spherical Earth and an axial tilt of 23.4°. Hmm, on second thought, I think I didn’t account for the orbit of the Earth around the sun, so this is not entirely accurate. Please fork this example and fix it!
-
-
Save benelsen/4597883 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| <!DOCTYPE html> | |
| <meta charset="utf-8"> | |
| <style> | |
| .night { | |
| stroke: steelblue; | |
| fill: steelblue; | |
| fill-opacity: .3; | |
| } | |
| </style> | |
| <body> | |
| <script src="http://d3js.org/d3.v3.min.js"></script> | |
| <script src="http://d3js.org/d3.geo.projection.v0.min.js"></script> | |
| <script src="http://d3js.org/topojson.v0.min.js"></script> | |
| <script> | |
| var width = 960, | |
| height = 500; | |
| var projection = d3.geo.cylindricalEqualArea() | |
| .parallel(38.5) | |
| .scale(196) | |
| .precision(.1); | |
| var circle = d3.geo.circle() | |
| .angle(89.9975572); | |
| var path = d3.geo.path() | |
| .projection(projection); | |
| var svg = d3.select("body").append("svg") | |
| .attr("width", width) | |
| .attr("height", height); | |
| d3.json("/d/4090846/world-50m.json", function(error, world) { | |
| svg.append("path") | |
| .datum(topojson.object(world, world.objects.land)) | |
| .attr("class", "land") | |
| .attr("d", path); | |
| var night = svg.append("path") | |
| .attr("class", "night") | |
| .attr("d", path); | |
| redraw(); | |
| setInterval(redraw, 250); | |
| function redraw() { | |
| var sunPos = getSolarWGSPosition( Date.now() ); | |
| var darknessAngle = 90 - Math.asin( (constants.meanRsun - constants.meanRearth) / sunPos.range ) * constants.rad2deg; | |
| night.datum(circle.origin([ -180+sunPos.lon, -sunPos.lat ]).angle(darknessAngle)).attr("d", path); | |
| } | |
| }); | |
| var getSolarWGSPosition = function getSolarWGSPosition(time) { | |
| var eci_pos = getSolarECI(time); | |
| var wgs_pos = ECItoWGS84(eci_pos, time); | |
| return {lat: wgs_pos.latitude, lon: wgs_pos.longitude, range: eci_pos.w}; | |
| }; | |
| var getSolarECI = function getSolarECI(time) { | |
| var mjd, year, T, M, L, e, C, O, Lsa, nu, R, eps; | |
| var jd_utc = unix2jd( time / 1000 ); | |
| mjd = jd_utc - 2415020.0; | |
| year = 1900 + mjd / 365.25; | |
| T = (mjd + constants.deltaUTCTT / 86400) / 36525.0; | |
| M = (( 358.47583 + ((35999.04975 * T) % 360) - (0.000150 + 0.0000033 * T) * (T*T) ) % 360) * constants.deg2rad; | |
| L = ((279.69668 + ((36000.76892 * T) % 360.0) + 0.0003025 * (T*T) ) % 360) * constants.deg2rad; | |
| e = 0.01675104 - (0.0000418 + 0.000000126 * T) * T; | |
| C = ((1.919460 - (0.004789 + 0.000014 * T) * T) * Math.sin(M) + (0.020094 - 0.000100 * T) * Math.sin(2 * M) + 0.000293 * Math.sin(3 * M)) * constants.deg2rad; | |
| O = ((259.18 - 1934.142 * T) % 360) * constants.deg2rad; | |
| Lsa = (L + C -((0.00569 - 0.00479 * Math.sin(O)) * constants.deg2rad)) % (2*Math.PI); | |
| nu = (M + C) % (2*Math.PI); | |
| R = 1.0000002 * (1.0 - (e*e)) / (1.0 + e * Math.cos(nu)); | |
| eps = ((23.452294 - (0.0130125 + (0.00000164 - 0.000000503 * T) * T) * T + 0.00256 * Math.cos(O)) * constants.deg2rad); | |
| R = constants.au * R; | |
| var x = R * Math.cos (Lsa); | |
| var y = R * Math.sin (Lsa) * Math.cos (eps); | |
| var z = R * Math.sin (Lsa) * Math.sin (eps); | |
| var w = R; | |
| return { x : x, y : y, z : z, w : w } | |
| }; | |
| var ECItoWGS84 = function ECItoWGS84(eci, time) { | |
| var jd_utc = unix2jd( time / 1000 ); | |
| var theta = Math.atan2(eci.y, eci.x); // radians | |
| var lon = (theta - thetaG_JD(jd_utc)) % (2*Math.PI); // radians | |
| var r = Math.sqrt( (eci.x*eci.x) + (eci.y*eci.y) ); | |
| var e2 = constants.f * (2 - constants.f); | |
| var lat = Math.atan2(eci.z, r); // radians | |
| var sin_phi, phi, c; | |
| do { | |
| phi = lat; | |
| sin_phi = Math.sin(phi); | |
| c = 1 / Math.sqrt(1 - e2 * (sin_phi*sin_phi)); | |
| lat = Math.atan2(eci.z + constants.eqRearth * c * e2 * sin_phi, r); | |
| } while (Math.abs(lat - phi) > 1e-10); | |
| var alt = r / Math.cos(lat) - constants.eqRearth * c; // kilometers | |
| if (lat > (Math.PI / 2)) { | |
| lat -= (2 * Math.PI); | |
| } | |
| if (lon < -Math.PI) { | |
| lon += 2*Math.PI; | |
| } | |
| return { | |
| latitude: lat * constants.rad2deg, | |
| longitude: lon * constants.rad2deg, | |
| altitude: alt, | |
| theta: theta * constants.rad2deg | |
| }; | |
| }; | |
| // ** Astronomical, Geodetic & Mathematical Constants ** // | |
| var constants = { | |
| au: 149597870700, // [m] Astronomical unit | |
| deltaUTCTT: 67.184, | |
| deg2rad: 0.017453292519943295, | |
| rad2deg: 57.29577951308232, | |
| omega_E: 1.00273790934, | |
| f: 1/298.257222101, // Earth, reciprocal of flattening IERS 2010 | |
| eqRearth: 6378.1366, // Equatorial radius | |
| meanRearth: 6371.0, // Mean radius | |
| meanRsun: 1.392684e8, // Mean radius sun | |
| omega: 7.292115e-5, // Nominal mean angular vel. of Earth rotation | |
| }; | |
| // Converts a UNIX timestamp to JD (Julian Date) | |
| var unix2jd = function unix2jd(timestamp) { | |
| return (timestamp / 86400.0) + 2440587.5; | |
| }; | |
| var thetaG_JD = function thetaG_JD(jd) { | |
| // Reference: The 1992 Astronomical Almanac, page B6. | |
| var UT, TU, GMST; | |
| UT = (jd + 0.5) % 1; | |
| jd = jd - UT; | |
| TU = (jd - 2451545.0) / 36525; | |
| GMST = 24110.54841 + TU * (8640184.812866 + TU * (0.093104 - TU * 6.2e-6)); | |
| GMST = (GMST + 86400 * constants.omega_E * UT) % 86400; | |
| return (2*Math.PI) * GMST / 86400; | |
| }; | |
| </script> |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment