Loading...
Searching...
No Matches
geodesic.c
Go to the documentation of this file.
#include <math.h>
#include <grass/gis.h>
#include "pi.h"
/*
* This code is preliminary. I don't know if it is even
* correct.
*/
/*
* From "Map Projections" by Peter Richardus and Ron K. Alder, 1972
* (526.8 R39m in Map & Geography Library)
* page 19, formula 2.17
*
* Formula is the equation of a geodesic from (lat1,lon1) to (lat2,lon2)
* Input is lon, output is lat (all in degrees)
*
* Note formula only works if 0 < abs(lon2-lon1) < 180
* If lon1 == lon2 then geodesic is the merdian lon1
* (and the formula will fail)
* if lon2-lon1=180 then the geodesic is either meridian lon1 or lon2
*/
/* TODO:
*
* integrate code from raster/r.surf.idw/ll.c
*/
static void adjust_lat(double *);
static void adjust_lon(double *);
static struct state {
double A, B;
} state;
static struct state *st = &state;
int G_begin_geodesic_equation(double lon1, double lat1, double lon2,
double lat2)
{
double sin21, tan1, tan2;
adjust_lon(&lon1);
adjust_lon(&lon2);
adjust_lat(&lat1);
adjust_lat(&lat2);
if (lon1 > lon2) {
double temp;
temp = lon1;
lon1 = lon2;
lon2 = temp;
temp = lat1;
lat1 = lat2;
lat2 = temp;
}
if (lon1 == lon2) {
st->A = st->B = 0.0;
return 0;
}
lon1 = Radians(lon1);
lon2 = Radians(lon2);
lat1 = Radians(lat1);
lat2 = Radians(lat2);
sin21 = sin(lon2 - lon1);
tan1 = tan(lat1);
tan2 = tan(lat2);
st->A = (tan2 * cos(lon1) - tan1 * cos(lon2)) / sin21;
st->B = (tan2 * sin(lon1) - tan1 * sin(lon2)) / sin21;
return 1;
}
/* only works if lon1 < lon < lon2 */
double G_geodesic_lat_from_lon(double lon)
{
adjust_lon(&lon);
lon = Radians(lon);
return Degrees(atan(st->A * sin(lon) - st->B * cos(lon)));
}
static void adjust_lon(double *lon)
{
while (*lon > 180.0)
*lon -= 360.0;
while (*lon < -180.0)
*lon += 360.0;
}
static void adjust_lat(double *lat)
{
if (*lat > 90.0)
*lat = 90.0;
if (*lat < -90.0)
*lat = -90.0;
}
G_begin_geodesic_equation
int G_begin_geodesic_equation(double lon1, double lat1, double lon2, double lat2)
Definition geodesic.c:38
G_geodesic_lat_from_lon
double G_geodesic_lat_from_lon(double lon)
Definition geodesic.c:78
pi.h
Degrees
#define Degrees(x)
Definition pi.h:7
Radians
#define Radians(x)
Definition pi.h:6