#!/usr/bin/perl
# search-grid-points.pl - search area grid waypoint computations
# by Ian Kluft
# Stratofox Aerospace Tracking Team - http://www.stratofox.org/
# for Stratofox JCS08 search, Nov 2008

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

# configuration parameters
my $waypoint_prefix = "JCS08-"; # per mission prefix: JCS = "Japanese CanSat"
my $description = "Stratofox JCS08 search grid";
my @point_start = ( 40.84253229206456, -119.080207030237 );
my $ft_per_increment = 1584;
my $h1_heading_deg = 330; # degrees
my $h2_heading_deg = 270; # degrees

# conversions on some configuration parameters
my $rad_per_increment = dist_ft2rad( $ft_per_increment );
my $h1_heading = deg2rad( $h1_heading_deg ); #radians
my $h2_heading = deg2rad( $h2_heading_deg ); #radians

# conversion: distance in feet to radians
sub dist_ft2rad
{
	my $ft = shift;
	my $km = $ft / 3280.8399; # ft/km
	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
}

# project a waypoint in the search area
# shape of a parallelogram with sides at headings 030 (NNE) and 090 (east)
# sides are on heading 1 (030 degrees) and heading 2 (090 degrees)
# increments are 0-10
# each increment is 1584 ft so that 10 of them is 3 miles
sub project_waypoint
{
	my $h1_inc = shift;
	my $h2_inc = shift;
	
	# compute intermediate point on the first side of parallelogram
	my ( $lat_r1, $lon_r1 ) = gc_waypoint (
		deg2rad( $point_start[0]), deg2rad( $point_start[1]), 
		$h1_heading, $rad_per_increment * $h1_inc );

	# compute final projected waypoint in search area
	my ( $lat_r2, $lon_r2 ) = gc_waypoint (
		 $lat_r1, $lon_r1,
		 $h2_heading, $rad_per_increment * $h2_inc );

	# convert radians to degrees
	my $lat = rad2deg( $lat_r2 );
	my $lon = rad2deg( $lon_r2 );
	return ( $lat, $lon );
}

# process command-line parameters
my $use_gpsbabel_output=0;
GetOptions( "gpsbabel" => \$use_gpsbabel_output );

# main loop - project a 10x10 grid of waypoints for the search area
my ( $h1, $h2 ); # loop indices
if ( $use_gpsbabel_output ) {
	print "Name, Latitude, Longitude, Description\n";
}
for ( $h1 = 0; $h1 <= 10; $h1++ ) {
	for ( $h2 = 0; $h2 <= 10; $h2++ ) {
		my ( $lat, $lon ) = project_waypoint ( $h1, $h2 );
		if ( $use_gpsbabel_output ) {
			printf "%s%c%c,%8.6f,%8.6f,%s %c%c\n",
				$waypoint_prefix, ord("A")+$h1, ord("A")+$h2,
				$lat, $lon,
				$description, ord("A")+$h1, ord("A")+$h2;
				
		} else {
			printf "%s%c%c\t%12.6f\t%12.6f\n",
				$waypoint_prefix, ord("A")+$h1, ord("A")+$h2, 
				$lat, $lon;
		}
	}
}