/*

   xlog - GTK+ logging program for amateur radio operators
   Copyright (C) 2001 - 2007 Joop Stakenborg <pg4i@amsat.org>

   This file is part of xlog.

   Xlog 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 3 of the License, or
   (at your option) any later version.

   Xlog 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.

   You should have received a copy of the GNU General Public License
   along with xlog.  If not, see <http://www.gnu.org/licenses/>.

*/

/*
 * gc.c - originally by N5OWK in the Public Domain, adopted for xlog October 2001 
 */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <gtk/gtk.h>

#include "support.h"
#include "utils.h"
#include "gc.h"

#define RADIAN (180.0 / M_PI)

struct
{
	gdouble miles;		/* arc length for 1 degree, various units of measure */
	gchar *text;
}

Units[] =
{
	{60.0, "m"},
	{111.2, "km"},
};

/* Error routine */
static void
err (gint type)
{
	switch (type)
	{
	case 1:
		update_statusbar (_("Latitude Out of Range (90N to 90S)"));
		break;
	case 2:
		update_statusbar (_("Longitude Out of Range (180W to 180E)"));
		break;
	case 3:
		update_statusbar (_("Minutes Out of Range (0 to 59)"));
	}
}

/* Convert Degrees and Minutes to Decimal */
static gdouble
dm2dec (gdouble n)
{
	gdouble t;

	t = (gint) n;
	n -= t;
	n /= .60;
	if (n >= 1.0)
	{
		err (3);
		return (-1);
	}
	return (n + t);
}

/* Parse the input line dd(.mm)[NnSs]ddd(.mm)[EeWw] */
static gint
parse (gchar * s, gdouble * lat, gdouble * lon)
{
	gchar *i, *t = NULL, *e;
	gboolean first = TRUE, north = FALSE, east = FALSE;

	if (g_strrstr (s, _("N"))) north = TRUE;
	if (g_strrstr (s, _("E"))) east = TRUE;


	e = s + strlen (s);
	for (i = s; i < e; ++i)
	{
		if (g_ascii_isalpha (*i))
		{
			if (first)
			{
				*i = '\0';
				t = i + 1;
				first = FALSE;
				if (north)
					*lat = atof (s);
				else
					*lat = -atof (s);
			}
			else
			{
				*i = '\0';
				if (east)
					*lon = -atof (t);
				else
					*lon = atof (t);
			}
		}
	}

	*lat = dm2dec (*lat);
	*lon = dm2dec (*lon);

	if (*lat > 90.0 || *lat < -90.0)
	{
		err (1);
		return (-1);
	}

	if (*lon > 180.0 || *lon < -180.0)
	{
		err (2);
		return (-1);
	}

	/* Prevent ACOS() Domain Error */
	if (*lat == 90.0) *lat = 89.9;
	if (*lat == -90.0) *lat = -89.9;

	return (0);
}

gchar *
greatcircle (gint units, gchar * qth, gdouble DEST_Lat, gdouble DEST_Long)
{
	gdouble tmp, arc, cosaz, az, azsp, azlp, distsp, distlp;
	gdouble QTH_Lat = 0, QTH_Long = 0, Delta_Long;
	gint parseresult = 0;
	gchar *result = NULL;

	parseresult = parse (g_strdup (qth), &QTH_Lat, &QTH_Long);
	if (parseresult == -1) return ("");

	QTH_Lat /= RADIAN;		/* Convert variables to Radians */
	QTH_Long /= RADIAN;
	DEST_Lat /= RADIAN;
	DEST_Long /= RADIAN;

	Delta_Long = DEST_Long - QTH_Long;

	tmp = (sin (QTH_Lat) * sin (DEST_Lat)) +
		(cos (QTH_Lat) * cos (DEST_Lat) * cos (Delta_Long));

	if (tmp > .999999)
	{
		/* Station points coincide, use an Omni! */
		return ("");
	}
	else if (tmp < -.999999)
	{

/*
 * points are antipodal, he's straight down.
 * So take 180 Degrees of arc times 60 nm,
 * and you get 10800 nm, or whatever units...
 */
		return ("");
	}
	else
	{
		arc = acos (tmp);

/*
 * One degree of arc is 60 Nautical miles
 * at the surface of the earth, 111.2 km, or 69.1 sm
 * This method is easier than the one in the handbook
 */

		/* Short Path */
		distsp = (Units[units].miles) * (arc * RADIAN);
		/* Long Path */
		distlp = ((Units[units].miles) * 360.0) - distsp;
	}

	cosaz = (sin (DEST_Lat) -
		(sin (QTH_Lat) * cos (arc))) / (sin (arc) * cos (QTH_Lat));

	if (cosaz > .999999)
		az = 0.0;
	else if (cosaz < -.999999)
		az = 180.0;
	else
		az = acos (cosaz) * RADIAN;

	/* Handbook had the test ">= 0.0" which looks backwards??			 */

	if (sin (Delta_Long) < 0.0)
	{
		azsp = az;
		azlp = 180.0 + az;
	}
	else
	{
		azsp = 360.0 - az;
		azlp = 180.0 - az;
	}

	/* Computations complete, show answer */
	result = g_strdup_printf
(_("\nShort Path: %03.0f deg, %.0f %s\nLong Path: %03.0f deg, %.0f %s\n"), 
		azsp, distsp, Units[units].text, azlp, distlp, Units[units].text);
	return (result);
}