#!/usr/bin/perl -w
# $Id: make_small-circle,v 1.5 2004-04-23 09:36:07-06 braup Exp braup $
# Copyright 2004 Bruce Raup  (http://cires.colorado.edu/~braup/)
# This software is distributed under the GNU General Public License.  See
# http://www.gnu.org/licenses/gpl.txt for the text of this license.
# In short:
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
($progname = $0)        =~ s!^.*/!!;  # get basename of program

use Getopt::Std;

$version = '$Revision: 1.5 $';
($version) = $version =~ /^\$Revision:\s*(\d+\.\d*)/;

$spacing = 100; # spacing between generated points on the circle, km
$deg_size = 111.195;        # number of km in a degree
$earth_radius = 6371;   # km

$usage = "$progname version $version

This program depends on, and is intended to be used with, GMT (Generic Mapping
Tools (http://gmt.soest.hawaii.edu/)).

This program prints a list of geographic coordinates lying on a (small or
great) circle centered on an input longitude/latitude, with a radius specified
either as an angle or a distance.  The points can then be plotted on a map
using GMT's psxy command.

Usage:
  $progname  -h            (prints this help message and exits)
OR
  $progname  [-v] -x pole_longitude -y pole_latitude [-s point_spacing]
                 [-a circle_size_degrees | -d circle_radius_km]
where
  -x    specifies the longitude of the pole (center of circle)
  -y    specifies the latitude of the pole (center of circle)
  -a    radius of circle as an angle, in degrees (90 for great circle)
  -d    radius of circle as distance, in km (1 degree = $deg_size km is assumed)
  -s    spacing of points on the circle (km) [$spacing]
  -v    specifies verbose mode

Output is printed to stdout.  Note that if an argument to an option is a
negative number, there should be no space between it and the option letter, or
else the negative number will be interpreted as another option.
";

$opt_h = $opt_v = $opt_x = $opt_y = $opt_a = $opt_d = $opt_s = '';   # to keep perl -w quiet
if (! getopts('hvx:y:a:d:s:') ) {
  die "$usage\n";
}

die "$usage\n" if $opt_h;

warn "$progname version $version\n" if $opt_v;

if ($opt_x eq '' || $opt_y eq '') {
  die "Must specify -x and -y.  Use \"$progname -h\" for help.\n";
}
if ($opt_a eq '' && $opt_d eq '') {
  die "Must specify -a or -d.  Use \"$progname -h\" for help.\n";
}

$plon = $opt_x;
$plat = $opt_y;
$circlon = $plon;
$spacing = $opt_s if $opt_s;
$angle = $opt_a if $opt_a;

if ($opt_a) {
  if ($opt_d) {
    die "Must specify -a or -d, not both.  Use \"$progname -h\" for help.\n";
  }
  if ($opt_a < 0 || $opt_a > 90) {
    die "circle_size_degrees needs to be between 0 and 90, inclusive\n";
  }
  $circlat = $plat < 0 ? $plat + $opt_a : $plat - $opt_a;
}

if ($opt_d) {
  if ($opt_a) {
    die "Must specify -a or -d, not both.  Use \"$progname -h\" for help.\n";
  }
  if ($opt_d < 0 || $opt_d > 11000) {
    die "circle_size_km needs to be between 0 and 11000, inclusive\n";
  }
  $angle = $opt_d / $deg_size;
  $circlat = $plat < 0 ? $plat + $angle : $plat - $angle;
}

# calculate perimeter of circle as 2 pi R sin(a)  (a in radians)
$perim = 6.28318530718 * $earth_radius * sin($angle*0.0174532925);

if ($opt_v) {
  warn "Pole longitude, latitude                = $plon, $plat\n";
  warn "Point on circle at zero longitude       = $circlon, $circlat\n";
}

system( "project -T$plon/$plat -C$circlon/$circlat -L0/$perim -G$spacing -Q" ) &&
   die "project had a problem:  $!";