#! /usr/bin/perl

$pi = 3.1415;
$piO2 = $pi / 2.0;

$boslat =  42.365;
$boslon = -71.100;

$boslatrad = ($boslat / 180.0) * $pi;

while(<>) {
  if(/^#/o) { print; next; }
  local($x,$y,$z,$r);
  local($lon,$lat) = split;
  $lon = $lon +0;
  $lat = $lat +0;

  # Rotate earth on axis so Boston is logitude 0.

  $lon = $lon - $boslon;
  $lon = -360 + $lon if $lon > 180;

  # Convert from degrees to radians.

  $lon = ($lon / 180.0) * $pi;
  $lat = ($lat / 180.0) * $pi;

  #
  # Determine x,y,z coordinates.
  # x is along equator, left to right, -1 to 1, W to E.
  # y is along meridian, down to up, -1 to 1, S to N.
  # z is inward, shallow to deep, 1 to -1, this side to that side.
  # So lat,lon 0,0 is (0,0,1).
  #

  $y = sin($lat);
  $r = cos($lat);
  $x = $r * sin($lon);
  $z = $r * cos($lon);

  #
  # Now rotate a latitude down to the equator.
  # Axis (x) is horizontal, parellel to view plane.
  # Boston's lat/lon 42/"0" is rotated to "0"/"0".
  #

  $rotr = sqrt(($y*$y)+($z*$z));
  $rota = atan2($y,$z);
  $rotanew = $rota - $boslatrad;
  $y = $rotr * sin($rotanew);
  $z = $rotr * cos($rotanew);

  #
  # Pointing down to points on the earth... calculate
  #  ro: counter-clockwise (from x) direction (E,N,W,S) of pointing,
  #  dn: declanation(?), 0 for level, 90 for straight down.
  #

  $ro = atan2($y,$x);
  $ror = sqrt(($x*$x) + ($y*$y));
  $dn = atan2((1 - $z),$ror);

  #
  # Generate a flat cartesian plot of the "pointing down".
  #  dn is corrected so down is in the middle of a 0 centered plot.
  #

  $polx = ($piO2 - $dn) * sin($ro);
  $poly = ($piO2 - $dn) * cos($ro);

#  print $lon," ",$lat,"\n"; # for testing
#  print $x," ",$y,"\n" if $z > 0; # looking at the earth
#  print $x," ",$z,"\n" if $y > 0; # and from the top
  print $poly," ",$polx,"\n"; # the main diagram
}