The current solar terminator (the non-sunlit part of the Earth) is shown in blue. (Origin: Ben Elsen’s fix of Mike Bostock’s original)
-
-
Save johan/4645501 to your computer and use it in GitHub Desktop.
Earth, night and day sides.
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