/*============================================================================ * * WCSLIB - an implementation of the FITS WCS proposal. * Copyright (C) 1995-2002, Mark Calabretta * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * * Correspondence concerning WCSLIB may be directed to: * Internet email: mcalabre@atnf.csiro.au * Postal address: Dr. Mark Calabretta, * Australia Telescope National Facility, * P.O. Box 76, * Epping, NSW, 2121, * AUSTRALIA * *============================================================================= * * C implementation of the spherical map projections recognized by the FITS * "World Coordinate System" (WCS) convention. * * Summary of routines * ------------------- * Each projection is implemented via separate functions for the forward, * *fwd(), and reverse, *rev(), transformation. * * Initialization routines, *set(), compute intermediate values from the * projection parameters but need not be called explicitly - see the * explanation of prj.flag below. * * prjset prjfwd prjrev Driver routines (see below). * * azpset azpfwd azprev AZP: zenithal/azimuthal perspective * szpset szpfwd szprev SZP: slant zenithal perspective * tanset tanfwd tanrev TAN: gnomonic * stgset stgfwd stgrev STG: stereographic * sinset sinfwd sinrev SIN: orthographic/synthesis * arcset arcfwd arcrev ARC: zenithal/azimuthal equidistant * zpnset zpnfwd zpnrev ZPN: zenithal/azimuthal polynomial * zeaset zeafwd zearev ZEA: zenithal/azimuthal equal area * airset airfwd airrev AIR: Airy * cypset cypfwd cyprev CYP: cylindrical perspective * ceaset ceafwd cearev CEA: cylindrical equal area * carset carfwd carrev CAR: Cartesian * merset merfwd merrev MER: Mercator * sflset sflfwd sflrev SFL: Sanson-Flamsteed * parset parfwd parrev PAR: parabolic * molset molfwd molrev MOL: Mollweide * aitset aitfwd aitrev AIT: Hammer-Aitoff * copset copfwd coprev COP: conic perspective * coeset coefwd coerev COE: conic equal area * codset codfwd codrev COD: conic equidistant * cooset coofwd coorev COO: conic orthomorphic * bonset bonfwd bonrev BON: Bonne * pcoset pcofwd pcorev PCO: polyconic * tscset tscfwd tscrev TSC: tangential spherical cube * cscset cscfwd cscrev CSC: COBE quadrilateralized spherical cube * qscset qscfwd qscrev QSC: quadrilateralized spherical cube * * * Driver routines; prjset(), prjfwd() & prjrev() * ---------------------------------------------- * A set of driver routines are available for use as a generic interface to * the specific projection routines. The interfaces to prjfwd() and prjrev() * are the same as those of the forward and reverse transformation routines * for the specific projections (see below). * * The interface to prjset() differs slightly from that of the initialization * routines for the specific projections and unlike them it must be invoked * explicitly to use prjfwd() and prjrev(). * * Given: * pcode[4] const char * WCS projection code. * * Given and/or returned: * prj prjprm* Projection parameters (see below). * * Function return value: * int Error status * 0: Success. * * * Initialization routine; *set() * ------------------------------ * Initializes members of a prjprm data structure which hold intermediate * values. Note that this routine need not be called directly; it will be * invoked by prjfwd() and prjrev() if the "flag" structure member is * anything other than a predefined magic value. * * Given and/or returned: * prj prjprm* Projection parameters (see below). * * Function return value: * int Error status * 0: Success. * 1: Invalid projection parameters. * * Forward transformation; *fwd() * ----------------------------- * Compute (x,y) coordinates in the plane of projection from native spherical * coordinates (phi,theta). * * Given: * phi, const double * theta Longitude and latitude of the projected point in * native spherical coordinates, in degrees. * * Given and returned: * prj prjprm* Projection parameters (see below). * * Returned: * x,y double* Projected coordinates. * * Function return value: * int Error status * 0: Success. * 1: Invalid projection parameters. * 2: Invalid value of (phi,theta). * * Reverse transformation; *rev() * ----------------------------- * Compute native spherical coordinates (phi,theta) from (x,y) coordinates in * the plane of projection. * * Given: * x,y const double * Projected coordinates. * * Given and returned: * prj prjprm* Projection parameters (see below). * * Returned: * phi, double* Longitude and latitude of the projected point in * theta native spherical coordinates, in degrees. * * Function return value: * int Error status * 0: Success. * 1: Invalid projection parameters. * 2: Invalid value of (x,y). * 1: Invalid projection parameters. * * Projection parameters * --------------------- * The prjprm struct consists of the following: * * int flag * This flag must be set to zero whenever any of p[10] or r0 are set * or changed. This signals the initialization routine to recompute * intermediaries. flag may also be set to -1 to disable strict bounds * checking for the AZP, SZP, TAN, SIN, ZPN, and COP projections. * * double r0 * r0; The radius of the generating sphere for the projection, a linear * scaling parameter. If this is zero, it will be reset to the default * value of 180/pi (the value for FITS WCS). * * double p[10] * The first 10 elements contain projection parameters which correspond * to the PROJPn keywords in FITS, so p[0] is PROJP0, and p[9] is * PROJP9. Many projections use p[1] (PROJP1) and some also use p[2] * (PROJP2). ZPN is the only projection which uses any of the others. * * The remaining members of the prjprm struct are maintained by the * initialization routines and should not be modified. This is done for the * sake of efficiency and to allow an arbitrary number of contexts to be * maintained simultaneously. * * char code[4] * Three-letter projection code. * * double phi0, theta0 * Native longitude and latitude of the reference point, in degrees. * * double w[10] * int n * Intermediate values derived from the projection parameters. * * int (*prjfwd)() * int (*prjrev)() * Pointers to the forward and reverse projection routines. * * Usage of the p[] array as it applies to each projection is described in * the prologue to each trio of projection routines. * * Argument checking * ----------------- * Forward routines: * * The values of phi and theta (the native longitude and latitude) * normally lie in the range [-180,180] for phi, and [-90,90] for theta. * However, all forward projections will accept any value of phi and will * not normalize it. * * The forward projection routines do not explicitly check that theta lies * within the range [-90,90]. They do check for any value of theta which * produces an invalid argument to the projection equations (e.g. leading * to division by zero). The forward routines for AZP, SZP, TAN, SIN, * ZPN, and COP also return error 2 if (phi,theta) corresponds to the * overlapped (far) side of the projection but also return the * corresponding value of (x,y). This strict bounds checking may be * relaxed by setting prj->flag to -1 (rather than 0) when these * projections are initialized. * * Reverse routines: * * Error checking on the projected coordinates (x,y) is limited to that * required to ascertain whether a solution exists. Where a solution does * exist no check is made that the value of phi and theta obtained lie * within the ranges [-180,180] for phi, and [-90,90] for theta. * * Accuracy * -------- * Closure to a precision of at least 1E-10 degree of longitude and latitude * has been verified for typical projection parameters on the 1 degree grid * of native longitude and latitude (to within 5 degrees of any latitude * where the projection may diverge). * * Author: Mark Calabretta, Australia Telescope National Facility * $Id: proj.c,v 1.1.1.1 2016/03/30 20:00:02 joye Exp $ *===========================================================================*/ #include #include #include #include "wcslib.h" int npcode = 26; char pcodes[26][4] = {"AZP", "SZP", "TAN", "STG", "SIN", "ARC", "ZPN", "ZEA", "AIR", "CYP", "CEA", "CAR", "MER", "COP", "COE", "COD", "COO", "SFL", "PAR", "MOL", "AIT", "BON", "PCO", "TSC", "CSC", "QSC"}; const int AZP = 101; const int SZP = 102; const int TAN = 103; const int STG = 104; const int SIN = 105; const int ARC = 106; const int ZPN = 107; const int ZEA = 108; const int AIR = 109; const int CYP = 201; const int CEA = 202; const int CAR = 203; const int MER = 204; const int SFL = 301; const int PAR = 302; const int MOL = 303; const int AIT = 401; const int COP = 501; const int COE = 502; const int COD = 503; const int COO = 504; const int BON = 601; const int PCO = 602; const int TSC = 701; const int CSC = 702; const int QSC = 703; /* Map error number to error message for each function. */ const char *prjset_errmsg[] = { 0, "Invalid projection parameters"}; const char *prjfwd_errmsg[] = { 0, "Invalid projection parameters", "Invalid value of (phi,theta)"}; const char *prjrev_errmsg[] = { 0, "Invalid projection parameters", "Invalid value of (x,y)"}; #define copysgn(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X)) #define copysgni(X, Y) ((Y) < 0 ? -abs(X) : abs(X)) /*==========================================================================*/ int prjset(pcode, prj) const char pcode[4]; struct prjprm *prj; { /* Set pointers to the forward and reverse projection routines. */ if (strcmp(pcode, "AZP") == 0) { azpset(prj); } else if (strcmp(pcode, "SZP") == 0) { szpset(prj); } else if (strcmp(pcode, "TAN") == 0) { tanset(prj); } else if (strcmp(pcode, "STG") == 0) { stgset(prj); } else if (strcmp(pcode, "SIN") == 0) { sinset(prj); } else if (strcmp(pcode, "ARC") == 0) { arcset(prj); } else if (strcmp(pcode, "ZPN") == 0) { zpnset(prj); } else if (strcmp(pcode, "ZEA") == 0) { zeaset(prj); } else if (strcmp(pcode, "AIR") == 0) { airset(prj); } else if (strcmp(pcode, "CYP") == 0) { cypset(prj); } else if (strcmp(pcode, "CEA") == 0) { ceaset(prj); } else if (strcmp(pcode, "CAR") == 0) { carset(prj); } else if (strcmp(pcode, "MER") == 0) { merset(prj); } else if (strcmp(pcode, "SFL") == 0) { sflset(prj); } else if (strcmp(pcode, "PAR") == 0) { parset(prj); } else if (strcmp(pcode, "MOL") == 0) { molset(prj); } else if (strcmp(pcode, "AIT") == 0) { aitset(prj); } else if (strcmp(pcode, "COP") == 0) { copset(prj); } else if (strcmp(pcode, "COE") == 0) { coeset(prj); } else if (strcmp(pcode, "COD") == 0) { codset(prj); } else if (strcmp(pcode, "COO") == 0) { cooset(prj); } else if (strcmp(pcode, "BON") == 0) { bonset(prj); } else if (strcmp(pcode, "PCO") == 0) { pcoset(prj); } else if (strcmp(pcode, "TSC") == 0) { tscset(prj); } else if (strcmp(pcode, "CSC") == 0) { cscset(prj); } else if (strcmp(pcode, "QSC") == 0) { qscset(prj); } else { /* Unrecognized projection code. */ return 1; } return 0; } /*--------------------------------------------------------------------------*/ int prjfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { return prj->prjfwd(phi, theta, prj, x, y); } /*--------------------------------------------------------------------------*/ int prjrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { return prj->prjrev(x, y, prj, phi, theta); } /*============================================================================ * AZP: zenithal/azimuthal perspective projection. * * Given: * prj->p[1] Distance parameter, mu in units of r0. * prj->p[2] Tilt angle, gamma in degrees. * * Given and/or returned: * prj->flag AZP, or -AZP if prj->flag is given < 0. * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "AZP" * prj->phi0 0.0 * prj->theta0 90.0 * prj->w[0] r0*(mu+1) * prj->w[1] tan(gamma) * prj->w[2] sec(gamma) * prj->w[3] cos(gamma) * prj->w[4] sin(gamma) * prj->w[5] asin(-1/mu) for |mu| >= 1, -90 otherwise * prj->w[6] mu*cos(gamma) * prj->w[7] 1 if |mu*cos(gamma)| < 1, 0 otherwise * prj->prjfwd Pointer to azpfwd(). * prj->prjrev Pointer to azprev(). *===========================================================================*/ int azpset(prj) struct prjprm *prj; { strcpy(prj->code, "AZP"); prj->flag = copysgni (AZP, prj->flag); prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) prj->r0 = R2D; prj->w[0] = prj->r0*(prj->p[1] + 1.0); if (prj->w[0] == 0.0) { return 1; } prj->w[3] = cosdeg (prj->p[2]); if (prj->w[3] == 0.0) { return 1; } prj->w[2] = 1.0/prj->w[3]; prj->w[4] = sindeg (prj->p[2]); prj->w[1] = prj->w[4] / prj->w[3]; if (fabs(prj->p[1]) > 1.0) { prj->w[5] = asindeg (-1.0/prj->p[1]); } else { prj->w[5] = -90.0; } prj->w[6] = prj->p[1] * prj->w[3]; prj->w[7] = (fabs(prj->w[6]) < 1.0) ? 1.0 : 0.0; prj->prjfwd = azpfwd; prj->prjrev = azprev; return 0; } /*--------------------------------------------------------------------------*/ int azpfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, b, cphi, cthe, r, s, t; if (abs(prj->flag) != AZP) { if (azpset(prj)) return 1; } cphi = cosdeg (phi); cthe = cosdeg (theta); s = prj->w[1]*cphi; t = (prj->p[1] + sindeg (theta)) + cthe*s; if (t == 0.0) { return 2; } r = prj->w[0]*cthe/t; *x = r*sindeg (phi); *y = -r*cphi*prj->w[2]; /* Bounds checking. */ if (prj->flag > 0) { /* Overlap. */ if (theta < prj->w[5]) { return 2; } /* Divergence. */ if (prj->w[7] > 0.0) { t = prj->p[1] / sqrt(1.0 + s*s); if (fabs(t) <= 1.0) { s = atandeg (-s); t = asindeg (t); a = s - t; b = s + t + 180.0; if (a > 90.0) a -= 360.0; if (b > 90.0) b -= 360.0; if (theta < ((a > b) ? a : b)) { return 2; } } } } return 0; } /*--------------------------------------------------------------------------*/ int azprev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double a, b, r, s, t, ycosg; const double tol = 1.0e-13; if (abs(prj->flag) != AZP) { if (azpset(prj)) return 1; } ycosg = y*prj->w[3]; r = sqrt(x*x + ycosg*ycosg); if (r == 0.0) { *phi = 0.0; *theta = 90.0; } else { *phi = atan2deg (x, -ycosg); s = r / (prj->w[0] + y*prj->w[4]); t = s*prj->p[1]/sqrt(s*s + 1.0); s = atan2deg (1.0, s); if (fabs(t) > 1.0) { t = copysgn (90.0,t); if (fabs(t) > 1.0+tol) { return 2; } } else { t = asindeg (t); } a = s - t; b = s + t + 180.0; if (a > 90.0) a -= 360.0; if (b > 90.0) b -= 360.0; *theta = (a > b) ? a : b; } return 0; } /*============================================================================ * SZP: slant zenithal perspective projection. * * Given: * prj->p[1] Distance of the point of projection from the centre of the * generating sphere, mu in units of r0. * prj->p[2] Native longitude, phi_c, and ... * prj->p[3] Native latitude, theta_c, on the planewards side of the * intersection of the line through the point of projection * and the centre of the generating sphere, phi_c in degrees. * * Given and/or returned: * prj->flag SZP, or -SZP if prj->flag is given < 0. * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "SZP" * prj->phi0 0.0 * prj->theta0 90.0 * prj->w[0] 1/r0 * prj->w[1] xp = -mu*cos(theta_c)*sin(phi_c) * prj->w[2] yp = mu*cos(theta_c)*cos(phi_c) * prj->w[3] zp = mu*sin(theta_c) + 1 * prj->w[4] r0*xp * prj->w[5] r0*yp * prj->w[6] r0*zp * prj->w[7] (zp - 1)^2 * prj->w[8] asin(1-zp) if |1 - zp| < 1, -90 otherwise * prj->prjfwd Pointer to szpfwd(). * prj->prjrev Pointer to szprev(). *===========================================================================*/ int szpset(prj) struct prjprm *prj; { strcpy(prj->code, "SZP"); prj->flag = copysgni (SZP, prj->flag); prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) prj->r0 = R2D; prj->w[0] = 1.0/prj->r0; prj->w[3] = prj->p[1] * sindeg (prj->p[3]) + 1.0; if (prj->w[3] == 0.0) { return 1; } prj->w[1] = -prj->p[1] * cosdeg (prj->p[3]) * sindeg (prj->p[2]); prj->w[2] = prj->p[1] * cosdeg (prj->p[3]) * cosdeg (prj->p[2]); prj->w[4] = prj->r0 * prj->w[1]; prj->w[5] = prj->r0 * prj->w[2]; prj->w[6] = prj->r0 * prj->w[3]; prj->w[7] = (prj->w[3] - 1.0) * prj->w[3] - 1.0; if (fabs(prj->w[3] - 1.0) < 1.0) { prj->w[8] = asindeg (1.0 - prj->w[3]); } else { prj->w[8] = -90.0; } prj->prjfwd = szpfwd; prj->prjrev = szprev; return 0; } /*--------------------------------------------------------------------------*/ int szpfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, b, cphi, cthe, s, sphi, t; if (abs(prj->flag) != SZP) { if (szpset(prj)) return 1; } cphi = cosdeg (phi); sphi = sindeg (phi); cthe = cosdeg (theta); s = 1.0 - sindeg (theta); t = prj->w[3] - s; if (t == 0.0) { return 2; } *x = (prj->w[6]*cthe*sphi - prj->w[4]*s)/t; *y = -(prj->w[6]*cthe*cphi + prj->w[5]*s)/t; /* Bounds checking. */ if (prj->flag > 0) { /* Divergence. */ if (theta < prj->w[8]) { return 2; } /* Overlap. */ if (fabs(prj->p[1]) > 1.0) { s = prj->w[1]*sphi - prj->w[2]*cphi; t = 1.0/sqrt(prj->w[7] + s*s); if (fabs(t) <= 1.0) { s = atan2deg (s, prj->w[3] - 1.0); t = asindeg (t); a = s - t; b = s + t + 180.0; if (a > 90.0) a -= 360.0; if (b > 90.0) b -= 360.0; if (theta < ((a > b) ? a : b)) { return 2; } } } } return 0; } /*--------------------------------------------------------------------------*/ int szprev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double a, b, c, d, r2, sth1, sth2, sthe, sxy, t, x1, xp, y1, yp, z; const double tol = 1.0e-13; if (abs(prj->flag) != SZP) { if (szpset(prj)) return 1; } xp = x*prj->w[0]; yp = y*prj->w[0]; r2 = xp*xp + yp*yp; x1 = (xp - prj->w[1])/prj->w[3]; y1 = (yp - prj->w[2])/prj->w[3]; sxy = xp*x1 + yp*y1; if (r2 < 1.0e-10) { /* Use small angle formula. */ z = r2/2.0; *theta = 90.0 - R2D*sqrt(r2/(1.0 + sxy)); } else { t = x1*x1 + y1*y1; a = t + 1.0; b = sxy - t; c = r2 - sxy - sxy + t - 1.0; d = b*b - a*c; /* Check for a solution. */ if (d < 0.0) { return 2; } d = sqrt(d); /* Choose solution closest to pole. */ sth1 = (-b + d)/a; sth2 = (-b - d)/a; sthe = (sth1 > sth2) ? sth1 : sth2; if (sthe > 1.0) { if (sthe-1.0 < tol) { sthe = 1.0; } else { sthe = (sth1 < sth2) ? sth1 : sth2; } } if (sthe < -1.0) { if (sthe+1.0 > -tol) { sthe = -1.0; } } if (sthe > 1.0 || sthe < -1.0) { return 2; } *theta = asindeg (sthe); z = 1.0 - sthe; } *phi = atan2deg (xp - x1*z, -(yp - y1*z)); return 0; } /*============================================================================ * TAN: gnomonic projection. * * Given and/or returned: * prj->flag TAN, or -TAN if prj->flag is given < 0. * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "TAN" * prj->phi0 0.0 * prj->theta0 90.0 * prj->prjfwd Pointer to tanfwd(). * prj->prjrev Pointer to tanrev(). *===========================================================================*/ int tanset(prj) struct prjprm *prj; { int k; strcpy(prj->code, "TAN"); prj->flag = copysgni (TAN, prj->flag); prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) prj->r0 = R2D; prj->prjfwd = tanfwd; prj->prjrev = tanrev; for (k = (MAXPV-1); k >= 0 && prj->ppv[k] == 0.0 && prj->ppv[k+MAXPV] == 0.0; k--); if (k < 0) k = 0; prj->npv = k; return 0; } /*--------------------------------------------------------------------------*/ int tanfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double r, s; double xp[2]; if (abs(prj->flag) != TAN) { if(tanset(prj)) return 1; } s = sindeg (theta); if (s <= 0.0) { return 2; } r = prj->r0*cosdeg (theta)/s; xp[0] = r*sindeg (phi); xp[1] = -r*cosdeg (phi); *x = prj->inv_x? poly_func(prj->inv_x, xp) : xp[0]; *y = prj->inv_y? poly_func(prj->inv_y, xp) : xp[1]; if (prj->flag > 0 && s < 0.0) { return 2; } return 0; } /*--------------------------------------------------------------------------*/ int tanrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double r; double xp; double yp; if (abs(prj->flag) != TAN) { if (tanset(prj)) return 1; } if (prj->npv) { raw_to_pv(prj, x,y, &xp, &yp); } else { xp = x; yp = y; } r = sqrt(xp*xp + yp*yp); if (r == 0.0) { *phi = 0.0; } else { *phi = atan2deg (xp, -yp); } *theta = atan2deg (prj->r0, r); return 0; } /*============================================================================ * STG: stereographic projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "STG" * prj->flag STG * prj->phi0 0.0 * prj->theta0 90.0 * prj->w[0] 2*r0 * prj->w[1] 1/(2*r0) * prj->prjfwd Pointer to stgfwd(). * prj->prjrev Pointer to stgrev(). *===========================================================================*/ int stgset(prj) struct prjprm *prj; { strcpy(prj->code, "STG"); prj->flag = STG; prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 360.0/PI; prj->w[1] = PI/360.0; } else { prj->w[0] = 2.0*prj->r0; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = stgfwd; prj->prjrev = stgrev; return 0; } /*--------------------------------------------------------------------------*/ int stgfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double r, s; if (prj->flag != STG) { if (stgset(prj)) return 1; } s = 1.0 + sindeg (theta); if (s == 0.0) { return 2; } r = prj->w[0]*cosdeg (theta)/s; *x = r*sindeg (phi); *y = -r*cosdeg (phi); return 0; } /*--------------------------------------------------------------------------*/ int stgrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double r; if (prj->flag != STG) { if (stgset(prj)) return 1; } r = sqrt(x*x + y*y); if (r == 0.0) { *phi = 0.0; } else { *phi = atan2deg (x, -y); } *theta = 90.0 - 2.0*atandeg (r*prj->w[1]); return 0; } /*============================================================================ * SIN: orthographic/synthesis projection. * * Given: * prj->p[1:2] Obliqueness parameters, xi and eta. * * Given and/or returned: * prj->flag SIN, or -SIN if prj->flag is given < 0. * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "SIN" * prj->phi0 0.0 * prj->theta0 90.0 * prj->w[0] 1/r0 * prj->w[1] xi**2 + eta**2 * prj->w[2] xi**2 + eta**2 + 1 * prj->w[3] xi**2 + eta**2 - 1 * prj->prjfwd Pointer to sinfwd(). * prj->prjrev Pointer to sinrev(). *===========================================================================*/ int sinset(prj) struct prjprm *prj; { strcpy(prj->code, "SIN"); prj->flag = copysgni (SIN, prj->flag); prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) prj->r0 = R2D; prj->w[0] = 1.0/prj->r0; prj->w[1] = prj->p[1]*prj->p[1] + prj->p[2]*prj->p[2]; prj->w[2] = prj->w[1] + 1.0; prj->w[3] = prj->w[1] - 1.0; prj->prjfwd = sinfwd; prj->prjrev = sinrev; return 0; } /*--------------------------------------------------------------------------*/ int sinfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double cphi, cthe, sphi, t, z; if (abs(prj->flag) != SIN) { if (sinset(prj)) return 1; } t = (90.0 - fabs(theta))*D2R; if (t < 1.0e-5) { if (theta > 0.0) { z = t*t/2.0; } else { z = 2.0 - t*t/2.0; } cthe = t; } else { z = 1.0 - sindeg (theta); cthe = cosdeg (theta); } cphi = cosdeg (phi); sphi = sindeg (phi); *x = prj->r0*(cthe*sphi + prj->p[1]*z); *y = -prj->r0*(cthe*cphi - prj->p[2]*z); /* Validate this solution. */ if (prj->flag > 0) { if (prj->w[1] == 0.0) { /* Orthographic projection. */ if (theta < 0.0) { return 2; } } else { /* "Synthesis" projection. */ t = -atandeg (prj->p[1]*sphi - prj->p[2]*cphi); if (theta < t) { return 2; } } } return 0; } /*--------------------------------------------------------------------------*/ int sinrev (x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { const double tol = 1.0e-13; double a, b, c, d, r2, sth1, sth2, sthe, sxy, x0, x1, xp, y0, y1, yp, z; if (abs(prj->flag) != SIN) { if (sinset(prj)) return 1; } /* Compute intermediaries. */ x0 = x*prj->w[0]; y0 = y*prj->w[0]; r2 = x0*x0 + y0*y0; if (prj->w[1] == 0.0) { /* Orthographic projection. */ if (r2 != 0.0) { *phi = atan2deg (x0, -y0); } else { *phi = 0.0; } if (r2 < 0.5) { *theta = acosdeg (sqrt(r2)); } else if (r2 <= 1.0) { *theta = asindeg (sqrt(1.0 - r2)); } else { return 2; } } else { /* "Synthesis" projection. */ x1 = prj->p[1]; y1 = prj->p[2]; sxy = x0*x1 + y0*y1; if (r2 < 1.0e-10) { /* Use small angle formula. */ z = r2/2.0; *theta = 90.0 - R2D*sqrt(r2/(1.0 + sxy)); } else { a = prj->w[2]; b = sxy - prj->w[1]; c = r2 - sxy - sxy + prj->w[3]; d = b*b - a*c; /* Check for a solution. */ if (d < 0.0) { return 2; } d = sqrt(d); /* Choose solution closest to pole. */ sth1 = (-b + d)/a; sth2 = (-b - d)/a; sthe = (sth1 > sth2) ? sth1 : sth2; if (sthe > 1.0) { if (sthe-1.0 < tol) { sthe = 1.0; } else { sthe = (sth1 < sth2) ? sth1 : sth2; } } if (sthe < -1.0) { if (sthe+1.0 > -tol) { sthe = -1.0; } } if (sthe > 1.0 || sthe < -1.0) { return 2; } *theta = asindeg (sthe); z = 1.0 - sthe; } xp = -y0 + prj->p[2]*z; yp = x0 - prj->p[1]*z; if (xp == 0.0 && yp == 0.0) { *phi = 0.0; } else { *phi = atan2deg (yp,xp); } } return 0; } /*============================================================================ * ARC: zenithal/azimuthal equidistant projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "ARC" * prj->flag ARC * prj->phi0 0.0 * prj->theta0 90.0 * prj->w[0] r0*(pi/180) * prj->w[1] (180/pi)/r0 * prj->prjfwd Pointer to arcfwd(). * prj->prjrev Pointer to arcrev(). *===========================================================================*/ int arcset(prj) struct prjprm *prj; { strcpy(prj->code, "ARC"); prj->flag = ARC; prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 1.0; prj->w[1] = 1.0; } else { prj->w[0] = prj->r0*D2R; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = arcfwd; prj->prjrev = arcrev; return 0; } /*--------------------------------------------------------------------------*/ int arcfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double r; if (prj->flag != ARC) { if (arcset(prj)) return 1; } r = prj->w[0]*(90.0 - theta); *x = r*sindeg (phi); *y = -r*cosdeg (phi); return 0; } /*--------------------------------------------------------------------------*/ int arcrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double r; if (prj->flag != ARC) { if (arcset(prj)) return 1; } r = sqrt(x*x + y*y); if (r == 0.0) { *phi = 0.0; } else { *phi = atan2deg (x, -y); } *theta = 90.0 - r*prj->w[1]; return 0; } /*============================================================================ * ZPN: zenithal/azimuthal polynomial projection. * * Given: * prj->p[0:9] Polynomial coefficients. * * Given and/or returned: * prj->flag ZPN, or -ZPN if prj->flag is given < 0. * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "ZPN" * prj->phi0 0.0 * prj->theta0 90.0 * prj->n Degree of the polynomial, N. * prj->w[0] Co-latitude of the first point of inflection (N > 2). * prj->w[1] Radius of the first point of inflection (N > 2). * prj->prjfwd Pointer to zpnfwd(). * prj->prjrev Pointer to zpnrev(). *===========================================================================*/ int zpnset(prj) struct prjprm *prj; { int i, j, k; double d, d1, d2, r, zd, zd1, zd2; const double tol = 1.0e-13; strcpy(prj->code, "ZPN"); prj->flag = copysgni (ZPN, prj->flag); prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) prj->r0 = R2D; /* Find the highest non-zero coefficient. */ for (k = 9; k >= 0 && prj->p[k] == 0.0; k--){ i = 0; } /* if (k < 0) return 1; */ /* if (k < 0 switch to ARC projection */ if (k < 0) { return (arcset (prj)); } prj->n = k; /* No negative derivative -> no point of inflection. */ zd = PI; /* Processing subroutines */ prj->prjfwd = zpnfwd; prj->prjrev = zpnrev; if (k >= 3) { /* Find the point of inflection closest to the pole. */ zd1 = 0.0; d1 = prj->p[1]; if (d1 <= 0.0) { return 1; } /* Find the point where the derivative first goes negative. */ for (i = 0; i < 180; i++) { zd2 = i*D2R; d2 = 0.0; for (j = k; j > 0; j--) { d2 = d2*zd2 + j*prj->p[j]; } if (d2 <= 0.0) break; zd1 = zd2; d1 = d2; } if (i == 180) { /* No negative derivative -> no point of inflection. */ zd = PI; } else { /* Find where the derivative is zero. */ for (i = 1; i <= 10; i++) { zd = zd1 - d1*(zd2-zd1)/(d2-d1); d = 0.0; for (j = k; j > 0; j--) { d = d*zd + j*prj->p[j]; } if (fabs(d) < tol) break; if (d < 0.0) { zd2 = zd; d2 = d; } else { zd1 = zd; d1 = d; } } } r = 0.0; for (j = k; j >= 0; j--) { r = r*zd + prj->p[j]; } prj->w[0] = zd; prj->w[1] = r; } return 0; } /*--------------------------------------------------------------------------*/ int zpnfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { int j; double r, s; if (abs(prj->flag) != ZPN) { if (zpnset(prj)) return 1; } s = (90.0 - theta)*D2R; r = 0.0; for (j = 9; j >= 0; j--) { r = r*s + prj->p[j]; } r = prj->r0*r; *x = r*sindeg (phi); *y = -r*cosdeg (phi); if (prj->flag > 0 && s > prj->w[0]) { return 2; } return 0; } /*--------------------------------------------------------------------------*/ int zpnrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { int i, j, k; double a, b, c, d, lambda, r, r1, r2, rt, zd, zd1, zd2; const double tol = 1.0e-13; if (abs(prj->flag) != ZPN) { if (zpnset(prj)) return 1; } k = prj->n; r = sqrt(x*x + y*y)/prj->r0; if (k < 1) { /* Constant - no solution. */ return 1; } else if (k == 1) { /* Linear. */ zd = (r - prj->p[0])/prj->p[1]; } else if (k == 2) { /* Quadratic. */ a = prj->p[2]; b = prj->p[1]; c = prj->p[0] - r; d = b*b - 4.0*a*c; if (d < 0.0) { return 2; } d = sqrt(d); /* Choose solution closest to pole. */ zd1 = (-b + d)/(2.0*a); zd2 = (-b - d)/(2.0*a); zd = (zd1zd2) ? zd1 : zd2; if (zd < 0.0) { if (zd < -tol) { return 2; } zd = 0.0; } else if (zd > PI) { if (zd > PI+tol) { return 2; } zd = PI; } } else { /* Higher order - solve iteratively. */ zd1 = 0.0; r1 = prj->p[0]; zd2 = prj->w[0]; r2 = prj->w[1]; if (r < r1) { if (r < r1-tol) { return 2; } zd = zd1; } else if (r > r2) { if (r > r2+tol) { return 2; } zd = zd2; } else { /* Disect the interval. */ for (j = 0; j < 100; j++) { lambda = (r2 - r)/(r2 - r1); if (lambda < 0.1) { lambda = 0.1; } else if (lambda > 0.9) { lambda = 0.9; } zd = zd2 - lambda*(zd2 - zd1); rt = 0.0; for (i = k; i >= 0; i--) { rt = (rt * zd) + prj->p[i]; } if (rt < r) { if (r-rt < tol) break; r1 = rt; zd1 = zd; } else { if (rt-r < tol) break; r2 = rt; zd2 = zd; } if (fabs(zd2-zd1) < tol) break; } } } if (r == 0.0) { *phi = 0.0; } else { *phi = atan2deg (x, -y); } *theta = 90.0 - zd*R2D; return 0; } /*============================================================================ * ZEA: zenithal/azimuthal equal area projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "ZEA" * prj->flag ZEA * prj->phi0 0.0 * prj->theta0 90.0 * prj->w[0] 2*r0 * prj->w[1] 1/(2*r0) * prj->prjfwd Pointer to zeafwd(). * prj->prjrev Pointer to zearev(). *===========================================================================*/ int zeaset(prj) struct prjprm *prj; { strcpy(prj->code, "ZEA"); prj->flag = ZEA; prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 360.0/PI; prj->w[1] = PI/360.0; } else { prj->w[0] = 2.0*prj->r0; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = zeafwd; prj->prjrev = zearev; return 0; } /*--------------------------------------------------------------------------*/ int zeafwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double r; if (prj->flag != ZEA) { if (zeaset(prj)) return 1; } r = prj->w[0]*sindeg ((90.0 - theta)/2.0); *x = r*sindeg (phi); *y = -r*cosdeg (phi); return 0; } /*--------------------------------------------------------------------------*/ int zearev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double r, s; const double tol = 1.0e-12; if (prj->flag != ZEA) { if (zeaset(prj)) return 1; } r = sqrt(x*x + y*y); if (r == 0.0) { *phi = 0.0; } else { *phi = atan2deg (x, -y); } s = r*prj->w[1]; if (fabs(s) > 1.0) { if (fabs(r - prj->w[0]) < tol) { *theta = -90.0; } else { return 2; } } else { *theta = 90.0 - 2.0*asindeg (s); } return 0; } /*============================================================================ * AIR: Airy's projection. * * Given: * prj->p[1] Latitude theta_b within which the error is minimized, in * degrees. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "AIR" * prj->flag AIR * prj->phi0 0.0 * prj->theta0 90.0 * prj->w[0] 2*r0 * prj->w[1] ln(cos(xi_b))/tan(xi_b)**2, where xi_b = (90-theta_b)/2 * prj->w[2] 1/2 - prj->w[1] * prj->w[3] 2*r0*prj->w[2] * prj->w[4] tol, cutoff for using small angle approximation, in * radians. * prj->w[5] prj->w[2]*tol * prj->w[6] (180/pi)/prj->w[2] * prj->prjfwd Pointer to airfwd(). * prj->prjrev Pointer to airrev(). *===========================================================================*/ int airset(prj) struct prjprm *prj; { const double tol = 1.0e-4; double cxi; strcpy(prj->code, "AIR"); prj->flag = AIR; prj->phi0 = 0.0; prj->theta0 = 90.0; if (prj->r0 == 0.0) prj->r0 = R2D; prj->w[0] = 2.0*prj->r0; if (prj->p[1] == 90.0) { prj->w[1] = -0.5; prj->w[2] = 1.0; } else if (prj->p[1] > -90.0) { cxi = cosdeg ((90.0 - prj->p[1])/2.0); prj->w[1] = log(cxi)*(cxi*cxi)/(1.0-cxi*cxi); prj->w[2] = 0.5 - prj->w[1]; } else { return 1; } prj->w[3] = prj->w[0] * prj->w[2]; prj->w[4] = tol; prj->w[5] = prj->w[2]*tol; prj->w[6] = R2D/prj->w[2]; prj->prjfwd = airfwd; prj->prjrev = airrev; return 0; } /*--------------------------------------------------------------------------*/ int airfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double cxi, r, txi, xi; if (prj->flag != AIR) { if (airset(prj)) return 1; } if (theta == 90.0) { r = 0.0; } else if (theta > -90.0) { xi = D2R*(90.0 - theta)/2.0; if (xi < prj->w[4]) { r = xi*prj->w[3]; } else { cxi = cosdeg ((90.0 - theta)/2.0); txi = sqrt(1.0-cxi*cxi)/cxi; r = -prj->w[0]*(log(cxi)/txi + prj->w[1]*txi); } } else { return 2; } *x = r*sindeg (phi); *y = -r*cosdeg (phi); return 0; } /*--------------------------------------------------------------------------*/ int airrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { int j; double cxi, lambda, r, r1, r2, rt, txi, x1, x2, xi; const double tol = 1.0e-12; if (prj->flag != AIR) { if (airset(prj)) return 1; } r = sqrt(x*x + y*y)/prj->w[0]; if (r == 0.0) { xi = 0.0; } else if (r < prj->w[5]) { xi = r*prj->w[6]; } else { /* Find a solution interval. */ x1 = 1.0; r1 = 0.0; for (j = 0; j < 30; j++) { x2 = x1/2.0; txi = sqrt(1.0-x2*x2)/x2; r2 = -(log(x2)/txi + prj->w[1]*txi); if (r2 >= r) break; x1 = x2; r1 = r2; } if (j == 30) return 2; for (j = 0; j < 100; j++) { /* Weighted division of the interval. */ lambda = (r2-r)/(r2-r1); if (lambda < 0.1) { lambda = 0.1; } else if (lambda > 0.9) { lambda = 0.9; } cxi = x2 - lambda*(x2-x1); txi = sqrt(1.0-cxi*cxi)/cxi; rt = -(log(cxi)/txi + prj->w[1]*txi); if (rt < r) { if (r-rt < tol) break; r1 = rt; x1 = cxi; } else { if (rt-r < tol) break; r2 = rt; x2 = cxi; } } if (j == 100) return 2; xi = acosdeg (cxi); } if (r == 0.0) { *phi = 0.0; } else { *phi = atan2deg (x, -y); } *theta = 90.0 - 2.0*xi; return 0; } /*============================================================================ * CYP: cylindrical perspective projection. * * Given: * prj->p[1] Distance of point of projection from the centre of the * generating sphere, mu, in units of r0. * prj->p[2] Radius of the cylinder of projection, lambda, in units of * r0. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "CYP" * prj->flag CYP * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*lambda*(pi/180) * prj->w[1] (180/pi)/(r0*lambda) * prj->w[2] r0*(mu + lambda) * prj->w[3] 1/(r0*(mu + lambda)) * prj->prjfwd Pointer to cypfwd(). * prj->prjrev Pointer to cyprev(). *===========================================================================*/ int cypset(prj) struct prjprm *prj; { strcpy(prj->code, "CYP"); prj->flag = CYP; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = prj->p[2]; if (prj->w[0] == 0.0) { return 1; } prj->w[1] = 1.0/prj->w[0]; prj->w[2] = R2D*(prj->p[1] + prj->p[2]); if (prj->w[2] == 0.0) { return 1; } prj->w[3] = 1.0/prj->w[2]; } else { prj->w[0] = prj->r0*prj->p[2]*D2R; if (prj->w[0] == 0.0) { return 1; } prj->w[1] = 1.0/prj->w[0]; prj->w[2] = prj->r0*(prj->p[1] + prj->p[2]); if (prj->w[2] == 0.0) { return 1; } prj->w[3] = 1.0/prj->w[2]; } prj->prjfwd = cypfwd; prj->prjrev = cyprev; return 0; } /*--------------------------------------------------------------------------*/ int cypfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double s; if (prj->flag != CYP) { if (cypset(prj)) return 1; } s = prj->p[1] + cosdeg (theta); if (s == 0.0) { return 2; } *x = prj->w[0]*phi; *y = prj->w[2]*sindeg (theta)/s; return 0; } /*--------------------------------------------------------------------------*/ int cyprev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double eta; if (prj->flag != CYP) { if (cypset(prj)) return 1; } *phi = x*prj->w[1]; eta = y*prj->w[3]; *theta = atan2deg (eta,1.0) + asindeg (eta*prj->p[1]/sqrt(eta*eta+1.0)); return 0; } /*============================================================================ * CEA: cylindrical equal area projection. * * Given: * prj->p[1] Square of the cosine of the latitude at which the * projection is conformal, lambda. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "CEA" * prj->flag CEA * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/180) * prj->w[1] (180/pi)/r0 * prj->w[2] r0/lambda * prj->w[3] lambda/r0 * prj->prjfwd Pointer to ceafwd(). * prj->prjrev Pointer to cearev(). *===========================================================================*/ int ceaset(prj) struct prjprm *prj; { strcpy(prj->code, "CEA"); prj->flag = CEA; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 1.0; prj->w[1] = 1.0; if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) { return 1; } prj->w[2] = prj->r0/prj->p[1]; prj->w[3] = prj->p[1]/prj->r0; } else { prj->w[0] = prj->r0*D2R; prj->w[1] = R2D/prj->r0; if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) { return 1; } prj->w[2] = prj->r0/prj->p[1]; prj->w[3] = prj->p[1]/prj->r0; } prj->prjfwd = ceafwd; prj->prjrev = cearev; return 0; } /*--------------------------------------------------------------------------*/ int ceafwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { if (prj->flag != CEA) { if (ceaset(prj)) return 1; } *x = prj->w[0]*phi; *y = prj->w[2]*sindeg (theta); return 0; } /*--------------------------------------------------------------------------*/ int cearev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double s; const double tol = 1.0e-13; if (prj->flag != CEA) { if (ceaset(prj)) return 1; } s = y*prj->w[3]; if (fabs(s) > 1.0) { if (fabs(s) > 1.0+tol) { return 2; } s = copysgn (1.0,s); } *phi = x*prj->w[1]; *theta = asindeg (s); return 0; } /*============================================================================ * CAR: Cartesian projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "CAR" * prj->flag CAR * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/180) * prj->w[1] (180/pi)/r0 * prj->prjfwd Pointer to carfwd(). * prj->prjrev Pointer to carrev(). *===========================================================================*/ int carset(prj) struct prjprm *prj; { strcpy(prj->code, "CAR"); prj->flag = CAR; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 1.0; prj->w[1] = 1.0; } else { prj->w[0] = prj->r0*D2R; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = carfwd; prj->prjrev = carrev; return 0; } /*--------------------------------------------------------------------------*/ int carfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { if (prj->flag != CAR) { if (carset(prj)) return 1; } *x = prj->w[0]*phi; *y = prj->w[0]*theta; return 0; } /*--------------------------------------------------------------------------*/ int carrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { if (prj->flag != CAR) { if (carset(prj)) return 1; } *phi = prj->w[1]*x; *theta = prj->w[1]*y; return 0; } /*============================================================================ * MER: Mercator's projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "MER" * prj->flag MER * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/180) * prj->w[1] (180/pi)/r0 * prj->prjfwd Pointer to merfwd(). * prj->prjrev Pointer to merrev(). *===========================================================================*/ int merset(prj) struct prjprm *prj; { strcpy(prj->code, "MER"); prj->flag = MER; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 1.0; prj->w[1] = 1.0; } else { prj->w[0] = prj->r0*D2R; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = merfwd; prj->prjrev = merrev; return 0; } /*--------------------------------------------------------------------------*/ int merfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { if (prj->flag != MER) { if (merset(prj)) return 1; } if (theta <= -90.0 || theta >= 90.0) { return 2; } *x = prj->w[0]*phi; *y = prj->r0*log(tandeg ((90.0+theta)/2.0)); return 0; } /*--------------------------------------------------------------------------*/ int merrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { if (prj->flag != MER) { if (merset(prj)) return 1; } *phi = x*prj->w[1]; *theta = 2.0*atandeg (exp(y/prj->r0)) - 90.0; return 0; } /*============================================================================ * SFL: Sanson-Flamsteed ("global sinusoid") projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "SFL" * prj->flag SFL * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/180) * prj->w[1] (180/pi)/r0 * prj->prjfwd Pointer to sflfwd(). * prj->prjrev Pointer to sflrev(). *===========================================================================*/ int sflset(prj) struct prjprm *prj; { strcpy(prj->code, "SFL"); prj->flag = SFL; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 1.0; prj->w[1] = 1.0; } else { prj->w[0] = prj->r0*D2R; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = sflfwd; prj->prjrev = sflrev; return 0; } /*--------------------------------------------------------------------------*/ int sflfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { if (prj->flag != SFL) { if (sflset(prj)) return 1; } *x = prj->w[0]*phi*cosdeg (theta); *y = prj->w[0]*theta; return 0; } /*--------------------------------------------------------------------------*/ int sflrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double w; if (prj->flag != SFL) { if (sflset(prj)) return 1; } w = cos(y/prj->r0); if (w == 0.0) { *phi = 0.0; } else { *phi = x*prj->w[1]/cos(y/prj->r0); } *theta = y*prj->w[1]; return 0; } /*============================================================================ * PAR: parabolic projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "PAR" * prj->flag PAR * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/180) * prj->w[1] (180/pi)/r0 * prj->w[2] pi*r0 * prj->w[3] 1/(pi*r0) * prj->prjfwd Pointer to parfwd(). * prj->prjrev Pointer to parrev(). *===========================================================================*/ int parset(prj) struct prjprm *prj; { strcpy(prj->code, "PAR"); prj->flag = PAR; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 1.0; prj->w[1] = 1.0; prj->w[2] = 180.0; prj->w[3] = 1.0/prj->w[2]; } else { prj->w[0] = prj->r0*D2R; prj->w[1] = 1.0/prj->w[0]; prj->w[2] = PI*prj->r0; prj->w[3] = 1.0/prj->w[2]; } prj->prjfwd = parfwd; prj->prjrev = parrev; return 0; } /*--------------------------------------------------------------------------*/ int parfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double s; if (prj->flag != PAR) { if (parset(prj)) return 1; } s = sindeg (theta/3.0); *x = prj->w[0]*phi*(1.0 - 4.0*s*s); *y = prj->w[2]*s; return 0; } /*--------------------------------------------------------------------------*/ int parrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double s, t; if (prj->flag != PAR) { if (parset(prj)) return 1; } s = y*prj->w[3]; if (s > 1.0 || s < -1.0) { return 2; } t = 1.0 - 4.0*s*s; if (t == 0.0) { if (x == 0.0) { *phi = 0.0; } else { return 2; } } else { *phi = prj->w[1]*x/t; } *theta = 3.0*asindeg (s); return 0; } /*============================================================================ * MOL: Mollweide's projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "MOL" * prj->flag MOL * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] sqrt(2)*r0 * prj->w[1] sqrt(2)*r0/90 * prj->w[2] 1/(sqrt(2)*r0) * prj->w[3] 90/r0 * prj->prjfwd Pointer to molfwd(). * prj->prjrev Pointer to molrev(). *===========================================================================*/ int molset(prj) struct prjprm *prj; { strcpy(prj->code, "MOL"); prj->flag = MOL; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) prj->r0 = R2D; prj->w[0] = SQRT2*prj->r0; prj->w[1] = prj->w[0]/90.0; prj->w[2] = 1.0/prj->w[0]; prj->w[3] = 90.0/prj->r0; prj->w[4] = 2.0/PI; prj->prjfwd = molfwd; prj->prjrev = molrev; return 0; } /*--------------------------------------------------------------------------*/ int molfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { int j; double gamma, resid, u, v, v0, v1; const double tol = 1.0e-13; if (prj->flag != MOL) { if (molset(prj)) return 1; } if (fabs(theta) == 90.0) { *x = 0.0; *y = copysgn (prj->w[0],theta); } else if (theta == 0.0) { *x = prj->w[1]*phi; *y = 0.0; } else { u = PI*sindeg (theta); v0 = -PI; v1 = PI; v = u; for (j = 0; j < 100; j++) { resid = (v - u) + sin(v); if (resid < 0.0) { if (resid > -tol) break; v0 = v; } else { if (resid < tol) break; v1 = v; } v = (v0 + v1)/2.0; } gamma = v/2.0; *x = prj->w[1]*phi*cos(gamma); *y = prj->w[0]*sin(gamma); } return 0; } /*--------------------------------------------------------------------------*/ int molrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double s, y0, z; const double tol = 1.0e-12; if (prj->flag != MOL) { if (molset(prj)) return 1; } y0 = y/prj->r0; s = 2.0 - y0*y0; if (s <= tol) { if (s < -tol) { return 2; } s = 0.0; if (fabs(x) > tol) { return 2; } *phi = 0.0; } else { s = sqrt(s); *phi = prj->w[3]*x/s; } z = y*prj->w[2]; if (fabs(z) > 1.0) { if (fabs(z) > 1.0+tol) { return 2; } z = copysgn (1.0,z) + y0*s/PI; } else { z = asin(z)*prj->w[4] + y0*s/PI; } if (fabs(z) > 1.0) { if (fabs(z) > 1.0+tol) { return 2; } z = copysgn (1.0,z); } *theta = asindeg (z); return 0; } /*============================================================================ * AIT: Hammer-Aitoff projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "AIT" * prj->flag AIT * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] 2*r0**2 * prj->w[1] 1/(2*r0)**2 * prj->w[2] 1/(4*r0)**2 * prj->w[3] 1/(2*r0) * prj->prjfwd Pointer to aitfwd(). * prj->prjrev Pointer to aitrev(). *===========================================================================*/ int aitset(prj) struct prjprm *prj; { strcpy(prj->code, "AIT"); prj->flag = AIT; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) prj->r0 = R2D; prj->w[0] = 2.0*prj->r0*prj->r0; prj->w[1] = 1.0/(2.0*prj->w[0]); prj->w[2] = prj->w[1]/4.0; prj->w[3] = 1.0/(2.0*prj->r0); prj->prjfwd = aitfwd; prj->prjrev = aitrev; return 0; } /*--------------------------------------------------------------------------*/ int aitfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double cthe, w; if (prj->flag != AIT) { if (aitset(prj)) return 1; } cthe = cosdeg (theta); w = sqrt(prj->w[0]/(1.0 + cthe*cosdeg (phi/2.0))); *x = 2.0*w*cthe*sindeg (phi/2.0); *y = w*sindeg (theta); return 0; } /*--------------------------------------------------------------------------*/ int aitrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double s, u, xp, yp, z; const double tol = 1.0e-13; if (prj->flag != AIT) { if (aitset(prj)) return 1; } u = 1.0 - x*x*prj->w[2] - y*y*prj->w[1]; if (u < 0.0) { if (u < -tol) { return 2; } u = 0.0; } z = sqrt(u); s = z*y/prj->r0; if (fabs(s) > 1.0) { if (fabs(s) > 1.0+tol) { return 2; } s = copysgn (1.0,s); } xp = 2.0*z*z - 1.0; yp = z*x*prj->w[3]; if (xp == 0.0 && yp == 0.0) { *phi = 0.0; } else { *phi = 2.0*atan2deg (yp, xp); } *theta = asindeg (s); return 0; } /*============================================================================ * COP: conic perspective projection. * * Given: * prj->p[1] sigma = (theta2+theta1)/2 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the * latitudes of the standard parallels, in degrees. * * Given and/or returned: * prj->flag COP, or -COP if prj->flag is given < 0. * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "COP" * prj->phi0 0.0 * prj->theta0 sigma * prj->w[0] C = sin(sigma) * prj->w[1] 1/C * prj->w[2] Y0 = r0*cos(delta)*cot(sigma) * prj->w[3] r0*cos(delta) * prj->w[4] 1/(r0*cos(delta) * prj->w[5] cot(sigma) * prj->prjfwd Pointer to copfwd(). * prj->prjrev Pointer to coprev(). *===========================================================================*/ int copset(prj) struct prjprm *prj; { strcpy(prj->code, "COP"); prj->flag = copysgni (COP, prj->flag); prj->phi0 = 0.0; prj->theta0 = prj->p[1]; if (prj->r0 == 0.0) prj->r0 = R2D; prj->w[0] = sindeg (prj->p[1]); if (prj->w[0] == 0.0) { return 1; } prj->w[1] = 1.0/prj->w[0]; prj->w[3] = prj->r0*cosdeg (prj->p[2]); if (prj->w[3] == 0.0) { return 1; } prj->w[4] = 1.0/prj->w[3]; prj->w[5] = 1.0/tandeg (prj->p[1]); prj->w[2] = prj->w[3]*prj->w[5]; prj->prjfwd = copfwd; prj->prjrev = coprev; return 0; } /*--------------------------------------------------------------------------*/ int copfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, r, s, t; if (abs(prj->flag) != COP) { if (copset(prj)) return 1; } t = theta - prj->p[1]; s = cosdeg (t); if (s == 0.0) { return 2; } a = prj->w[0]*phi; r = prj->w[2] - prj->w[3]*sindeg (t)/s; *x = r*sindeg (a); *y = prj->w[2] - r*cosdeg (a); if (prj->flag > 0 && r*prj->w[0] < 0.0) { return 2; } return 0; } /*--------------------------------------------------------------------------*/ int coprev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double a, dy, r; if (abs(prj->flag) != COP) { if (copset(prj)) return 1; } dy = prj->w[2] - y; r = sqrt(x*x + dy*dy); if (prj->p[1] < 0.0) r = -r; if (r == 0.0) { a = 0.0; } else { a = atan2deg (x/r, dy/r); } *phi = a*prj->w[1]; *theta = prj->p[1] + atandeg (prj->w[5] - r*prj->w[4]); return 0; } /*============================================================================ * COE: conic equal area projection. * * Given: * prj->p[1] sigma = (theta2+theta1)/2 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the * latitudes of the standard parallels, in degrees. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "COE" * prj->flag COE * prj->phi0 0.0 * prj->theta0 sigma * prj->w[0] C = (sin(theta1) + sin(theta2))/2 * prj->w[1] 1/C * prj->w[2] Y0 = chi*sqrt(psi - 2C*sindeg (sigma)) * prj->w[3] chi = r0/C * prj->w[4] psi = 1 + sin(theta1)*sin(theta2) * prj->w[5] 2C * prj->w[6] (1 + sin(theta1)*sin(theta2))*(r0/C)**2 * prj->w[7] C/(2*r0**2) * prj->w[8] chi*sqrt(psi + 2C) * prj->prjfwd Pointer to coefwd(). * prj->prjrev Pointer to coerev(). *===========================================================================*/ int coeset(prj) struct prjprm *prj; { double theta1, theta2; strcpy(prj->code, "COE"); prj->flag = COE; prj->phi0 = 0.0; prj->theta0 = prj->p[1]; if (prj->r0 == 0.0) prj->r0 = R2D; theta1 = prj->p[1] - prj->p[2]; theta2 = prj->p[1] + prj->p[2]; prj->w[0] = (sindeg (theta1) + sindeg (theta2))/2.0; if (prj->w[0] == 0.0) { return 1; } prj->w[1] = 1.0/prj->w[0]; prj->w[3] = prj->r0/prj->w[0]; prj->w[4] = 1.0 + sindeg (theta1)*sindeg (theta2); prj->w[5] = 2.0*prj->w[0]; prj->w[6] = prj->w[3]*prj->w[3]*prj->w[4]; prj->w[7] = 1.0/(2.0*prj->r0*prj->w[3]); prj->w[8] = prj->w[3]*sqrt(prj->w[4] + prj->w[5]); prj->w[2] = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*sindeg (prj->p[1])); prj->prjfwd = coefwd; prj->prjrev = coerev; return 0; } /*--------------------------------------------------------------------------*/ int coefwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, r; if (prj->flag != COE) { if (coeset(prj)) return 1; } a = phi*prj->w[0]; if (theta == -90.0) { r = prj->w[8]; } else { r = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*sindeg (theta)); } *x = r*sindeg (a); *y = prj->w[2] - r*cosdeg (a); return 0; } /*--------------------------------------------------------------------------*/ int coerev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double a, dy, r, w; const double tol = 1.0e-12; if (prj->flag != COE) { if (coeset(prj)) return 1; } dy = prj->w[2] - y; r = sqrt(x*x + dy*dy); if (prj->p[1] < 0.0) r = -r; if (r == 0.0) { a = 0.0; } else { a = atan2deg (x/r, dy/r); } *phi = a*prj->w[1]; if (fabs(r - prj->w[8]) < tol) { *theta = -90.0; } else { w = (prj->w[6] - r*r)*prj->w[7]; if (fabs(w) > 1.0) { if (fabs(w-1.0) < tol) { *theta = 90.0; } else if (fabs(w+1.0) < tol) { *theta = -90.0; } else { return 2; } } else { *theta = asindeg (w); } } return 0; } /*============================================================================ * COD: conic equidistant projection. * * Given: * prj->p[1] sigma = (theta2+theta1)/2 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the * latitudes of the standard parallels, in degrees. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "COD" * prj->flag COD * prj->phi0 0.0 * prj->theta0 sigma * prj->w[0] C = r0*sin(sigma)*sin(delta)/delta * prj->w[1] 1/C * prj->w[2] Y0 = delta*cot(delta)*cot(sigma) * prj->w[3] Y0 + sigma * prj->prjfwd Pointer to codfwd(). * prj->prjrev Pointer to codrev(). *===========================================================================*/ int codset(prj) struct prjprm *prj; { strcpy(prj->code, "COD"); prj->flag = COD; prj->phi0 = 0.0; prj->theta0 = prj->p[1]; if (prj->r0 == 0.0) prj->r0 = R2D; if (prj->p[2] == 0.0) { prj->w[0] = prj->r0*sindeg (prj->p[1])*D2R; } else { prj->w[0] = prj->r0*sindeg (prj->p[1])*sindeg (prj->p[2])/prj->p[2]; } if (prj->w[0] == 0.0) { return 1; } prj->w[1] = 1.0/prj->w[0]; prj->w[2] = prj->r0*cosdeg (prj->p[2])*cosdeg (prj->p[1])/prj->w[0]; prj->w[3] = prj->w[2] + prj->p[1]; prj->prjfwd = codfwd; prj->prjrev = codrev; return 0; } /*--------------------------------------------------------------------------*/ int codfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, r; if (prj->flag != COD) { if (codset(prj)) return 1; } a = prj->w[0]*phi; r = prj->w[3] - theta; *x = r*sindeg (a); *y = prj->w[2] - r*cosdeg (a); return 0; } /*--------------------------------------------------------------------------*/ int codrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double a, dy, r; if (prj->flag != COD) { if (codset(prj)) return 1; } dy = prj->w[2] - y; r = sqrt(x*x + dy*dy); if (prj->p[1] < 0.0) r = -r; if (r == 0.0) { a = 0.0; } else { a = atan2deg (x/r, dy/r); } *phi = a*prj->w[1]; *theta = prj->w[3] - r; return 0; } /*============================================================================ * COO: conic orthomorphic projection. * * Given: * prj->p[1] sigma = (theta2+theta1)/2 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the * latitudes of the standard parallels, in degrees. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "COO" * prj->flag COO * prj->phi0 0.0 * prj->theta0 sigma * prj->w[0] C = ln(cos(theta2)/cos(theta1))/ln(tan(tau2)/tan(tau1)) * where tau1 = (90 - theta1)/2 * tau2 = (90 - theta2)/2 * prj->w[1] 1/C * prj->w[2] Y0 = psi*tan((90-sigma)/2)**C * prj->w[3] psi = (r0*cos(theta1)/C)/tan(tau1)**C * prj->w[4] 1/psi * prj->prjfwd Pointer to coofwd(). * prj->prjrev Pointer to coorev(). *===========================================================================*/ int cooset(prj) struct prjprm *prj; { double cos1, cos2, tan1, tan2, theta1, theta2; strcpy(prj->code, "COO"); prj->flag = COO; prj->phi0 = 0.0; prj->theta0 = prj->p[1]; if (prj->r0 == 0.0) prj->r0 = R2D; theta1 = prj->p[1] - prj->p[2]; theta2 = prj->p[1] + prj->p[2]; tan1 = tandeg ((90.0 - theta1)/2.0); cos1 = cosdeg (theta1); if (theta1 == theta2) { prj->w[0] = sindeg (theta1); } else { tan2 = tandeg ((90.0 - theta2)/2.0); cos2 = cosdeg (theta2); prj->w[0] = log(cos2/cos1)/log(tan2/tan1); } if (prj->w[0] == 0.0) { return 1; } prj->w[1] = 1.0/prj->w[0]; prj->w[3] = prj->r0*(cos1/prj->w[0])/pow(tan1,prj->w[0]); if (prj->w[3] == 0.0) { return 1; } prj->w[2] = prj->w[3]*pow(tandeg ((90.0 - prj->p[1])/2.0),prj->w[0]); prj->w[4] = 1.0/prj->w[3]; prj->prjfwd = coofwd; prj->prjrev = coorev; return 0; } /*--------------------------------------------------------------------------*/ int coofwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, r; if (prj->flag != COO) { if (cooset(prj)) return 1; } a = prj->w[0]*phi; if (theta == -90.0) { if (prj->w[0] < 0.0) { r = 0.0; } else { return 2; } } else { r = prj->w[3]*pow(tandeg ((90.0 - theta)/2.0),prj->w[0]); } *x = r*sindeg (a); *y = prj->w[2] - r*cosdeg (a); return 0; } /*--------------------------------------------------------------------------*/ int coorev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double a, dy, r; if (prj->flag != COO) { if (cooset(prj)) return 1; } dy = prj->w[2] - y; r = sqrt(x*x + dy*dy); if (prj->p[1] < 0.0) r = -r; if (r == 0.0) { a = 0.0; } else { a = atan2deg (x/r, dy/r); } *phi = a*prj->w[1]; if (r == 0.0) { if (prj->w[0] < 0.0) { *theta = -90.0; } else { return 2; } } else { *theta = 90.0 - 2.0*atandeg (pow(r*prj->w[4],prj->w[1])); } return 0; } /*============================================================================ * BON: Bonne's projection. * * Given: * prj->p[1] Bonne conformal latitude, theta1, in degrees. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "BON" * prj->flag BON * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[1] r0*pi/180 * prj->w[2] Y0 = r0*(cot(theta1) + theta1*pi/180) * prj->prjfwd Pointer to bonfwd(). * prj->prjrev Pointer to bonrev(). *===========================================================================*/ int bonset(prj) struct prjprm *prj; { strcpy(prj->code, "BON"); prj->flag = BON; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[1] = 1.0; prj->w[2] = prj->r0*cosdeg (prj->p[1])/sindeg (prj->p[1]) + prj->p[1]; } else { prj->w[1] = prj->r0*D2R; prj->w[2] = prj->r0*(cosdeg (prj->p[1])/sindeg (prj->p[1]) + prj->p[1]*D2R); } prj->prjfwd = bonfwd; prj->prjrev = bonrev; return 0; } /*--------------------------------------------------------------------------*/ int bonfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, r; if (prj->p[1] == 0.0) { /* Sanson-Flamsteed. */ return sflfwd(phi, theta, prj, x, y); } if (prj->flag != BON) { if (bonset(prj)) return 1; } r = prj->w[2] - theta*prj->w[1]; a = prj->r0*phi*cosdeg (theta)/r; *x = r*sindeg (a); *y = prj->w[2] - r*cosdeg (a); return 0; } /*--------------------------------------------------------------------------*/ int bonrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double a, cthe, dy, r; if (prj->p[1] == 0.0) { /* Sanson-Flamsteed. */ return sflrev(x, y, prj, phi, theta); } if (prj->flag != BON) { if (bonset(prj)) return 1; } dy = prj->w[2] - y; r = sqrt(x*x + dy*dy); if (prj->p[1] < 0.0) r = -r; if (r == 0.0) { a = 0.0; } else { a = atan2deg (x/r, dy/r); } *theta = (prj->w[2] - r)/prj->w[1]; cthe = cosdeg (*theta); if (cthe == 0.0) { *phi = 0.0; } else { *phi = a*(r/prj->r0)/cthe; } return 0; } /*============================================================================ * PCO: polyconic projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "PCO" * prj->flag PCO * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/180) * prj->w[1] 1/r0 * prj->w[2] 2*r0 * prj->prjfwd Pointer to pcofwd(). * prj->prjrev Pointer to pcorev(). *===========================================================================*/ int pcoset(prj) struct prjprm *prj; { strcpy(prj->code, "PCO"); prj->flag = PCO; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 1.0; prj->w[1] = 1.0; prj->w[2] = 360.0/PI; } else { prj->w[0] = prj->r0*D2R; prj->w[1] = 1.0/prj->w[0]; prj->w[2] = 2.0*prj->r0; } prj->prjfwd = pcofwd; prj->prjrev = pcorev; return 0; } /*--------------------------------------------------------------------------*/ int pcofwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { double a, cthe, cotthe, sthe; if (prj->flag != PCO) { if (pcoset(prj)) return 1; } cthe = cosdeg (theta); sthe = sindeg (theta); a = phi*sthe; if (sthe == 0.0) { *x = prj->w[0]*phi; *y = 0.0; } else { cotthe = cthe/sthe; *x = prj->r0*cotthe*sindeg (a); *y = prj->r0*(cotthe*(1.0 - cosdeg (a)) + theta*D2R); } return 0; } /*--------------------------------------------------------------------------*/ int pcorev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { int j; double f, fneg, fpos, lambda, tanthe, theneg, thepos, w, xp, xx, ymthe, yp; const double tol = 1.0e-12; if (prj->flag != PCO) { if (pcoset(prj)) return 1; } w = fabs(y*prj->w[1]); if (w < tol) { *phi = x*prj->w[1]; *theta = 0.0; } else if (fabs(w-90.0) < tol) { *phi = 0.0; *theta = copysgni (90.0,y); } else { /* Iterative solution using weighted division of the interval. */ if (y > 0.0) { thepos = 90.0; } else { thepos = -90.0; } theneg = 0.0; xx = x*x; ymthe = y - prj->w[0]*thepos; fpos = xx + ymthe*ymthe; fneg = -999.0; for (j = 0; j < 64; j++) { if (fneg < -100.0) { /* Equal division of the interval. */ *theta = (thepos+theneg)/2.0; } else { /* Weighted division of the interval. */ lambda = fpos/(fpos-fneg); if (lambda < 0.1) { lambda = 0.1; } else if (lambda > 0.9) { lambda = 0.9; } *theta = thepos - lambda*(thepos-theneg); } /* Compute the residue. */ ymthe = y - prj->w[0]*(*theta); tanthe = tandeg (*theta); f = xx + ymthe*(ymthe - prj->w[2]/tanthe); /* Check for convergence. */ if (fabs(f) < tol) break; if (fabs(thepos-theneg) < tol) break; /* Redefine the interval. */ if (f > 0.0) { thepos = *theta; fpos = f; } else { theneg = *theta; fneg = f; } } xp = prj->r0 - ymthe*tanthe; yp = x*tanthe; if (xp == 0.0 && yp == 0.0) { *phi = 0.0; } else { *phi = atan2deg (yp, xp)/sindeg (*theta); } } return 0; } /*============================================================================ * TSC: tangential spherical cube projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "TSC" * prj->flag TSC * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/4) * prj->w[1] (4/pi)/r0 * prj->prjfwd Pointer to tscfwd(). * prj->prjrev Pointer to tscrev(). *===========================================================================*/ int tscset(prj) struct prjprm *prj; { strcpy(prj->code, "TSC"); prj->flag = TSC; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 45.0; prj->w[1] = 1.0/45.0; } else { prj->w[0] = prj->r0*PI/4.0; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = tscfwd; prj->prjrev = tscrev; return 0; } /*--------------------------------------------------------------------------*/ int tscfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { int face; double cthe, l, m, n, rho; double x0 = 0.0; double y0 = 0.0; double xf = 0.0; double yf = 0.0; const double tol = 1.0e-12; if (prj->flag != TSC) { if (tscset(prj)) return 1; } cthe = cosdeg (theta); l = cthe*cosdeg (phi); m = cthe*sindeg (phi); n = sindeg (theta); face = 0; rho = n; if (l > rho) { face = 1; rho = l; } if (m > rho) { face = 2; rho = m; } if (-l > rho) { face = 3; rho = -l; } if (-m > rho) { face = 4; rho = -m; } if (-n > rho) { face = 5; rho = -n; } if (face == 0) { xf = m/rho; yf = -l/rho; x0 = 0.0; y0 = 2.0; } else if (face == 1) { xf = m/rho; yf = n/rho; x0 = 0.0; y0 = 0.0; } else if (face == 2) { xf = -l/rho; yf = n/rho; x0 = 2.0; y0 = 0.0; } else if (face == 3) { xf = -m/rho; yf = n/rho; x0 = 4.0; y0 = 0.0; } else if (face == 4) { xf = l/rho; yf = n/rho; x0 = 6.0; y0 = 0.0; } else if (face == 5) { xf = m/rho; yf = l/rho; x0 = 0.0; y0 = -2.0; } if (fabs(xf) > 1.0) { if (fabs(xf) > 1.0+tol) { return 2; } xf = copysgn (1.0,xf); } if (fabs(yf) > 1.0) { if (fabs(yf) > 1.0+tol) { return 2; } yf = copysgn (1.0,yf); } *x = prj->w[0]*(xf + x0); *y = prj->w[0]*(yf + y0); return 0; } /*--------------------------------------------------------------------------*/ int tscrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { double l, m, n, xf, yf; if (prj->flag != TSC) { if (tscset(prj)) return 1; } xf = x*prj->w[1]; yf = y*prj->w[1]; /* Check bounds. */ if (fabs(xf) <= 1.0) { if (fabs(yf) > 3.0) return 2; } else { if (fabs(xf) > 7.0) return 2; if (fabs(yf) > 1.0) return 2; } /* Map negative faces to the other side. */ if (xf < -1.0) xf += 8.0; /* Determine the face. */ if (xf > 5.0) { /* face = 4 */ xf = xf - 6.0; m = -1.0/sqrt(1.0 + xf*xf + yf*yf); l = -m*xf; n = -m*yf; } else if (xf > 3.0) { /* face = 3 */ xf = xf - 4.0; l = -1.0/sqrt(1.0 + xf*xf + yf*yf); m = l*xf; n = -l*yf; } else if (xf > 1.0) { /* face = 2 */ xf = xf - 2.0; m = 1.0/sqrt(1.0 + xf*xf + yf*yf); l = -m*xf; n = m*yf; } else if (yf > 1.0) { /* face = 0 */ yf = yf - 2.0; n = 1.0/sqrt(1.0 + xf*xf + yf*yf); l = -n*yf; m = n*xf; } else if (yf < -1.0) { /* face = 5 */ yf = yf + 2.0; n = -1.0/sqrt(1.0 + xf*xf + yf*yf); l = -n*yf; m = -n*xf; } else { /* face = 1 */ l = 1.0/sqrt(1.0 + xf*xf + yf*yf); m = l*xf; n = l*yf; } if (l == 0.0 && m == 0.0) { *phi = 0.0; } else { *phi = atan2deg (m, l); } *theta = asindeg (n); return 0; } /*============================================================================ * CSC: COBE quadrilateralized spherical cube projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "CSC" * prj->flag CSC * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/4) * prj->w[1] (4/pi)/r0 * prj->prjfwd Pointer to cscfwd(). * prj->prjrev Pointer to cscrev(). *===========================================================================*/ int cscset(prj) struct prjprm *prj; { strcpy(prj->code, "CSC"); prj->flag = CSC; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 45.0; prj->w[1] = 1.0/45.0; } else { prj->w[0] = prj->r0*PI/4.0; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = cscfwd; prj->prjrev = cscrev; return 0; } /*--------------------------------------------------------------------------*/ int cscfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { int face; double cthe, eta, l, m, n, rho, xi; const float tol = 1.0e-7; float a, a2, a2b2, a4, ab, b, b2, b4, ca2, cb2; float x0 = 0.0; float y0 = 0.0; float xf = 0.0; float yf = 0.0; const float gstar = 1.37484847732; const float mm = 0.004869491981; const float gamma = -0.13161671474; const float omega1 = -0.159596235474; const float d0 = 0.0759196200467; const float d1 = -0.0217762490699; const float c00 = 0.141189631152; const float c10 = 0.0809701286525; const float c01 = -0.281528535557; const float c11 = 0.15384112876; const float c20 = -0.178251207466; const float c02 = 0.106959469314; if (prj->flag != CSC) { if (cscset(prj)) return 1; } cthe = cosdeg (theta); l = cthe*cosdeg (phi); m = cthe*sindeg (phi); n = sindeg (theta); face = 0; rho = n; if (l > rho) { face = 1; rho = l; } if (m > rho) { face = 2; rho = m; } if (-l > rho) { face = 3; rho = -l; } if (-m > rho) { face = 4; rho = -m; } if (-n > rho) { face = 5; rho = -n; } if (face == 0) { xi = m; eta = -l; x0 = 0.0; y0 = 2.0; } else if (face == 1) { xi = m; eta = n; x0 = 0.0; y0 = 0.0; } else if (face == 2) { xi = -l; eta = n; x0 = 2.0; y0 = 0.0; } else if (face == 3) { xi = -m; eta = n; x0 = 4.0; y0 = 0.0; } else if (face == 4) { xi = l; eta = n; x0 = 6.0; y0 = 0.0; } else if (face == 5) { xi = m; eta = l; x0 = 0.0; y0 = -2.0; } a = xi/rho; b = eta/rho; a2 = a*a; b2 = b*b; ca2 = 1.0 - a2; cb2 = 1.0 - b2; /* Avoid floating underflows. */ ab = fabs(a*b); a4 = (a2 > 1.0e-16) ? a2*a2 : 0.0; b4 = (b2 > 1.0e-16) ? b2*b2 : 0.0; a2b2 = (ab > 1.0e-16) ? a2*b2 : 0.0; xf = a*(a2 + ca2*(gstar + b2*(gamma*ca2 + mm*a2 + cb2*(c00 + c10*a2 + c01*b2 + c11*a2b2 + c20*a4 + c02*b4)) + a2*(omega1 - ca2*(d0 + d1*a2)))); yf = b*(b2 + cb2*(gstar + a2*(gamma*cb2 + mm*b2 + ca2*(c00 + c10*b2 + c01*a2 + c11*a2b2 + c20*b4 + c02*a4)) + b2*(omega1 - cb2*(d0 + d1*b2)))); if (fabs(xf) > 1.0) { if (fabs(xf) > 1.0+tol) { return 2; } xf = copysgn (1.0,xf); } if (fabs(yf) > 1.0) { if (fabs(yf) > 1.0+tol) { return 2; } yf = copysgn (1.0,yf); } *x = prj->w[0]*(x0 + xf); *y = prj->w[0]*(y0 + yf); return 0; } /*--------------------------------------------------------------------------*/ int cscrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { int face; double l = 0.0; double m = 0.0; double n = 0.0; float a, b, xf, xx, yf, yy, z0, z1, z2, z3, z4, z5, z6; const float p00 = -0.27292696; const float p10 = -0.07629969; const float p20 = -0.22797056; const float p30 = 0.54852384; const float p40 = -0.62930065; const float p50 = 0.25795794; const float p60 = 0.02584375; const float p01 = -0.02819452; const float p11 = -0.01471565; const float p21 = 0.48051509; const float p31 = -1.74114454; const float p41 = 1.71547508; const float p51 = -0.53022337; const float p02 = 0.27058160; const float p12 = -0.56800938; const float p22 = 0.30803317; const float p32 = 0.98938102; const float p42 = -0.83180469; const float p03 = -0.60441560; const float p13 = 1.50880086; const float p23 = -0.93678576; const float p33 = 0.08693841; const float p04 = 0.93412077; const float p14 = -1.41601920; const float p24 = 0.33887446; const float p05 = -0.63915306; const float p15 = 0.52032238; const float p06 = 0.14381585; if (prj->flag != CSC) { if (cscset(prj)) return 1; } xf = x*prj->w[1]; yf = y*prj->w[1]; /* Check bounds. */ if (fabs(xf) <= 1.0) { if (fabs(yf) > 3.0) return 2; } else { if (fabs(xf) > 7.0) return 2; if (fabs(yf) > 1.0) return 2; } /* Map negative faces to the other side. */ if (xf < -1.0) xf += 8.0; /* Determine the face. */ if (xf > 5.0) { face = 4; xf = xf - 6.0; } else if (xf > 3.0) { face = 3; xf = xf - 4.0; } else if (xf > 1.0) { face = 2; xf = xf - 2.0; } else if (yf > 1.0) { face = 0; yf = yf - 2.0; } else if (yf < -1.0) { face = 5; yf = yf + 2.0; } else { face = 1; } xx = xf*xf; yy = yf*yf; z0 = p00 + xx*(p10 + xx*(p20 + xx*(p30 + xx*(p40 + xx*(p50 + xx*(p60)))))); z1 = p01 + xx*(p11 + xx*(p21 + xx*(p31 + xx*(p41 + xx*(p51))))); z2 = p02 + xx*(p12 + xx*(p22 + xx*(p32 + xx*(p42)))); z3 = p03 + xx*(p13 + xx*(p23 + xx*(p33))); z4 = p04 + xx*(p14 + xx*(p24)); z5 = p05 + xx*(p15); z6 = p06; a = z0 + yy*(z1 + yy*(z2 + yy*(z3 + yy*(z4 + yy*(z5 + yy*z6))))); a = xf + xf*(1.0 - xx)*a; z0 = p00 + yy*(p10 + yy*(p20 + yy*(p30 + yy*(p40 + yy*(p50 + yy*(p60)))))); z1 = p01 + yy*(p11 + yy*(p21 + yy*(p31 + yy*(p41 + yy*(p51))))); z2 = p02 + yy*(p12 + yy*(p22 + yy*(p32 + yy*(p42)))); z3 = p03 + yy*(p13 + yy*(p23 + yy*(p33))); z4 = p04 + yy*(p14 + yy*(p24)); z5 = p05 + yy*(p15); z6 = p06; b = z0 + xx*(z1 + xx*(z2 + xx*(z3 + xx*(z4 + xx*(z5 + xx*z6))))); b = yf + yf*(1.0 - yy)*b; if (face == 0) { n = 1.0/sqrt(a*a + b*b + 1.0); l = -b*n; m = a*n; } else if (face == 1) { l = 1.0/sqrt(a*a + b*b + 1.0); m = a*l; n = b*l; } else if (face == 2) { m = 1.0/sqrt(a*a + b*b + 1.0); l = -a*m; n = b*m; } else if (face == 3) { l = -1.0/sqrt(a*a + b*b + 1.0); m = a*l; n = -b*l; } else if (face == 4) { m = -1.0/sqrt(a*a + b*b + 1.0); l = -a*m; n = -b*m; } else if (face == 5) { n = -1.0/sqrt(a*a + b*b + 1.0); l = -b*n; m = -a*n; } if (l == 0.0 && m == 0.0) { *phi = 0.0; } else { *phi = atan2deg (m, l); } *theta = asindeg (n); return 0; } /*============================================================================ * QSC: quadrilaterilized spherical cube projection. * * Given and/or returned: * prj->r0 r0; reset to 180/pi if 0. * * Returned: * prj->code "QSC" * prj->flag QSC * prj->phi0 0.0 * prj->theta0 0.0 * prj->w[0] r0*(pi/4) * prj->w[1] (4/pi)/r0 * prj->prjfwd Pointer to qscfwd(). * prj->prjrev Pointer to qscrev(). *===========================================================================*/ int qscset(prj) struct prjprm *prj; { strcpy(prj->code, "QSC"); prj->flag = QSC; prj->phi0 = 0.0; prj->theta0 = 0.0; if (prj->r0 == 0.0) { prj->r0 = R2D; prj->w[0] = 45.0; prj->w[1] = 1.0/45.0; } else { prj->w[0] = prj->r0*PI/4.0; prj->w[1] = 1.0/prj->w[0]; } prj->prjfwd = qscfwd; prj->prjrev = qscrev; return 0; } /*--------------------------------------------------------------------------*/ int qscfwd(phi, theta, prj, x, y) const double phi, theta; struct prjprm *prj; double *x, *y; { int face; double cthe, l, m, n, omega, p, rho, rhu, t, tau; double xi = 0.0; double eta = 0.0; double x0 = 0.0; double y0 = 0.0; double xf = 0.0; double yf = 0.0; const double tol = 1.0e-12; if (prj->flag != QSC) { if (qscset(prj)) return 1; } if (fabs(theta) == 90.0) { *x = 0.0; *y = copysgn (2.0*prj->w[0],theta); return 0; } cthe = cosdeg (theta); l = cthe*cosdeg (phi); m = cthe*sindeg (phi); n = sindeg (theta); face = 0; rho = n; if (l > rho) { face = 1; rho = l; } if (m > rho) { face = 2; rho = m; } if (-l > rho) { face = 3; rho = -l; } if (-m > rho) { face = 4; rho = -m; } if (-n > rho) { face = 5; rho = -n; } rhu = 1.0 - rho; if (face == 0) { xi = m; eta = -l; if (rhu < 1.0e-8) { /* Small angle formula. */ t = (90.0 - theta)*D2R; rhu = t*t/2.0; } x0 = 0.0; y0 = 2.0; } else if (face == 1) { xi = m; eta = n; if (rhu < 1.0e-8) { /* Small angle formula. */ t = theta*D2R; p = fmod(phi,360.0); if (p < -180.0) p += 360.0; if (p > 180.0) p -= 360.0; p *= D2R; rhu = (p*p + t*t)/2.0; } x0 = 0.0; y0 = 0.0; } else if (face == 2) { xi = -l; eta = n; if (rhu < 1.0e-8) { /* Small angle formula. */ t = theta*D2R; p = fmod(phi,360.0); if (p < -180.0) p += 360.0; p = (90.0 - p)*D2R; rhu = (p*p + t*t)/2.0; } x0 = 2.0; y0 = 0.0; } else if (face == 3) { xi = -m; eta = n; if (rhu < 1.0e-8) { /* Small angle formula. */ t = theta*D2R; p = fmod(phi,360.0); if (p < 0.0) p += 360.0; p = (180.0 - p)*D2R; rhu = (p*p + t*t)/2.0; } x0 = 4.0; y0 = 0.0; } else if (face == 4) { xi = l; eta = n; if (rhu < 1.0e-8) { /* Small angle formula. */ t = theta*D2R; p = fmod(phi,360.0); if (p > 180.0) p -= 360.0; p *= (90.0 + p)*D2R; rhu = (p*p + t*t)/2.0; } x0 = 6; y0 = 0.0; } else if (face == 5) { xi = m; eta = l; if (rhu < 1.0e-8) { /* Small angle formula. */ t = (90.0 + theta)*D2R; rhu = t*t/2.0; } x0 = 0.0; y0 = -2; } if (xi == 0.0 && eta == 0.0) { xf = 0.0; yf = 0.0; } else if (-xi >= fabs(eta)) { omega = eta/xi; tau = 1.0 + omega*omega; xf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+tau))); yf = (xf/15.0)*(atandeg (omega) - asindeg (omega/sqrt(tau+tau))); } else if (xi >= fabs(eta)) { omega = eta/xi; tau = 1.0 + omega*omega; xf = sqrt(rhu/(1.0-1.0/sqrt(1.0+tau))); yf = (xf/15.0)*(atandeg (omega) - asindeg (omega/sqrt(tau+tau))); } else if (-eta > fabs(xi)) { omega = xi/eta; tau = 1.0 + omega*omega; yf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+tau))); xf = (yf/15.0)*(atandeg (omega) - asindeg (omega/sqrt(tau+tau))); } else if (eta > fabs(xi)) { omega = xi/eta; tau = 1.0 + omega*omega; yf = sqrt(rhu/(1.0-1.0/sqrt(1.0+tau))); xf = (yf/15.0)*(atandeg (omega) - asindeg (omega/sqrt(tau+tau))); } if (fabs(xf) > 1.0) { if (fabs(xf) > 1.0+tol) { return 2; } xf = copysgn (1.0,xf); } if (fabs(yf) > 1.0) { if (fabs(yf) > 1.0+tol) { return 2; } yf = copysgn (1.0,yf); } *x = prj->w[0]*(xf + x0); *y = prj->w[0]*(yf + y0); return 0; } /*--------------------------------------------------------------------------*/ int qscrev(x, y, prj, phi, theta) const double x, y; struct prjprm *prj; double *phi, *theta; { int direct, face; double omega, rho, rhu, tau, xf, yf, w; double l = 0.0; double m = 0.0; double n = 0.0; const double tol = 1.0e-12; if (prj->flag != QSC) { if (qscset(prj)) return 1; } xf = x*prj->w[1]; yf = y*prj->w[1]; /* Check bounds. */ if (fabs(xf) <= 1.0) { if (fabs(yf) > 3.0) return 2; } else { if (fabs(xf) > 7.0) return 2; if (fabs(yf) > 1.0) return 2; } /* Map negative faces to the other side. */ if (xf < -1.0) xf += 8.0; /* Determine the face. */ if (xf > 5.0) { face = 4; xf = xf - 6.0; } else if (xf > 3.0) { face = 3; xf = xf - 4.0; } else if (xf > 1.0) { face = 2; xf = xf - 2.0; } else if (yf > 1.0) { face = 0; yf = yf - 2.0; } else if (yf < -1.0) { face = 5; yf = yf + 2.0; } else { face = 1; } direct = (fabs(xf) > fabs(yf)); if (direct) { if (xf == 0.0) { omega = 0.0; tau = 1.0; rho = 1.0; rhu = 0.0; } else { w = 15.0*yf/xf; omega = sindeg (w)/(cosdeg (w) - SQRT2INV); tau = 1.0 + omega*omega; rhu = xf*xf*(1.0 - 1.0/sqrt(1.0 + tau)); rho = 1.0 - rhu; } } else { if (yf == 0.0) { omega = 0.0; tau = 1.0; rho = 1.0; rhu = 0.0; } else { w = 15.0*xf/yf; omega = sindeg (w)/(cosdeg (w) - SQRT2INV); tau = 1.0 + omega*omega; rhu = yf*yf*(1.0 - 1.0/sqrt(1.0 + tau)); rho = 1.0 - rhu; } } if (rho < -1.0) { if (rho < -1.0-tol) { return 2; } rho = -1.0; rhu = 2.0; w = 0.0; } else { w = sqrt(rhu*(2.0-rhu)/tau); } if (face == 0) { n = rho; if (direct) { m = w; if (xf < 0.0) m = -m; l = -m*omega; } else { l = w; if (yf > 0.0) l = -l; m = -l*omega; } } else if (face == 1) { l = rho; if (direct) { m = w; if (xf < 0.0) m = -m; n = m*omega; } else { n = w; if (yf < 0.0) n = -n; m = n*omega; } } else if (face == 2) { m = rho; if (direct) { l = w; if (xf > 0.0) l = -l; n = -l*omega; } else { n = w; if (yf < 0.0) n = -n; l = -n*omega; } } else if (face == 3) { l = -rho; if (direct) { m = w; if (xf > 0.0) m = -m; n = -m*omega; } else { n = w; if (yf < 0.0) n = -n; m = -n*omega; } } else if (face == 4) { m = -rho; if (direct) { l = w; if (xf < 0.0) l = -l; n = l*omega; } else { n = w; if (yf < 0.0) n = -n; l = n*omega; } } else if (face == 5) { n = -rho; if (direct) { m = w; if (xf < 0.0) m = -m; l = m*omega; } else { l = w; if (yf < 0.0) l = -l; m = l*omega; } } if (l == 0.0 && m == 0.0) { *phi = 0.0; } else { *phi = atan2deg (m, l); } *theta = asindeg (n); return 0; } /* This routine comes from E. Bertin sextractor-2.8.6 */ int raw_to_pv(struct prjprm *prj, double x, double y, double *xo, double *yo) { int k; double *a,*b, r,r3,r5,r7,xy,x2,x3,x4,x5,x6,x7,y2,y3,y4,y5,y6,y7,xp,yp; k=prj->npv; a = prj->ppv+MAXPV; /* Latitude comes first for compatibility */ b = prj->ppv; /* Longitude */ xp = *(a++); xp += *(a++)*x; yp = *(b++); yp += *(b++)*y; if (!--k) goto poly_end; xp += *(a++)*y; yp += *(b++)*x; if (!--k) goto poly_end; r = sqrt(x*x + y*y); xp += *(a++)*r; yp += *(b++)*r; if (!--k) goto poly_end; xp += *(a++)*(x2=x*x); yp += *(b++)*(y2=y*y); if (!--k) goto poly_end; xp += *(a++)*(xy=x*y); yp += *(b++)*xy; if (!--k) goto poly_end; xp += *(a++)*y2; yp += *(b++)*x2; if (!--k) goto poly_end; xp += *(a++)*(x3=x*x2); yp += *(b++)*(y3=y*y2); if (!--k) goto poly_end; xp += *(a++)*x2*y; yp += *(b++)*y2*x; if (!--k) goto poly_end; xp += *(a++)*x*y2; yp += *(b++)*y*x2; if (!--k) goto poly_end; xp += *(a++)*y3; yp += *(b++)*x3; if (!--k) goto poly_end; xp += *(a++)*(r3=r*r*r); yp += *(b++)*r3; if (!--k) goto poly_end; xp += *(a++)*(x4=x2*x2); yp += *(b++)*(y4=y2*y2); if (!--k) goto poly_end; xp += *(a++)*x3*y; yp += *(b++)*y3*x; if (!--k) goto poly_end; xp += *(a++)*x2*y2; yp += *(b++)*x2*y2; if (!--k) goto poly_end; xp += *(a++)*x*y3; yp += *(b++)*y*x3; if (!--k) goto poly_end; xp += *(a++)*y4; yp += *(b++)*x4; if (!--k) goto poly_end; xp += *(a++)*(x5=x4*x); yp += *(b++)*(y5=y4*y); if (!--k) goto poly_end; xp += *(a++)*x4*y; yp += *(b++)*y4*x; if (!--k) goto poly_end; xp += *(a++)*x3*y2; yp += *(b++)*y3*x2; if (!--k) goto poly_end; xp += *(a++)*x2*y3; yp += *(b++)*y2*x3; if (!--k) goto poly_end; xp += *(a++)*x*y4; yp += *(b++)*y*x4; if (!--k) goto poly_end; xp += *(a++)*y5; yp += *(b++)*x5; if (!--k) goto poly_end; xp += *(a++)*(r5=r3*r*r); yp += *(b++)*r5; if (!--k) goto poly_end; xp += *(a++)*(x6=x5*x); yp += *(b++)*(y6=y5*y); if (!--k) goto poly_end; xp += *(a++)*x5*y; yp += *(b++)*y5*x; if (!--k) goto poly_end; xp += *(a++)*x4*y2; yp += *(b++)*y4*x2; if (!--k) goto poly_end; xp += *(a++)*x3*y3; yp += *(b++)*y3*x3; if (!--k) goto poly_end; xp += *(a++)*x2*y4; yp += *(b++)*y2*x4; if (!--k) goto poly_end; xp += *(a++)*x*y5; yp += *(b++)*y*x5; if (!--k) goto poly_end; xp += *(a++)*y6; yp += *(b++)*x6; if (!--k) goto poly_end; xp += *(a++)*(x7=x6*x); yp += *(b++)*(y7=y6*y); if (!--k) goto poly_end; xp += *(a++)*x6*y; yp += *(b++)*y6*x; if (!--k) goto poly_end; xp += *(a++)*x5*y2; yp += *(b++)*y5*x2; if (!--k) goto poly_end; xp += *(a++)*x4*y3; yp += *(b++)*y4*x3; if (!--k) goto poly_end; xp += *(a++)*x3*y4; yp += *(b++)*y3*x4; if (!--k) goto poly_end; xp += *(a++)*x2*y5; yp += *(b++)*y2*x5; if (!--k) goto poly_end; xp += *(a++)*x*y6; yp += *(b++)*y*x6; if (!--k) goto poly_end; xp += *(a++)*y7; yp += *(b++)*x7; if (!--k) goto poly_end; xp += *a*(r7=r5*r*r); yp += *b*r7; poly_end: *xo = xp; *yo = yp; return 0; } /* Dec 20 1999 Doug Mink - Change cosd() and sind() to cosdeg() and sindeg() * Dec 20 1999 Doug Mink - Include wcslib.h, which includes proj.h, wcsmath.h * Dec 20 1999 Doug Mink - Define copysign only if it is not defined * Dec 20 1999 Doug Mink - tanfwd() returns error if s<=0.0, not only if s==0.0 * * Jun 2 2000 Doug Mink - include stdlib.h to get abs() * * Feb 15 2001 Doug Mink - update zearev() for WCSLIB 2.6 * Sep 19 2001 Doug Mink - Make above changes for WCSLIB 2.7 * * Mar 15 2002 Doug Mink - Make above changes for WCSLIB 2.8.2 * * Feb 3 2003 Doug Mink - Use locally defined copysgn() and copysgni(), * not copysign() * Apr 1 2003 Doug Mink - include string.h for strcpy() and strcmp() * * Mar 14 2011 Doug Mink - If no coefficients in ZPN, make ARC * Mar 14 2011 Doug Mink - Add Emmanuel Bertin's TAN polynomial from Ed Los */