/*
 *   Copyright (C) 1999-2000 by Jonathan Naylor G4KLX
 *
 *   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.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <stdio.h>
#include <string.h>
#include <math.h>

#include <gtk/gtk.h>

#include "global.h"

double RAD(double angle)
{
	return (angle / 180.0) * PI;
}

double DEG(double angle)
{
	return (angle / PI) * 180.0;
}

int Convert_Locator(char *buffer, double *latitude, double *longitude)
{
	if (strlen(buffer) != 6 && strlen(buffer) != 8)
		return FALSE;

	if (buffer[0] < 'A' || buffer[0] > 'R' ||
	    buffer[1] < 'A' || buffer[1] > 'R' ||
	    buffer[2] < '0' || buffer[2] > '9' ||
	    buffer[3] < '0' || buffer[3] > '9' ||
	    buffer[4] < 'A' || buffer[4] > 'X' ||
	    buffer[5] < 'A' || buffer[5] > 'X')
		return FALSE;

	if (strlen(buffer) == 8 &&
	   (buffer[6] < '0' || buffer[6] > '9' ||
	    buffer[7] < '0' || buffer[7] > '9'))
		return FALSE;

	*longitude = -180.0 + FIELD_WIDTH   * (buffer[0] - 'A') +
			      SQUARE_WIDTH  * (buffer[2] - '0') +
			      SUB_WIDTH     * (buffer[4] - 'A');

	*latitude  = -90.0  + FIELD_HEIGHT  * (buffer[1] - 'A') +
			      SQUARE_HEIGHT * (buffer[3] - '0') +
			      SUB_HEIGHT    * (buffer[5] - 'A');

	if (strlen(buffer) == 8) {
		*longitude += SUB_SUB_WIDTH  * (buffer[6] - '0') + SUB_SUB_WIDTH / 2.0;
		*latitude  += SUB_SUB_HEIGHT * (buffer[7] - '0') + SUB_SUB_HEIGHT / 2.0;
	} else {
		*longitude += SUB_SUB_WIDTH * 5.0;
		*latitude  += SUB_SUB_HEIGHT * 5.0;
	}

	return TRUE;
}

int Convert_Angles(double Latitude, double Longitude, char *buffer)
{
	int ILatitude, ILongitude;

	if (Latitude > 90.0 || Latitude < -90.0)
		return FALSE;

	if (Longitude > 180.0 || Longitude < -180.0)
		return FALSE;

	Latitude  += 90.0;
	Longitude += 180.0;

	ILongitude = (int)(Longitude / FIELD_WIDTH);
	buffer[0]  = ILongitude + 'A';
	Longitude -= (double)ILongitude * FIELD_WIDTH;

	ILongitude = (int)(Longitude / SQUARE_WIDTH);
	buffer[2]  = ILongitude + '0';
	Longitude -= (double)ILongitude * SQUARE_WIDTH;

	ILongitude = (int)(Longitude / SUB_WIDTH);
	buffer[4]  = ILongitude + 'A';

	ILatitude = (int)(Latitude / FIELD_HEIGHT);
	buffer[1] = ILatitude + 'A';
	Latitude -= (double)ILatitude * FIELD_HEIGHT;

	ILatitude = (int)(Latitude / SQUARE_HEIGHT);
	buffer[3] = ILatitude + '0';
	Latitude -= (double)ILatitude * SQUARE_HEIGHT;

	ILatitude = (int)(Latitude / SUB_HEIGHT);
	buffer[5] = ILatitude + 'A';

	buffer[6] = '\0';

	return TRUE;
}

/*	Full distance calculation on elliptical earth (non-optimised).
	John Morris, GM4ANB, 1994.
*/
	
void Calc_Distance_Bearing(double lat1, double lon1, double lat2, double lon2, double *bearing, double *distance)
{
	double e, si, co, th, ca, a1, a2, l1, l2, ep;	/* Locals */
	double hangle;
	int rev = FALSE;

	lat1 = RAD(lat1);
	lat2 = RAD(lat2);
	lon1 = RAD(lon1);
	lon2 = RAD(lon2);

	e = 1.0 - (RB * RB) / (RA * RA);		/* Eccentricity squared */

	/* Convert latitudes to geocentric */
	lat1 = atan((RB * RB) / (RA * RA) * tan(lat1));
	lat2 = atan((RB * RB) / (RA * RA) * tan(lat2));

	/* Calculate central angle */
	ca = acos(sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon2 - lon1));

	/* Calculate angles along great ellipse from equator */
	if (lat1 > lat2) {
		double t;

		t = lat2; 
		lat2 = lat1;
		lat1 = t;

		rev = TRUE;
	}

	si = sin(lat1) * sin(ca);
	co = sin(lat2) - sin(lat1) * cos(ca);
	a1 = atan2(si, co);
	a2 = a1 + ca;

	/* Calculate ellipticity (squared) of great circle */
	th = sin(lat2) / sin(a2);
	ep = e * th * th;

	/* Calculate great circle distance from equator for each */
	l1 = a1 - ep * a1 / 4.0 + ep * sin(2.0 * a1) / 8.0;
	l2 = a2 - ep * a2 / 4.0 + ep * sin(2.0 * a2) / 8.0;
	hangle = l2 - l1;

	*distance = RA * hangle;	/* Station distance is difference */ 
	*bearing  = DEG(acos((sin(lat2) - sin(lat1) * cos(hangle)) / (cos(lat1) * sin(hangle))));

	if (lon1 > lon2 && !rev)
		*bearing = 360.0 - *bearing;
	else if (lon1 > lon2 && rev)
		*bearing = 180.0 + *bearing;
	else if (lon1 < lon2 && rev)
		*bearing = 180.0 - *bearing;
}

