/Proj4/geod_for.c
C | 107 lines | 106 code | 1 blank | 0 comment | 23 complexity | 92f137bd6e14fa48b50ed8ed04bc33f6 MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)geod_for.c 4.6 95/09/23 GIE REL"; 3#endif 4# include "projects.h" 5# include "geodesic.h" 6# define MERI_TOL 1e-9 7 static double 8th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1; 9 static int 10merid, signS; 11 void 12geod_pre(void) { 13 al12 = adjlon(al12); /* reduce to +- 0-PI */ 14 signS = fabs(al12) > HALFPI ? 1 : 0; 15 th1 = ellipse ? atan(onef * tan(phi1)) : phi1; 16 costh1 = cos(th1); 17 sinth1 = sin(th1); 18 if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) { 19 sina12 = 0.; 20 cosa12 = fabs(al12) < HALFPI ? 1. : -1.; 21 M = 0.; 22 } else { 23 cosa12 = cos(al12); 24 M = costh1 * sina12; 25 } 26 N = costh1 * cosa12; 27 if (ellipse) { 28 if (merid) { 29 c1 = 0.; 30 c2 = f4; 31 D = 1. - c2; 32 D *= D; 33 P = c2 / D; 34 } else { 35 c1 = geod_f * M; 36 c2 = f4 * (1. - M * M); 37 D = (1. - c2)*(1. - c2 - c1 * M); 38 P = (1. + .5 * c1 * M) * c2 / D; 39 } 40 } 41 if (merid) s1 = HALFPI - th1; 42 else { 43 s1 = (fabs(M) >= 1.) ? 0. : acos(M); 44 s1 = sinth1 / sin(s1); 45 s1 = (fabs(s1) >= 1.) ? 0. : acos(s1); 46 } 47} 48 void 49geod_for(void) { 50 double d,sind,u,V,X,ds,cosds,sinds,ss,de; 51 ss = 0.; 52 53 if (ellipse) { 54 d = geod_S / (D * geod_a); 55 if (signS) d = -d; 56 u = 2. * (s1 - d); 57 V = cos(u + d); 58 X = c2 * c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.); 59 ds = d + X - 2. * P * V * (1. - 2. * P * cos(u)) * sind; 60 ss = s1 + s1 - ds; 61 } else { 62 ds = geod_S / geod_a; 63 if (signS) ds = - ds; 64 } 65 cosds = cos(ds); 66 sinds = sin(ds); 67 if (signS) sinds = - sinds; 68 al21 = N * cosds - sinth1 * sinds; 69 if (merid) { 70 phi2 = atan( tan(HALFPI + s1 - ds) / onef); 71 if (al21 > 0.) { 72 al21 = PI; 73 if (signS) 74 de = PI; 75 else { 76 phi2 = - phi2; 77 de = 0.; 78 } 79 } else { 80 al21 = 0.; 81 if (signS) { 82 phi2 = - phi2; 83 de = 0; 84 } else 85 de = PI; 86 } 87 } else { 88 al21 = atan(M / al21); 89 if (al21 > 0) 90 al21 += PI; 91 if (al12 < 0.) 92 al21 -= PI; 93 al21 = adjlon(al21); 94 phi2 = atan(-(sinth1 * cosds + N * sinds) * sin(al21) / 95 (ellipse ? onef * M : M)); 96 de = atan2(sinds * sina12 , 97 (costh1 * cosds - sinth1 * sinds * cosa12)); 98 if (ellipse) 99 if (signS) 100 de += c1 * ((1. - c2) * ds + 101 c2 * sinds * cos(ss)); 102 else 103 de -= c1 * ((1. - c2) * ds - 104 c2 * sinds * cos(ss)); 105 } 106 lam2 = adjlon( lam1 + de ); 107}