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