int Calc_Scatter_Point(double lat1, double long1, double bear1S, double lat2, double long2, double bear2S, double *slat, double *slong, double *dist1S, double *dist2S)
{
	double bear12, bear21, dist12, dist21;
	double lat1r, long1r, bear1Sr, bear12r;
	double lat2r, long2r, bear2Sr, bear21r;
	double bear1dr, bear2dr, dist12r, dist21r;
	double dist1Sr, dist2Sr, slatr, slongr;
	
	double za, zg, zb;
	double ye;

	Calc_Distance_Bearing(lat1, long1, lat2, long2, &bear12, &dist12);
	Calc_Distance_Bearing(lat2, long2, lat1, long1, &bear21, &dist21);

	lat1r = RAD(lat1); long1r = RAD(long1); bear1Sr = RAD(bear1S);
	lat2r = RAD(lat2); long2r = RAD(long2); bear2Sr = RAD(bear2S);

	bear12r = RAD(bear12); bear21r = RAD(bear21);
	dist12r = dist12 / RA;  dist21r = dist21 / RA;

	bear1dr = fabs(bear12r - bear1Sr);
	bear2dr = fabs(bear21r - bear2Sr);

	if (bear1dr == 0.0 || bear1dr == PI)
		return FALSE;

	if (bear1dr >= PI)
		bear1dr = (2.0 * PI) - bear1dr;

	if (bear2dr >= PI)
		bear2dr = (2.0 * PI) - bear2dr;

	zg = -cos(bear1dr) * cos(bear2dr) + sin(bear1dr) * sin(bear2dr) * cos(dist12r);
	ye = atan2(sqrt(1.0 - pow(zg, 2.0)), zg);
	if (ye < 0.0) ye += PI;

	za = (cos(bear2dr) + cos(bear1dr) * cos(ye)) / (sin(bear1dr) * sin(ye));
	dist1Sr = atan2(sqrt(1.0 - pow(za, 2.0)), za);
	if (dist1Sr < 0.0) dist1Sr += PI;

	zb = (cos(bear1dr) + cos(bear2dr) * cos(ye)) / (sin(bear2dr) * sin(ye));
	dist2Sr = atan2(sqrt(1.0 - pow(zb, 2.0)), zb);
	if (dist2Sr < 0.0) dist2Sr += PI;

	if (bear1Sr == 0.0) {
		slongr = long1r + dist1Sr;
		slatr  = lat1r;
	} else if (bear1Sr == PI) {
		slongr = long1r - dist1Sr;
		slatr  = lat1r;
	} else if (bear2Sr == 0.0) {
		slongr = long2r + dist2Sr;
		slatr  = lat2r;
	} else if (bear2Sr == PI) {
		slongr = long2r - dist2Sr;
		slatr  = lat2r;
	} else {
		double ix, s1, zc, xm, zf, zh;

		if (bear1Sr < PI)
			ix = bear1Sr;
		else if (bear1Sr > PI)
			ix = (2.0 * PI) - bear1Sr;
		else
			return FALSE;

		s1 = (PI / 2.0) - long1r;
		
		zc = cos(dist1Sr) * cos(s1) + sin(dist1Sr) * sin(s1) * cos(ix);
		xm = atan2(sqrt(1.0 - pow(zc, 2.0)), zc);

		zf = sin(dist1Sr) * sin(ix) / sin(xm);
		zh = atan2(zf, sqrt(1.0 - pow(zf, 2.0)));
		
		slongr = (PI / 2.0) - xm;

		if (bear1Sr < PI)
			slatr = lat1r + zh;
		else
			slatr = lat1r - zh;
	}

	*slat  = DEG(slatr);
	*slong = DEG(slongr);

	*dist1S = dist1Sr * RA;
	*dist2S = dist2Sr * RA;

	return TRUE;
}