diff options
Diffstat (limited to 'navit/tools/gpx2navit_txt/src/geod_for.c')
| -rw-r--r-- | navit/tools/gpx2navit_txt/src/geod_for.c | 138 |
1 files changed, 68 insertions, 70 deletions
diff --git a/navit/tools/gpx2navit_txt/src/geod_for.c b/navit/tools/gpx2navit_txt/src/geod_for.c index ebff3bcb4..30fa829e3 100644 --- a/navit/tools/gpx2navit_txt/src/geod_for.c +++ b/navit/tools/gpx2navit_txt/src/geod_for.c @@ -26,102 +26,100 @@ static const char SCCSID[] = # 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) -{ +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.; + sina12 = 0.; + cosa12 = fabs(al12) < HALFPI ? 1. : -1.; + M = 0.; } else { - cosa12 = cos(al12); - M = costh1 * sina12; + 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) { + 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; + s1 = HALFPI - th1; else { - s1 = (fabs(M) >= 1.) ? 0. : acos(M); - s1 = sinth1 / sin(s1); - s1 = (fabs(s1) >= 1.) ? 0. : acos(s1); + s1 = (fabs(M) >= 1.) ? 0. : acos(M); + s1 = sinth1 / sin(s1); + s1 = (fabs(s1) >= 1.) ? 0. : acos(s1); } } -void geod_for(void) -{ +void geod_for(void) { double d, sind, u, V, X, ds, cosds, sinds, ss = 0, de; 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; + 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; + ds = geod_S / geod_a; + if (signS) + ds = -ds; } cosds = cos(ds); sinds = sin(ds); if (signS) - sinds = -sinds; + 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; - } + 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)); - } + 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); } |