initial commit

1 files changed, 146 insertions, 0 deletions

diff --git a/gridcal.c b/gridcal.c
new file mode 100644
index 0000000..97e8cf9
--- /dev/null
+++ b/gridcal.c
@@ -0,0 +1,146 @@
+/**
+*
+* 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);
+	}
+}
+