geotaxis/geotaxis.c

#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <errno.h>

struct coord {
    double lat;
    double lon;
};

const char *usage_str = "Usage: %s LATITUDE LONGITUDE DISTANCE BEARING\n\
\n\
Calculate destination from starting LATITUDE and LONGITUDE and going DISTANCE in BEARING direction.\n\
\n\
Arguments:\n\
\t* BEARING is degrees deasil (i.e. clockwise) from North.\n\
\t* DIRECTION is in kilometres.\n\
\t* All arguments will be parsed as doubles, i.e. use 30.5 not 30°30'0\".\n";

void usage(FILE *out, const char *whoami) {
    fprintf(out, usage_str, whoami);
}

double parse_arg(const char *arg) {
    char *endptr;
    errno = 0;
    double val = strtod(arg, &endptr);
    if(errno != 0) {
        fprintf(stderr, "Error %d: %s\n(While attempting to parse '%s')\n", errno, strerror(errno), arg);
        exit(1);
    }
    if(endptr == arg) {
        fprintf(stderr, "Error while parsing '%s': not a number?\n", arg);
        exit(1);
    }
    return val;
}

double deg2rad(double deg) {
    return deg*M_PI/180.0;
}

double rad2deg(double rad) {
    return rad*180.0/M_PI;
}

void forward_haversine(struct coord *dest, const struct coord *origin, const double D, const double b) {
    /*
        Based on:
        * https://stackoverflow.com/a/7835325/10210798
        * https://www.movable-type.co.uk/scripts/latlong.html

        This function expects and returns radians.
    */
    
    // Save calculating the same thing more than once
    double d = D/6378.1;
    double cos_lat_1 = cos(origin->lat);
    double sin_lat_1 = sin(origin->lat);
    double cos_d = cos(d);
    double sin_d = sin(d);

    dest->lat = asin(sin_lat_1*cos_d+cos_lat_1*sin_d*cos(b));
    dest->lon = origin->lon + atan2(sin(b)*sin_d*cos_lat_1, cos_d-sin_lat_1*sin(dest->lat));
}

int main(int argc, char **argv) {
    if(argc != 5) {
        usage(stderr, argv[0]);
        exit(1);
    }
    struct coord origin, dest;
    origin.lat = deg2rad(parse_arg(argv[1]));
    origin.lon = deg2rad(parse_arg(argv[2]));
    double distance = parse_arg(argv[3]);
    double bearing = deg2rad(parse_arg(argv[4]));
    forward_haversine(&dest, &origin, distance, bearing);
    printf("%f %f\n", rad2deg(dest.lat), rad2deg(dest.lon));
    return 0;
}