diff options
| -rw-r--r-- | ChangeLog | 4 | ||||
| -rw-r--r-- | src/PJ_aitoff.c | 88 |
2 files changed, 89 insertions, 3 deletions
@@ -1,4 +1,8 @@ 2015-02-16 Howard Butler <howard@hobu.co> + * src/PJ_aitoff.c: New Winkel Tripel Inverse #250 + from Drazen Tutic + +2015-02-16 Howard Butler <howard@hobu.co> * src/pj_init.c: Fix #237 warning about initialization ordering due to setlocale diff --git a/src/PJ_aitoff.c b/src/PJ_aitoff.c index c8d2b56b..45f77d2d 100644 --- a/src/PJ_aitoff.c +++ b/src/PJ_aitoff.c @@ -36,8 +36,8 @@ PJ_CVSID("$Id$"); -PROJ_HEAD(aitoff, "Aitoff") "\n\tMisc Sph, no inv."; -PROJ_HEAD(wintri, "Winkel Tripel") "\n\tMisc Sph, no inv.\n\tlat_1"; +PROJ_HEAD(aitoff, "Aitoff") "\n\tMisc Sph"; +PROJ_HEAD(wintri, "Winkel Tripel") "\n\tMisc Sph\n\tlat_1"; FORWARD(s_forward); /* spheroid */ double c, d; @@ -53,10 +53,92 @@ FORWARD(s_forward); /* spheroid */ } return (xy); } + +/*********************************************************************************** +* +* Inverse functions added by Drazen Tutic and Lovro Gradiser based on paper: +* +* I.Özbug Biklirici and Cengizhan Ipbüker. A General Algorithm for the Inverse +* Transformation of Map Projections Using Jacobian Matrices. In Proceedings of the +* Third International Symposium Mathematical & Computational Applications, +* pages 175{182, Turkey, September 2002. +* +* Expected accuracy is defined by EPSILON = 1e-12. Should be appropriate for +* most applications of Aitoff and Winkel Tripel projections. +* +* Longitudes of 180W and 180E can be mixed in solution obtained. +* +* Inverse for Aitoff projection in poles is undefined, longitude value of 0 is assumed. +* +* Contact : dtutic@geof.hr +* Date: 2014-10-31 +* +************************************************************************************/ + +INVERSE(s_inverse); /* sphere */ + int iter, MAXITER = 20, round = 0, MAXROUND = 10; + double EPSILON = 1e-12, D, C, f1, f2, f1p, f1l, f2p, f2l, dp, dl, sl, sp, cp, cl, x, y; + + if ((fabs(xy.x) < EPSILON) && (fabs(xy.y) < EPSILON )) { lp.phi = 0.; lp.lam = 0.; return (lp); } + + /* intial values for Newton-Raphson method */ + lp.phi = xy.y; lp.lam = xy.x; + do { + iter = 0; + do { + sl = sin(lp.lam * 0.5); cl = cos(lp.lam * 0.5); + sp = sin(lp.phi); cp = cos(lp.phi); + D = cp * cl; + C = 1. - D * D; + D = acos(D) / pow(C, 1.5); + f1 = 2. * D * C * cp * sl; + f2 = D * C * sp; + f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl); + f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp; + f2p = sp * sp * cl / C + D * sl * sl * cp; + f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl); + if (P->mode) { /* Winkel Tripel */ + f1 = 0.5 * (f1 + lp.lam * P->cosp1); + f2 = 0.5 * (f2 + lp.phi); + f1p *= 0.5; + f1l = 0.5 * (f1l + P->cosp1); + f2p = 0.5 * (f2p + 1.); + f2l *= 0.5; + } + f1 -= xy.x; f2 -= xy.y; + dl = (f2 * f1p - f1 * f2p) / (dp = f1p * f2l - f2p * f1l); + dp = (f1 * f2l - f2 * f1l) / dp; + while (dl > M_PI) dl -= M_PI; /* set to interval [-M_PI, M_PI] */ + while (dl < -M_PI) dl += M_PI; /* set to interval [-M_PI, M_PI] */ + lp.phi -= dp; lp.lam -= dl; + } while ((dp > EPSILON || dl > EPSILON) && (iter++ < MAXITER)); + + if (lp.phi > M_PI_2) lp.phi -= 2.*(lp.phi-M_PI_2); /* correct if symmetrical solution for Aitoff */ + if (lp.phi < -M_PI_2) lp.phi -= 2.*(lp.phi+M_PI_2); /* correct if symmetrical solution for Aitoff */ + if ((fabs(fabs(lp.phi) - M_PI_2) < EPSILON) && (!P->mode)) lp.lam = 0.; /* if pole in Aitoff, return longitude of 0 */ + + /* calculate x,y coordinates with solution obtained */ + if((D = acos(cos(lp.phi) * cos(C = 0.5 * lp.lam)))) {/* Aitoff */ + x = 2. * D * cos(lp.phi) * sin(C) * (y = 1. / sin(D)); + y *= D * sin(lp.phi); + } else + x = y = 0.; + if (P->mode) { /* Winkel Tripel */ + x = (x + lp.lam * P->cosp1) * 0.5; + y = (y + lp.phi) * 0.5; + } + /* if too far from given values of x,y, repeat with better approximation of phi,lam */ + } while (((fabs(xy.x-x) > EPSILON) || (fabs(xy.y-y) > EPSILON)) && (round++ < MAXROUND)); + + if (iter == MAXITER && round == MAXROUND) fprintf(stderr, "Max number of iterations reached! Last iteration increments dlat.=%f and dlon.=%f\n", dp, dl); + + return (lp); +} + FREEUP; if (P) pj_dalloc(P); } static PJ * setup(PJ *P) { - P->inv = 0; + P->inv = s_inverse; P->fwd = s_forward; P->es = 0.; return P; |