# ------------------------- Location functions.

# Parse a location and calculate the latitude, longitude and
# precision of the location. The function returns 
# ($e, $n, $ne, $nn, $type) where $e is the longitude, $n is
# the latitude, $ne,$nn is the size of the square in which the point
# could lie. $type returns the type of locator entered, QRA, NGR, LATLONG
# or UNPARSED.
sub parse
{
	local ($e, $n, $ne, $nn, $type, @n, @e);

	# Make things easier for searches and substitutions.
	$_ = $_[0];

	# Remove all white space, and uppercase all
	# letters.
	s/\s*//g;
	tr/a-z/A-Z/;

	# Is it a QRA locator? (e.g. IO82HN)
	if (/^[A-R]{2}[0-9]{2}[A-X]{2}$/)
	{
		($e, $n, $ne, $nn) = &qra2latlong($_);
		$type = "QRA";
	}
	# Is it an NGR? (e.g. SO097988)
	elsif  ((/^[A-HJ-Z]{2}[0-9]{2}$/) ||
			(/^[A-HJ-Z]{2}[0-9]{4}$/) ||
			(/^[A-HJ-Z]{2}[0-9]{6}$/) ||
			(/^[A-HJ-Z]{2}[0-9]{8}$/) ||
			(/^[A-HJ-Z]{2}[0-9]{10}$/))
	{
		($e, $n, $ne, $nn) = &ngr2latlong($_);
		$type = "NGR";
	}
	# Latitude Longitude in decimal form (e.g. 52.5N4.4W)
	elsif (/^([0-9]+\.?[0-9]*)([NS])([0-9]+\.?[0-9]*)([EW])$/)
	{
		# Take N/S, E/W into account.
		$n = $1 * ($2 eq "S" ? -1 : 1);
		$e = $3 * ($4 eq "W" ? -1 : 1);

		return (0, 0, 0, 0, "UNPARSED")
			if (($n > 90) || ($n < -90) || ($e > 180) || ($e < -180));

		# Extract the precision: Zero all digits, remove any
		# trailing point and convert the last digit to a 1.
		$nn = $1;
		$ne = $3;

		$nn =~ s/[0-9]/0/g;
		$nn =~ s/\.$//;
		$nn =~ s/0$/1/;

		$ne =~ s/[0-9]/0/g;
		$ne =~ s/\.$//;
		$ne =~ s/0$/1/;

		$type = "LATLONG";
	}
	# Latitude Longitude in DMS form (e.g. 52'12'14N14'4'4W)
	elsif (/^([0-9\.\']+)([NS])([0-9\.\']+)([EW])$/)
	{
		$n = $1;
		$dn= $2;
		$e = $3;
		$de= $4;

		return (0, 0, 0, 0, "UNPARSED") if ($n =~ /\'\'/);
		$n =~ s/\'$//;
		return (0, 0, 0, 0, "UNPARSED") if ($n =~ /\..*\'/);
		@n = split('\'', $n);
		return (0, 0, 0, 0, "UNPARSED") if ($#n > 2);
		$nn= $n[$#n];
		$nn=~ s/[0-9]/0/g;
		$nn=~ s/\.$//;
		$nn=~ s/0$/1/;
		$nn= $nn*((1/60)**($#n));
		$n = ($n[0] + $n[1]/60 + $n[2]/3600)*($dn eq "S" ? -1 : 1);

		return (0, 0, 0, 0, "UNPARSED") if ($e =~ /\'\'/);
		$e =~ s/\'$//;
		return (0, 0, 0, 0, "UNPARSED") if ($e =~ /\..*\'/);
		@e = split('\'', $e);
		return (0, 0, 0, 0, "UNPARSED") if ($#e > 2);
		$ne= $e[$#e];
		$ne=~ s/[0-9]/0/g;
		$ne=~ s/\.$//;
		$ne=~ s/0$/1/;
		$ne= $ne*((1/60)**($#n));
		$e = ($e[0] + $e[1]/60 + $e[2]/3600)*($de eq "W" ? -1 : 1);

		return (0, 0, 0, 0, "UNPARSED")
			if (($n > 90) || ($n < -90) || ($e > 180) || ($e < -180));

		$type = "LATLONG";
	}
	else
	{
		return (0, 0, 0, 0, "UNPARSED");
	}

	return ($e, $n, $ne, $nn, $type);
}

# ------------------------- QRA Locator functions.

# QRA locators - example: IO82HN
# I = 20s of degrees east of 180W [A-R]
# O = 10s of degrees north of 90S [A-R]
# 8 = 2s of degrees east of I     [0-9]
# 2 = 1s of degrees north of O    [0-9]
# H = 5s of minutes east of I8    [A-X]
# N = 2.5s of minutes north of O2 [A-X]

# Return the latitude and longitude, and the size of the square in 
# degrees.
sub qra2latlong
{
	local (@l) = unpack('cccccc', $_[0]);
	local ($n, $e);

	# Calculate the latitude and longitude.
	$e = (20*($l[0]-65) + 2*($l[2]-48) + 5*($l[4]-65)/60)  - 180;
	$n = (10*($l[1]-65) +   ($l[3]-48) + 5*($l[5]-65)/120) - 90;

	($e, $n, 5/60, 5/120);
}

# Return the QRA locator square in which the specified point (degrees)
# lies.
sub latlong2qra
{
	local ($e, $n) = @_;

	$e += 180;
	$l[0] = int($e/20);		$e = $e - 20*$l[0];
	$l[2] = int($e/2);		$e = $e - 2 *$l[2];
	$l[4] = int($e*60/5);

	$n += 90;
	$l[1] = int($n/10);		$n = $n - 10*$l[1];
	$l[3] = int($n);		$n = $n -    $l[3];
	$l[5] = int($n*120/5);

	sprintf("%c%c%c%c%c%c", $l[0]+65, $l[1]+65, $l[2]+48, $l[3]+48,
							$l[4]+65, $l[5]+65);
}

# ------------------------- NGR functions.

# NGRs - example: SO 097 988
# S = 500km square.
# O = 100km square.

# 0 = 10s of km east of S
# 9 = 1s of km east of S0
# 7 = 0.1s of km east of S09

# 9 = 10s of km north of O
# 8 = 1s of km north of O9
# 8 = 0.1s of km north of O98

# Return the latitude and longitude, and the size of the square in 
# degrees.
sub ngr2latlong
{
	local ($s500, $s100, @ngr) = split('', $_[0]);
	local ($eastings, $northings, $ss);

	# Decode the 500km square.
	$s500 = unpack('c', $s500) - 65;
	$s500-- if ($s500 > 8);
	$northings = 500000 * (3 - int($s500 / 5));
	$eastings  = 500000 * (($s500 % 5) - 2);

	# Decode the 100km square.
	$s100 = unpack('c', $s100) - 65;
	$s100-- if ($s100 > 8);
	$northings += 100000 * (4 - int($s100 / 5));
	$eastings  += 100000 * ($s100 % 5);

	# Figure out the precision, or return an error.
	$ss = 10**((9-$#ngr)/2);
	if ($ss ne int($ss))
	{
		die "I can only cope with 2, 4, 6, 8 and 10 figure NGRs\n";
	}

	# Extract the figures from the string to complete the
	# eastings and northings calculation.
	$eastings  += (substr($_[0], 2, (1+$#ngr)/2) * $ss);
	$northings += (substr($_[0], 2+(1+$#ngr)/2, (1+$#ngr)/2) * $ss);

	# Now convert these figures to latitude and longitude.
	($e, $n) = &en2latlong($eastings, $northings);

	# Also convert the other corner of the square into latitude
	# and longitude to get a precision.
	($ne, $nn) = &en2latlong($eastings+$ss, $northings+$ss);

	return ($e, $n, $ne-$e, $nn-$n);
}

# Convert latitude and longitude to NGR.
sub latlong2ngr
{
	local ($e, $n) = @_;

	# Find the eastings and northings in metres.
	($e, $n) = &latlong2en($e, $n);

	# Move the figure from relative to the origin at SV000000 to 
	# relative to the grid bottom left at VV000000.
	$e += 1000000;
	$n += 500000;

	# Make sure it lies within the grid.
	if (($e < 0) || ($n < 0) || ($n >= 2500000) || ($e >= 2500000))
	{
		return "-- --- ---";
	}

	# Find the 500km square.
	$s500 = int($e/500000) + 5*(4-int($n/500000));
	$s500++ if ($s500 > 7);
	$e = $e % 500000;
	$n = $n % 500000;

	# Find the 100km square.
	$s100 = int($e/100000) + 5*(4-int($n/100000));
	$s100++ if ($s100 > 7);
	$e = $e % 100000;
	$n = $n % 100000;

	# Create the 6 figure reference.
	sprintf("%c%c %03d %03d", 
				$s500+65, $s100+65, 0.5+$e/100, 0.5+$n/100);
}

# ------------------------- NGR functions (2).

# Decode metres East and North to latitude and longitude.
# Hacked from routines in _Amateur Radio Software_ by John Morris.
sub en2latlong
{
	local ($e, $n) = @_;
	local ($t1, $t2);

	$t1 = ($n/1000+5548.79) / 6371.28;
	$t2 = 2*atan2(exp(($e/1000 - 400) / 6389.70), 1);

	$e  = &deg * atan2(-cos($t2) / (cos($t1)*sin($t2)), 1) - 2;
	$n  = sin($t2)*sin($t1);
	$n  = &deg * atan2($n / sqrt(1-$n*$n), 1);

	($e, $n);
}

# Encode latitude and longitude to metres East and North.
# Also hacked, from the same book.
sub latlong2en
{
	local ($e, $n) = @_;
	local ($e1, $n1);

	$e = &rad*($e+2);
	$n = &rad*$n;

	$e1 = cos($n) * sin($e);
	$e1 = &pi/4 + atan2($e1/sqrt(1-$e1*$e1), 1)/2;
	$e1 = 6389.70 * log(sin($e1)/cos($e1)) + 400;
	$n1 = 6371.28 * atan2(sin($n)/(cos($n)*cos($e)), 1) - 5548.79;

	($e1*1000, $n1*1000);
}

# ------------------------- DMS pretty print.
sub latlong2dms
{
	local ($e, $n) = @_;
	local ($d, $m, $s, $p, $dms);

	($p, $n) = ($n < 0) ? ('S', -$n) : ('N', $n);
	$n = $n + 1/7200;
	$d = int($n);
	$n = 60 * ($n - $d);
	$m = int($n);
	$n = 60 * ($n - $m);
	$s = int($n);
	$dms = sprintf("%2d'%02d'%02d'%s", $d, $m, $s, $p);

	($p, $e) = ($e < 0) ? ('W', -$e) : ('E', $e);
	$e = $e + 1/7200;
	$d = int($e);
	$e = 60 * ($e - $d);
	$m = int($e);
	$e = 60 * ($e - $m);
	$s = int($e);
	$dms = $dms.sprintf(" %2d'%02d'%02d'%s", $d, $m, $s, $p);

	$dms;
}

# ------------------------- Distance and bearing

sub distance
{
	local ($he, $hn, $e, $n) = @_;
	local ($si, $co, $ca, $dx);

	# Convert to radians.
	$hn = &rad*$hn;
	$he = &rad*$he;
	$n  = &rad*$n;
	$e  = &rad*$e;

	# Calculate the distance.
	$co = cos($he-$e)*cos($hn)*cos($n) + sin($hn)*sin($n);
	$ca = &acos($co);
	$dx = 6367*$ca;

	# If the distance is reasonable, calculate the
	# bearing. Otherwise return it as zero.
	if ($dx > 1e-3)
	{
		$si = sin($e-$he)*cos($n)*cos($hn);
		$co = sin($n) - sin($hn)*cos($ca);
		$az = atan2(($si > 0 ? $si : -$si), ($co > 0 ? $co : -$co));
		$az = &pi - $az if ($co < 0);
		$az = -$az if ($si < 0);
		$az = $az + 2*&pi if ($az < 0);
		$az = &deg*$az;
	}
	else
	{
		$az = 0;
	}

	($dx, $az);
}

# ------------------------- Distance and bearing

sub pi { 3.141592653589793; }
sub deg { 180 / &pi; }
sub rad { &pi / 180; }

sub acos
{
	local ($ac, $co) = (0, $_[0]);

	$ac = atan2(sqrt(1-$co*$co), $co);
	$ac = &pi-$ac if ($co < 0);

	return $ac;
}

sub asin
{
	local ($si) = $_[0];

	return atan2($si, sqrt(1-$si*$si));
}

1;