#define PROJ_PARMS__ \ double n; \ double rho_c; \ double rho_0; \ double sig; \ double c1, c2; \ int type; #define PJ_LIB__ #include #define EULER 0 #define MURD1 1 #define MURD2 2 #define MURD3 3 #define PCONIC 4 #define TISSOT 5 #define VITK1 6 #define EPS10 1.e-10 #define EPS 1e-10 #define LINE2 "\n\tConic, Sph\n\tlat_1= and lat_2=" PROJ_HEAD(tissot, "Tissot") LINE2; PROJ_HEAD(murd1, "Murdoch I") LINE2; PROJ_HEAD(murd2, "Murdoch II") LINE2; PROJ_HEAD(murd3, "Murdoch III") LINE2; PROJ_HEAD(euler, "Euler") LINE2; PROJ_HEAD(pconic, "Perspective Conic") LINE2; PROJ_HEAD(vitk1, "Vitkovsky I") LINE2; /* get common factors for simple conics */ static int phi12(PJ *P, double *del) { double p1, p2; int err = 0; if (!pj_param(P->params, "tlat_1").i || !pj_param(P->params, "tlat_2").i) { err = -41; } else { p1 = pj_param(P->params, "rlat_1").f; p2 = pj_param(P->params, "rlat_2").f; *del = 0.5 * (p2 - p1); P->sig = 0.5 * (p2 + p1); err = (fabs(*del) < EPS || fabs(P->sig) < EPS) ? -42 : 0; *del = *del; } return err; } FORWARD(s_forward); /* spheroid */ double rho; switch (P->type) { case MURD2: rho = P->rho_c + tan(P->sig - lp.phi); break; case PCONIC: rho = P->c2 * (P->c1 - tan(lp.phi)); break; default: rho = P->rho_c - lp.phi; break; } xy.x = rho * sin( lp.lam *= P->n ); xy.y = P->rho_0 - rho * cos(lp.lam); return (xy); } INVERSE(s_inverse); /* ellipsoid & spheroid */ double rho; rho = hypot(xy.x, xy.y = P->rho_0 - xy.y); if (P->n < 0.) { rho = - rho; xy.x = - xy.x; xy.y = - xy.y; } lp.lam = atan2(xy.x, xy.y) / P->n; switch (P->type) { case PCONIC: lp.phi = atan(P->c1 - rho / P->c2) + P->sig; break; case MURD2: lp.phi = P->sig - atan(rho - P->rho_c); break; default: lp.phi = P->rho_c - rho; } return (lp); } FREEUP; if (P) pj_dalloc(P); } static PJ * setup(PJ *P) { double del, cs; int i; if( (i = phi12(P, &del)) ) E_ERROR(i); switch (P->type) { case TISSOT: P->n = sin(P->sig); cs = cos(del); P->rho_c = P->n / cs + cs / P->n; P->rho_0 = sqrt((P->rho_c - 2 * sin(P->phi0))/P->n); break; case MURD1: P->rho_c = sin(del)/(del * tan(P->sig)) + P->sig; P->rho_0 = P->rho_c - P->phi0; P->n = sin(P->sig); break; case MURD2: P->rho_c = (cs = sqrt(cos(del))) / tan(P->sig); P->rho_0 = P->rho_c + tan(P->sig - P->phi0); P->n = sin(P->sig) * cs; break; case MURD3: P->rho_c = del / (tan(P->sig) * tan(del)) + P->sig; P->rho_0 = P->rho_c - P->phi0; P->n = sin(P->sig) * sin(del) * tan(del) / (del * del); break; case EULER: P->n = sin(P->sig) * sin(del) / del; del *= 0.5; P->rho_c = del / (tan(del) * tan(P->sig)) + P->sig; P->rho_0 = P->rho_c - P->phi0; break; case PCONIC: P->n = sin(P->sig); P->c2 = cos(del); P->c1 = 1./tan(P->sig); if (fabs(del = P->phi0 - P->sig) - EPS10 >= HALFPI) E_ERROR(-43); P->rho_0 = P->c2 * (P->c1 - tan(del)); break; case VITK1: P->n = (cs = tan(del)) * sin(P->sig) / del; P->rho_c = del / (cs * tan(P->sig)) + P->sig; P->rho_0 = P->rho_c - P->phi0; break; } P->inv = s_inverse; P->fwd = s_forward; P->es = 0; return (P); } ENTRY0(euler) P->type = EULER; ENDENTRY(setup(P)) ENTRY0(tissot) P->type = TISSOT; ENDENTRY(setup(P)) ENTRY0(murd1) P->type = MURD1; ENDENTRY(setup(P)) ENTRY0(murd2) P->type = MURD2; ENDENTRY(setup(P)) ENTRY0(murd3) P->type = MURD3; ENDENTRY(setup(P)) ENTRY0(pconic) P->type = PCONIC; ENDENTRY(setup(P)) ENTRY0(vitk1) P->type = VITK1; ENDENTRY(setup(P))