/**
*
* Get two grid locators and compute the great-circle distance
* between the centers of the two corresponding squares
*
*
* Usage: gridcal <grid1> [<grid2> [<pwr>]]
*
*/

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

#define MIN(a,b) ((a)>(b)?(b):(a))

double R = 6371;

typedef struct {
	double lat;
	double lon;
} latlon_t;


void usage(char *argv0){

	printf("Usage: %s <grid1> <grid2> [<pwr>]", argv0);
	exit(1);
}

int check_input(char *s){
	int S = strlen(s);
#ifdef __DEBUG__
	fprintf(stderr, "locator: %s (length: %d)\n", s, S);
#endif
	if (S>6 || S<4 || S%2 == 1 
	    || tolower(s[0]) < 'a' || tolower(s[0]) > 'r' 
	    || tolower(s[1]) < 'a' || tolower(s[1]) > 'r'
	    || s[2] < '0' || s[2] > '9' 
	    || s[3] < '0' || s[3] > '9'
	    || (S > 4 && (tolower(s[4]) < 'a' || tolower(s[4]) > 'x'
	    || tolower(s[5]) < 'a' || tolower(s[5]) > 'x'))) {
		printf("Invalid locator: %s\n", s);
		exit(1); 
	}
}


/**
*
* This function takes a valid locator and returns the coordinates
* of the centre of the corresponding square in the latlon_t structure
*
*/

void grid_to_latlon(char *s, latlon_t *c){

	c->lat = c->lon = 0;
	/* longitude first */	
 	c->lon += (tolower(s[0])-'a')*20.0 - 180;	
	c->lon += (s[2]-'0')*2.0;
	/* we always consider the centre of the square, not the corner...*/
	if (strlen(s)>4){
		c->lon += (tolower(s[4]) - 'a') * 1.0/12 + 1.0/24;
	}
	else {
		c->lon += 1.0;
	}

	/* then latitude */
 	c->lat += (tolower(s[1])-'a')*10.0 - 90;	
	c->lat += (s[3]-'0')*1.0;
	/* we always consider the centre of the square, not the corner...*/
	if (strlen(s)>4){
		c->lat += (tolower(s[5]) - 'a') * 1.0/24 + 1.0/48;
	}
	else {
		c->lat += 0.5;
	}
	c->lon = c->lon/180 * M_PI;
	c->lat = c->lat/90 * M_PI_2;
}

double heversine_dist(latlon_t *c1, latlon_t *c2){

	double d=0.0, d1=0, d2=0;

	d1 = sin(0.5*(c2->lat - c1->lat));
	d1 = d1 * d1;
	
	d2 = sin(0.5*(c2->lon - c1->lon));
	d2 = d2 * d2;
	d2 = d2 * cos(c1->lat) * cos(c2->lat);

	d = 2 * R * asin(sqrt(d1 + d2));
	return d;
}


int main(int argc, char *argv[]){
	
	char loc1[7], loc2[7];
	int n1, n2;
	latlon_t c1, c2;
	double d;
	double pwr;

	if (argc < 2){
		usage(argv[0]);
	}
	if (argc == 2){
		/* We only convert a grid locator to lat-long */
		n1=MIN(strlen(argv[1]),6);
		strncpy(loc1, argv[1], n1);
		loc1[n1] = '\0';
		check_input(loc1);
		grid_to_latlon(loc1, &c1);
		printf("%s %g %g\n", loc1, c1.lat, c1.lon);
		exit(0);
	}
	
	n1=MIN(strlen(argv[1]),6);
	n2=MIN(strlen(argv[2]),6);
	strncpy(loc1, argv[1], n1);
	strncpy(loc2, argv[2], n2);
	loc1[n1] = '\0';
	loc2[n2] = '\0';
	check_input(loc1);
	check_input(loc2);

	grid_to_latlon(loc1, &c1);
	grid_to_latlon(loc2, &c2);
	
	d = heversine_dist(&c1, &c2);
	if (argc > 3){
		pwr = atof(argv[3]);
		printf("%s %s %g %g %g %g %g\n", \
		       loc1, loc2, d, d/1.60934, pwr, d/pwr, d/(1.60934*pwr));
	}
	else {
		printf("%s %s %g %g\n", loc1, loc2, d, d/1.60934);
	}
}