[BR-Crater] KML for the ring and code that generated it (was: estimate for crater outline...)

[Ian Kluft]

Here's the KML file used to overlay the circle on the Google Map.  Obviously
you can also use this in Google Earth.
   http://ian.kluft.com/blackrock/impact-crater/br-crater-est-circle-82km.kml

I wrote the following perl program to generate that KML file as I experimented
with parameters that would fit the terrain.

------------------------------------------------------------------------
#!/usr/bin/perl
# br-est-circle-gen.pl - Black Rock circular feature estimate KML generation
# by Ian Kluft

use strict;
use POSIX qw( fmod );
use Math::Trig qw( asin acos deg2rad rad2deg :pi :great_circle );

# configuration parameters
my @midpoint = ( 40.900, -118.940 );
my $circle_points = 100;
my $diameter_km = 82;

# conversion: distance in km to radians
sub dist_km2rad
{
	my $km = shift;
	return $km/6371.0; # divide by FAI standard Earth radius in km
}

# compute lat/lon for a point given the course and distance
# from Aviation Formulary http://williams.best.vwh.net/avform.htm
sub gc_waypoint
{
	my $lat1 = shift; # latitude (radians)
	my $lon1 = shift; # longitude (radians)
	my $tc = shift; # true course (radians)
	my $d = shift; # distance (radians)

	my $lat = asin(sin($lat1)*cos($d)+cos($lat1)*sin($d)*cos($tc));
	my $dlon = atan2(sin($tc)*sin($d)*cos($lat1),cos($d)-sin($lat1)
		*sin($lat));

	my $lon=fmod($lon1-$dlon+pi,2*pi) - pi;

	return ( $lat, $lon ); # lat/lon in radians
}

# main loop - project points around the circle and feed it to gpsbabel
if ( ! open CSVPIPE, "|gpsbabel -t -i unicsv -f - -o 'kml,lines=1,points=0,line_width=25,line_color=2000ffff,labels=0,trackdata=1' -F -" ) {
	die "$0: failed to open pipe: $!\n";
}

print CSVPIPE "Latitude, Longitude\n";

my $i;
for ( $i = 0; $i < $circle_points; $i++ ) {
	my ( $lat, $lon ) = gc_waypoint (
		deg2rad($midpoint[0]), deg2rad($midpoint[1]),
		$i * 2.0 * pi / 100.0, dist_km2rad( $diameter_km / 2 ));
	printf CSVPIPE "%10.8f,%10.8f\n", rad2deg($lat), rad2deg($lon);
}
close CSVPIPE;