/*============================================================================
*
* 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 2.20 2002/04/03 01:25:29 mcalabre Exp $
*===========================================================================*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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 = (zd1<zd2) ? zd1 : zd2;
if (zd < -tol) zd = (zd1>zd2) ? 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 = copysgn (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=0.0, l, m, n, rho, xi=0.0;
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
*/