Skip to content

Commit

Permalink
Merge 8066dcd into 6c16367
Browse files Browse the repository at this point in the history
  • Loading branch information
cffk committed Aug 16, 2015
2 parents 6c16367 + 8066dcd commit 31f18b5
Show file tree
Hide file tree
Showing 5 changed files with 232 additions and 199 deletions.
23 changes: 14 additions & 9 deletions man/man1/geod.1
Expand Up @@ -30,16 +30,22 @@ file[s]
]
file[s]
.SH DESCRIPTION
.I Geod
.I geod
(direct) and
.I invgeod
(inverse)
perform geodesic (\*(lqGreat Circle\*(rq) computations for determining
perform geodesic ("Great Circle") computations for determining
latitude, longitude and back azimuth of a terminus point
given a initial point latitude, longitude, azimuth and distance (direct) or
the forward and back azimuths and distance between an initial and
terminus point latitudes and longitudes (inverse). The results are
accurate to round off for |\fIf\fR| < 1/50, where \fIf\fR is flattening.
.B invgeod
may not be available on all platforms; in this case call
.B geod
with the
.B \-I
option.
.PP
The following command-line options can appear in any order:
.TP
Expand Down Expand Up @@ -120,12 +126,11 @@ input for the inverse mode and respective forward and back
azimuth from the initial and terminus points are output along
with the distance between the points.
.PP
Input geographic coordinates
(latitude and longitude) and azimuthal data must be in DMS format and input
distance data must be in units consistent with the ellipsoid
major axis or sphere radius units.
Output geographic coordinates will be in DMS
(if the
Input geographic coordinates (latitude and longitude) and azimuthal data
must be in decimal degrees or DMS format and input distance data must be
in units consistent with the ellipsoid major axis or sphere radius
units. The latitude must lie in the range [-90d,90d]. Output
geographic coordinates will be in DMS (if the
.B \-f
switch is not employed) to 0.001"
with trailing, zero-valued minute-second fields deleted.
Expand Down Expand Up @@ -220,4 +225,4 @@ The \fIonline geodesic bibliography\fR,
.br
http://geographiclib.sf.net/geodesic-papers/biblio.html
.SH HOME PAGE
http://proj.osgeo.org
https://github.com/OSGeo/proj.4/wiki
6 changes: 3 additions & 3 deletions man/man3/geodesic.3
Expand Up @@ -38,7 +38,7 @@ measure angles (latitudes, longitudes, and azimuths) in degrees, unlike
the rest of the \fBproj\fR library, which uses radians. The
documentation for this library is included in geodesic.h. A formatted
version of the documentation is available at
http://geographiclib.sf.net/1.43/C
http://geographiclib.sf.net/1.44/C
.SH EXAMPLE
The following program reads in lines with the coordinates for two points
in decimal degrees (\fIlat1\fR, \fIlon1\fR, \fIlat2\fR, \fIlon2\fR) and
Expand Down Expand Up @@ -72,7 +72,7 @@ libproj.a \- library of projections and support procedures
.SH SEE ALSO
Full online documentation for \fBgeodesic(3)\fR,
.br
http://geographiclib.sf.net/1.43/C
http://geographiclib.sf.net/1.44/C
.PP
.B geod(1)
.PP
Expand All @@ -94,4 +94,4 @@ The \fIonline geodesic bibliography\fR,
.br
http://geographiclib.sf.net/geodesic-papers/biblio.html
.SH HOME PAGE
http://proj.osgeo.org
https://github.com/OSGeo/proj.4/wiki
29 changes: 17 additions & 12 deletions src/geod_interface.c
@@ -1,34 +1,39 @@
#include "projects.h"
#include "geod_interface.h"

/* DEG_IN is a crock to work around the problem that dmstor.c uses the wrong
* value for pi/180 (namely .0174532925199433 which is an inaccurately
* truncated version of DEG_TO_RAD).
*/
#define DEG_IN .0174532925199433
#define DEG_OUT DEG_TO_RAD;

void geod_ini(void) {
geod_init(&GlobalGeodesic, geod_a, geod_f);
}

void geod_pre(void) {
double
degree = PI/180,
lat1 = phi1 / degree, lon1 = lam1 /degree, azi1 = al12 / degree;
geod_lineinit(&GlobalGeodesicLine, &GlobalGeodesic,
lat1, lon1, azi1, 0U);
lat1 = phi1 / DEG_IN, lon1 = lam1 / DEG_IN, azi1 = al12 / DEG_IN;
geod_lineinit(&GlobalGeodesicLine, &GlobalGeodesic, lat1, lon1, azi1, 0U);
}

void geod_for(void) {
double degree = PI/180, s12 = geod_S, lat2, lon2, azi2;
double
s12 = geod_S, lat2, lon2, azi2;
geod_position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
phi2 = lat2 * degree;
lam2 = lon2 * degree;
al21 = azi2 * degree;
phi2 = lat2 * DEG_OUT;
lam2 = lon2 * DEG_OUT;
al21 = azi2 * DEG_OUT;
}

void geod_inv(void) {
double
degree = PI / 180,
lat1 = phi1 / degree, lon1 = lam1 / degree,
lat2 = phi2 / degree, lon2 = lam2 / degree,
lat1 = phi1 / DEG_IN, lon1 = lam1 / DEG_IN,
lat2 = phi2 / DEG_IN, lon2 = lam2 / DEG_IN,
azi1, azi2, s12;
geod_inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
al12 = azi1 * degree; al21 = azi2 * degree; geod_S = s12;
al12 = azi1 * DEG_OUT; al21 = azi2 * DEG_OUT; geod_S = s12;
}

0 comments on commit 31f18b5

Please sign in to comment.