benelsen/README.md

## README.md

      
        

          Raw
        
        
          
            
          
          
            
              README.md
            
          
        
      
    
  
    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!

  

  

## index.html
<!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>

## thumbnail.png

      
        

          Raw
        
        
          
            
          
          
            
              thumbnail.png
	<!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) * (TT) ) % 360) constants.deg2rad;
	L = ((279.69668 + ((36000.76892 * T) % 360.0) + 0.0003025 * (TT) ) % 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 - (ee)) / (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.xeci.x) + (eci.yeci.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 (2Math.PI) GMST / 86400;
	};

	</script>