#! /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
}