/Proj4/geod_for.c
C | 107 lines | 106 code | 1 blank | 0 comment | 23 complexity | 92f137bd6e14fa48b50ed8ed04bc33f6 MD5 | raw file
- #ifndef lint
- static const char SCCSID[]="@(#)geod_for.c 4.6 95/09/23 GIE REL";
- #endif
- # include "projects.h"
- # include "geodesic.h"
- # define MERI_TOL 1e-9
- static double
- th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1;
- static int
- merid, signS;
- void
- geod_pre(void) {
- al12 = adjlon(al12); /* reduce to +- 0-PI */
- signS = fabs(al12) > HALFPI ? 1 : 0;
- th1 = ellipse ? atan(onef * tan(phi1)) : phi1;
- costh1 = cos(th1);
- sinth1 = sin(th1);
- if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {
- sina12 = 0.;
- cosa12 = fabs(al12) < HALFPI ? 1. : -1.;
- M = 0.;
- } else {
- cosa12 = cos(al12);
- M = costh1 * sina12;
- }
- N = costh1 * cosa12;
- if (ellipse) {
- if (merid) {
- c1 = 0.;
- c2 = f4;
- D = 1. - c2;
- D *= D;
- P = c2 / D;
- } else {
- c1 = geod_f * M;
- c2 = f4 * (1. - M * M);
- D = (1. - c2)*(1. - c2 - c1 * M);
- P = (1. + .5 * c1 * M) * c2 / D;
- }
- }
- if (merid) s1 = HALFPI - th1;
- else {
- s1 = (fabs(M) >= 1.) ? 0. : acos(M);
- s1 = sinth1 / sin(s1);
- s1 = (fabs(s1) >= 1.) ? 0. : acos(s1);
- }
- }
- void
- geod_for(void) {
- double d,sind,u,V,X,ds,cosds,sinds,ss,de;
- ss = 0.;
-
- if (ellipse) {
- d = geod_S / (D * geod_a);
- if (signS) d = -d;
- u = 2. * (s1 - d);
- V = cos(u + d);
- X = c2 * c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.);
- ds = d + X - 2. * P * V * (1. - 2. * P * cos(u)) * sind;
- ss = s1 + s1 - ds;
- } else {
- ds = geod_S / geod_a;
- if (signS) ds = - ds;
- }
- cosds = cos(ds);
- sinds = sin(ds);
- if (signS) sinds = - sinds;
- al21 = N * cosds - sinth1 * sinds;
- if (merid) {
- phi2 = atan( tan(HALFPI + s1 - ds) / onef);
- if (al21 > 0.) {
- al21 = PI;
- if (signS)
- de = PI;
- else {
- phi2 = - phi2;
- de = 0.;
- }
- } else {
- al21 = 0.;
- if (signS) {
- phi2 = - phi2;
- de = 0;
- } else
- de = PI;
- }
- } else {
- al21 = atan(M / al21);
- if (al21 > 0)
- al21 += PI;
- if (al12 < 0.)
- al21 -= PI;
- al21 = adjlon(al21);
- phi2 = atan(-(sinth1 * cosds + N * sinds) * sin(al21) /
- (ellipse ? onef * M : M));
- de = atan2(sinds * sina12 ,
- (costh1 * cosds - sinth1 * sinds * cosa12));
- if (ellipse)
- if (signS)
- de += c1 * ((1. - c2) * ds +
- c2 * sinds * cos(ss));
- else
- de -= c1 * ((1. - c2) * ds -
- c2 * sinds * cos(ss));
- }
- lam2 = adjlon( lam1 + de );
- }