diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2018-01-09 19:28:07 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2018-01-09 19:28:07 (GMT) |
commit | 3bfbbb90d4d6ebd44cc5eac38af24d1f78318a58 (patch) | |
tree | f278119398ae5d67a6a338705a76db420f6b8f7e /ast/proj.c | |
parent | 1332d38f2805d986ea130e43218c0d2e870b4dc1 (diff) | |
download | blt-3bfbbb90d4d6ebd44cc5eac38af24d1f78318a58.zip blt-3bfbbb90d4d6ebd44cc5eac38af24d1f78318a58.tar.gz blt-3bfbbb90d4d6ebd44cc5eac38af24d1f78318a58.tar.bz2 |
update ast 8.6.2
Diffstat (limited to 'ast/proj.c')
-rw-r--r-- | ast/proj.c | 4840 |
1 files changed, 4840 insertions, 0 deletions
diff --git a/ast/proj.c b/ast/proj.c new file mode 100644 index 0000000..8b92c21 --- /dev/null +++ b/ast/proj.c @@ -0,0 +1,4840 @@ +/*============================================================================ +* +* 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 Library 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 Library +* General Public License for more details. +* +* You should have received a copy of the GNU Library General Public License +* along with this library; if not, write to the Free Software Foundation, +* Inc., 51 Franklin Street,Fifth Floor, Boston, MA 02110-1301, 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 +* +* +*============================================================================= +* +* This version of proj.c is based on the version in wcslib-2.9, but has +* been modified in the following ways by the Starlink project (e-mail: +* ussc@star.rl.ac.uk): +* - The copysign macro is now always defined within this file +* instead of only being defined if the COPYSIGN macro has previously +* been defined. +* - Sine values which are slightly larger than 1.0 are now treated +* as 1.0 in function astCYPrev. +* - The maximum number of projection parameters has been changed from +* 10 to 100. +* - The maximum number of projection parameters is given by the +* WCSLIB_MXPAR macro (defined in proj.h) instead of being hard-wired. +* - The names of all functions and structures have been chanegd to avoid +* clashes with wcslib. This involves adding "Ast" or "ast" at the +* front and changing the capitalisation. +* - Include string.h (for strcpy and strcmp prototypes). +* - Include stdlib.h (for abs prototype). +* - Comment out declarations of npcode and pcodes variables (they +* are not needed by AST) in order to avoid clash with similar names +* in other modules imported as part of other software systems (e.g. +* SkyCat). +* - astZPNfwd: Loop from prj->n to zero, not from MAXPAR to zero. +* - astZPNfwd: Only return "2" if prj->n is larger than 2. +* - Lots of variables are initialised to null values in order to +* avoid "use of uninitialised variable" messages from compilers which +* are not clever enough to work out that the uninitialised variable is +* not in fact ever used. +* - Use dynamic rather than static memory for the parameter arrays in +* the AstPrjPrm structure.Override astGetObjSize. This is to +* reduce the in-memory size of a WcsMap. +* - HPX and XPH projections included from a more recent version of WCSLIB, +* and modified to use scalar instead of vector positions +* - The expressions for xc in astHPXrev and phic in astHPXfwd have +* been conditioned differently to the WCSLIB code in order to improve +* accuracy of the floor function for arguments very slightly below an +* integer value. + +*============================================================================= +* +* 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. +* +* astPRJset astPRJfwd astPRJrev Driver routines (see below). +* +* astAZPset astAZPfwd astAZPrev AZP: zenithal/azimuthal perspective +* astSZPset astSZPfwd astSZPrev SZP: slant zenithal perspective +* astTANset astTANfwd astTANrev TAN: gnomonic +* astSTGset astSTGfwd astSTGrev STG: stereographic +* astSINset astSINfwd astSINrev SIN: orthographic/synthesis +* astARCset astARCfwd astARCrev ARC: zenithal/azimuthal equidistant +* astZPNset astZPNfwd astZPNrev ZPN: zenithal/azimuthal polynomial +* astZEAset astZEAfwd astZEArev ZEA: zenithal/azimuthal equal area +* astAIRset astAIRfwd astAIRrev AIR: Airy +* astCYPset astCYPfwd astCYPrev CYP: cylindrical perspective +* astCEAset astCEAfwd astCEArev CEA: cylindrical equal area +* astCARset astCARfwd astCARrev CAR: Cartesian +* astMERset astMERfwd astMERrev MER: Mercator +* astSFLset astSFLfwd astSFLrev SFL: Sanson-Flamsteed +* astPARset astPARfwd astPARrev PAR: parabolic +* astMOLset astMOLfwd astMOLrev MOL: Mollweide +* astAITset astAITfwd astAITrev AIT: Hammer-Aitoff +* astCOPset astCOPfwd astCOPrev COP: conic perspective +* astCOEset astCOEfwd astCOErev COE: conic equal area +* astCODset astCODfwd astCODrev COD: conic equidistant +* astCOOset astCOOfwd astCOOrev COO: conic orthomorphic +* astBONset astBONfwd astBONrev BON: Bonne +* astPCOset astPCOfwd astPCOrev PCO: polyconic +* astTSCset astTSCfwd astTSCrev TSC: tangential spherical cube +* astCSCset astCSCfwd astCSCrev CSC: COBE quadrilateralized spherical cube +* astQSCset astQSCfwd astQSCrev QSC: quadrilateralized spherical cube +* astHPXset astHPXfwd astHPXrev HPX: HEALPix projection +* astXPHset astXPHfwd astXPHrev XPH: HEALPix polar, aka "butterfly" +* +* +* Driver routines; astPRJset(), astPRJfwd() & astPRJrev() +* ---------------------------------------------- +* A set of driver routines are available for use as a generic interface to +* the specific projection routines. The interfaces to astPRJfwd() and astPRJrev() +* are the same as those of the forward and reverse transformation routines +* for the specific projections (see below). +* +* The interface to astPRJset() differs slightly from that of the initialization +* routines for the specific projections and unlike them it must be invoked +* explicitly to use astPRJfwd() and astPRJrev(). +* +* Given: +* pcode[4] const char +* WCS projection code. +* +* Given and/or returned: +* prj AstPrjPrm* Projection parameters (see below). +* +* Function return value: +* int Error status +* 0: Success. +* +* +* Initialization routine; *set() +* ------------------------------ +* Initializes members of a AstPrjPrm data structure which hold intermediate +* values. Note that this routine need not be called directly; it will be +* invoked by astPRJfwd() and astPRJrev() if the "flag" structure member is +* anything other than a predefined magic value. +* +* Given and/or returned: +* prj AstPrjPrm* 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 AstPrjPrm* 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 AstPrjPrm* 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 AstPrjPrm struct consists of the following: +* +* int flag +* This flag must be set to zero whenever any of p[] 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[] +* Contains the projection parameters associated with the +* longitude axis. +* +* The remaining members of the AstPrjPrm 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 (*astPRJfwd)() +* int (*astPRJrev)() +* 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$ +*===========================================================================*/ + +/* Set the name of the module we are implementing. This indicates to + the header files that define class interfaces that they should make + "protected" symbols available. NB, this module is not a proper AST + class, but it defines this macro sanyway in order to get the protected + symbols defined in memory.h */ + +#include <math.h> +#include <string.h> +#include <stdlib.h> +#include "wcsmath.h" +#include "wcstrig.h" +#include "memory.h" +#include "proj.h" + +/* Following variables are not needed in AST and are commented out to + avoid name clashes with other software systems (e.g. SkyCat) which + defines them. + +int npcode = 28; +char pcodes[28][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", "HPX", "XPH"}; +*/ + +const int WCS__AZP = 101; +const int WCS__SZP = 102; +const int WCS__TAN = 103; +const int WCS__STG = 104; +const int WCS__SIN = 105; +const int WCS__ARC = 106; +const int WCS__ZPN = 107; +const int WCS__ZEA = 108; +const int WCS__AIR = 109; +const int WCS__CYP = 201; +const int WCS__CEA = 202; +const int WCS__CAR = 203; +const int WCS__MER = 204; +const int WCS__SFL = 301; +const int WCS__PAR = 302; +const int WCS__MOL = 303; +const int WCS__AIT = 401; +const int WCS__COP = 501; +const int WCS__COE = 502; +const int WCS__COD = 503; +const int WCS__COO = 504; +const int WCS__BON = 601; +const int WCS__PCO = 602; +const int WCS__TSC = 701; +const int WCS__CSC = 702; +const int WCS__QSC = 703; +const int WCS__HPX = 801; +const int WCS__XPH = 802; + +/* Map error number to error message for each function. */ +const char *astPRJset_errmsg[] = { + 0, + "Invalid projection parameters"}; + +const char *astPRJfwd_errmsg[] = { + 0, + "Invalid projection parameters", + "Invalid value of (phi,theta)"}; + +const char *astPRJrev_errmsg[] = { + 0, + "Invalid projection parameters", + "Invalid value of (x,y)"}; + + +#define copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X)) +#define icopysign(X, Y) ((Y) < 0.0 ? -abs(X) : abs(X)) + + + +/*==========================================================================*/ + +int astPRJset(pcode, prj) + +const char pcode[4]; +struct AstPrjPrm *prj; + +{ + /* Set pointers to the forward and reverse projection routines. */ + if (strcmp(pcode, "AZP") == 0) { + astAZPset(prj); + } else if (strcmp(pcode, "SZP") == 0) { + astSZPset(prj); + } else if (strcmp(pcode, "TAN") == 0) { + astTANset(prj); + } else if (strcmp(pcode, "STG") == 0) { + astSTGset(prj); + } else if (strcmp(pcode, "SIN") == 0) { + astSINset(prj); + } else if (strcmp(pcode, "ARC") == 0) { + astARCset(prj); + } else if (strcmp(pcode, "ZPN") == 0) { + astZPNset(prj); + } else if (strcmp(pcode, "ZEA") == 0) { + astZEAset(prj); + } else if (strcmp(pcode, "AIR") == 0) { + astAIRset(prj); + } else if (strcmp(pcode, "CYP") == 0) { + astCYPset(prj); + } else if (strcmp(pcode, "CEA") == 0) { + astCEAset(prj); + } else if (strcmp(pcode, "CAR") == 0) { + astCARset(prj); + } else if (strcmp(pcode, "MER") == 0) { + astMERset(prj); + } else if (strcmp(pcode, "SFL") == 0) { + astSFLset(prj); + } else if (strcmp(pcode, "PAR") == 0) { + astPARset(prj); + } else if (strcmp(pcode, "MOL") == 0) { + astMOLset(prj); + } else if (strcmp(pcode, "AIT") == 0) { + astAITset(prj); + } else if (strcmp(pcode, "COP") == 0) { + astCOPset(prj); + } else if (strcmp(pcode, "COE") == 0) { + astCOEset(prj); + } else if (strcmp(pcode, "COD") == 0) { + astCODset(prj); + } else if (strcmp(pcode, "COO") == 0) { + astCOOset(prj); + } else if (strcmp(pcode, "BON") == 0) { + astBONset(prj); + } else if (strcmp(pcode, "PCO") == 0) { + astPCOset(prj); + } else if (strcmp(pcode, "TSC") == 0) { + astTSCset(prj); + } else if (strcmp(pcode, "CSC") == 0) { + astCSCset(prj); + } else if (strcmp(pcode, "QSC") == 0) { + astQSCset(prj); + } else if (strcmp(pcode, "HPX") == 0) { + astHPXset(prj); + } else if (strcmp(pcode, "XPH") == 0) { + astXPHset(prj); + } else { + /* Unrecognized projection code. */ + return 1; + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astPRJfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + return prj->astPRJfwd(phi, theta, prj, x, y); +} + +/*--------------------------------------------------------------------------*/ + +int astPRJrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + return prj->astPRJrev(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->astPRJfwd Pointer to astAZPfwd(). +* prj->astPRJrev Pointer to astAZPrev(). +*===========================================================================*/ + +int astAZPset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "AZP"); + prj->flag = icopysign(WCS__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] = astCosd(prj->p[2]); + if (prj->w[3] == 0.0) { + return 1; + } + + prj->w[2] = 1.0/prj->w[3]; + prj->w[4] = astSind(prj->p[2]); + prj->w[1] = prj->w[4] / prj->w[3]; + + if (fabs(prj->p[1]) > 1.0) { + prj->w[5] = astASind(-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->astPRJfwd = astAZPfwd; + prj->astPRJrev = astAZPrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astAZPfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, b, cphi, cthe, r, s, t; + + if (abs(prj->flag) != WCS__AZP) { + if (astAZPset(prj)) return 1; + } + + cphi = astCosd(phi); + cthe = astCosd(theta); + + s = prj->w[1]*cphi; + t = (prj->p[1] + astSind(theta)) + cthe*s; + if (t == 0.0) { + return 2; + } + + r = prj->w[0]*cthe/t; + *x = r*astSind(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 = astATand(-s); + t = astASind(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 astAZPrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double a, b, r, s, t, ycosg; + const double tol = 1.0e-13; + + if (abs(prj->flag) != WCS__AZP) { + if (astAZPset(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 = astATan2d(x, -ycosg); + + s = r / (prj->w[0] + y*prj->w[4]); + t = s*prj->p[1]/sqrt(s*s + 1.0); + + s = astATan2d(1.0, s); + + if (fabs(t) > 1.0) { + t = copysign(90.0,t); + if (fabs(t) > 1.0+tol) { + return 2; + } + } else { + t = astASind(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->astPRJfwd Pointer to astSZPfwd(). +* prj->astPRJrev Pointer to astSZPrev(). +*===========================================================================*/ + +int astSZPset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "SZP"); + prj->flag = icopysign(WCS__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] * astSind(prj->p[3]) + 1.0; + if (prj->w[3] == 0.0) { + return 1; + } + + prj->w[1] = -prj->p[1] * astCosd(prj->p[3]) * astSind(prj->p[2]); + prj->w[2] = prj->p[1] * astCosd(prj->p[3]) * astCosd(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] = astASind(1.0 - prj->w[3]); + } else { + prj->w[8] = -90.0; + } + + prj->astPRJfwd = astSZPfwd; + prj->astPRJrev = astSZPrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astSZPfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, b, cphi, cthe, s, sphi, t; + + if (abs(prj->flag) != WCS__SZP) { + if (astSZPset(prj)) return 1; + } + + cphi = astCosd(phi); + sphi = astSind(phi); + cthe = astCosd(theta); + s = 1.0 - astSind(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 = astATan2d(s, prj->w[3] - 1.0); + t = astASind(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 astSZPrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *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) != WCS__SZP) { + if (astSZPset(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 = astASind(sthe); + + z = 1.0 - sthe; + } + + *phi = astATan2d(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->astPRJfwd Pointer to astTANfwd(). +* prj->astPRJrev Pointer to astTANrev(). +*===========================================================================*/ + +int astTANset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "TAN"); + prj->flag = icopysign(WCS__TAN, prj->flag); + prj->phi0 = 0.0; + prj->theta0 = 90.0; + + if (prj->r0 == 0.0) prj->r0 = R2D; + + prj->astPRJfwd = astTANfwd; + prj->astPRJrev = astTANrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astTANfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double r, s; + + if (abs(prj->flag) != WCS__TAN) { + if(astTANset(prj)) return 1; + } + + s = astSind(theta); + if (s == 0.0) { + return 2; + } + + r = prj->r0*astCosd(theta)/s; + *x = r*astSind(phi); + *y = -r*astCosd(phi); + + if (prj->flag > 0 && s < 0.0) { + return 2; + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astTANrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double r; + + if (abs(prj->flag) != WCS__TAN) { + if (astTANset(prj)) return 1; + } + + r = sqrt(x*x + y*y); + if (r == 0.0) { + *phi = 0.0; + } else { + *phi = astATan2d(x, -y); + } + *theta = astATan2d(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->astPRJfwd Pointer to astSTGfwd(). +* prj->astPRJrev Pointer to astSTGrev(). +*===========================================================================*/ + +int astSTGset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "STG"); + prj->flag = WCS__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->astPRJfwd = astSTGfwd; + prj->astPRJrev = astSTGrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astSTGfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double r, s; + + if (prj->flag != WCS__STG) { + if (astSTGset(prj)) return 1; + } + + s = 1.0 + astSind(theta); + if (s == 0.0) { + return 2; + } + + r = prj->w[0]*astCosd(theta)/s; + *x = r*astSind(phi); + *y = -r*astCosd(phi); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astSTGrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double r; + + if (prj->flag != WCS__STG) { + if (astSTGset(prj)) return 1; + } + + r = sqrt(x*x + y*y); + if (r == 0.0) { + *phi = 0.0; + } else { + *phi = astATan2d(x, -y); + } + *theta = 90.0 - 2.0*astATand(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->astPRJfwd Pointer to astSINfwd(). +* prj->astPRJrev Pointer to astSINrev(). +*===========================================================================*/ + +int astSINset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "SIN"); + prj->flag = icopysign(WCS__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->astPRJfwd = astSINfwd; + prj->astPRJrev = astSINrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astSINfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double cphi, cthe, sphi, t, z; + + if (abs(prj->flag) != WCS__SIN) { + if (astSINset(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 - astSind(theta); + cthe = astCosd(theta); + } + + cphi = astCosd(phi); + sphi = astSind(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 = -astATand(prj->p[1]*sphi - prj->p[2]*cphi); + if (theta < t) { + return 2; + } + } + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astSINrev (x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *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) != WCS__SIN) { + if (astSINset(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 = astATan2d(x0, -y0); + } else { + *phi = 0.0; + } + + if (r2 < 0.5) { + *theta = astACosd(sqrt(r2)); + } else if (r2 <= 1.0) { + *theta = astASind(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 = astASind(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 = astATan2d(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->astPRJfwd Pointer to astARCfwd(). +* prj->astPRJrev Pointer to astARCrev(). +*===========================================================================*/ + +int astARCset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "ARC"); + prj->flag = WCS__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->astPRJfwd = astARCfwd; + prj->astPRJrev = astARCrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astARCfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double r; + + if (prj->flag != WCS__ARC) { + if (astARCset(prj)) return 1; + } + + r = prj->w[0]*(90.0 - theta); + *x = r*astSind(phi); + *y = -r*astCosd(phi); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astARCrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double r; + + if (prj->flag != WCS__ARC) { + if (astARCset(prj)) return 1; + } + + r = sqrt(x*x + y*y); + if (r == 0.0) { + *phi = 0.0; + } else { + *phi = astATan2d(x, -y); + } + *theta = 90.0 - r*prj->w[1]; + + return 0; +} + +/*============================================================================ +* ZPN: zenithal/azimuthal polynomial projection. +* +* Given: +* prj->p[0:WCSLIB_MXPAR-1] 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->astPRJfwd Pointer to astZPNfwd(). +* prj->astPRJrev Pointer to astZPNrev(). +*===========================================================================*/ + +int astZPNset(prj) + +struct AstPrjPrm *prj; + +{ + int i, j, k, plen; + double d, d1, d2, r, zd, zd1, zd2; + const double tol = 1.0e-13; + + strcpy(prj->code, "ZPN"); + prj->flag = icopysign(WCS__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. */ + plen = astSizeOf( prj->p )/sizeof( double ); + for (k = plen-1; k >= 0 && prj->p[k] == 0.0; k--); + if (k < 0) return 1; + + prj->n = k; + + 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; + } + + prj->astPRJfwd = astZPNfwd; + prj->astPRJrev = astZPNrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astZPNfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + int j; + double r, s; + + if (abs(prj->flag) != WCS__ZPN) { + if (astZPNset(prj)) return 1; + } + + s = (90.0 - theta)*D2R; + + r = 0.0; + for (j = prj->n; j >= 0; j--) { + r = r*s + prj->p[j]; + } + r = prj->r0*r; + + *x = r*astSind(phi); + *y = -r*astCosd(phi); + + if (prj->flag > 0 && s > prj->w[0] && prj->n > 2 ) { + return 2; + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astZPNrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *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) != WCS__ZPN) { + if (astZPNset(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 = astATan2d(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->astPRJfwd Pointer to astZEAfwd(). +* prj->astPRJrev Pointer to astZEArev(). +*===========================================================================*/ + +int astZEAset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "ZEA"); + prj->flag = WCS__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->astPRJfwd = astZEAfwd; + prj->astPRJrev = astZEArev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astZEAfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double r; + + if (prj->flag != WCS__ZEA) { + if (astZEAset(prj)) return 1; + } + + r = prj->w[0]*astSind((90.0 - theta)/2.0); + *x = r*astSind(phi); + *y = -r*astCosd(phi); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astZEArev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double r, s; + const double tol = 1.0e-12; + + if (prj->flag != WCS__ZEA) { + if (astZEAset(prj)) return 1; + } + + r = sqrt(x*x + y*y); + if (r == 0.0) { + *phi = 0.0; + } else { + *phi = astATan2d(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*astASind(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->astPRJfwd Pointer to astAIRfwd(). +* prj->astPRJrev Pointer to astAIRrev(). +*===========================================================================*/ + +int astAIRset(prj) + +struct AstPrjPrm *prj; + +{ + const double tol = 1.0e-4; + double cxi; + + strcpy(prj->code, "AIR"); + prj->flag = WCS__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 = astCosd((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->astPRJfwd = astAIRfwd; + prj->astPRJrev = astAIRrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astAIRfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double cxi, r, txi, xi; + + if (prj->flag != WCS__AIR) { + if (astAIRset(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 = astCosd((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*astSind(phi); + *y = -r*astCosd(phi); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astAIRrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *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 != WCS__AIR) { + if (astAIRset(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 = astACosd(cxi); + } + + if (r == 0.0) { + *phi = 0.0; + } else { + *phi = astATan2d(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->astPRJfwd Pointer to astCYPfwd(). +* prj->astPRJrev Pointer to astCYPrev(). +*===========================================================================*/ + +int astCYPset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "CYP"); + prj->flag = WCS__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->astPRJfwd = astCYPfwd; + prj->astPRJrev = astCYPrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCYPfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double s; + + if (prj->flag != WCS__CYP) { + if (astCYPset(prj)) return 1; + } + + s = prj->p[1] + astCosd(theta); + if (s == 0.0) { + return 2; + } + + *x = prj->w[0]*phi; + *y = prj->w[2]*astSind(theta)/s; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCYPrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double eta; + double a; + const double tol = 1.0e-13; + + if (prj->flag != WCS__CYP) { + if (astCYPset(prj)) return 1; + } + + *phi = x*prj->w[1]; + eta = y*prj->w[3]; + + a = eta*prj->p[1]/sqrt(eta*eta+1.0); + if( fabs( a ) < 1.0 ) { + *theta = astATan2d(eta,1.0) + astASind( a ); + + } else if( fabs( a ) < 1.0 + tol ) { + if( a > 0.0 ){ + *theta = astATan2d(eta,1.0) + 90.0; + } else { + *theta = astATan2d(eta,1.0) - 90.0; + } + + } else { + return 2; + } + + 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->astPRJfwd Pointer to astCEAfwd(). +* prj->astPRJrev Pointer to astCEArev(). +*===========================================================================*/ + +int astCEAset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "CEA"); + prj->flag = WCS__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->astPRJfwd = astCEAfwd; + prj->astPRJrev = astCEArev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCEAfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + if (prj->flag != WCS__CEA) { + if (astCEAset(prj)) return 1; + } + + *x = prj->w[0]*phi; + *y = prj->w[2]*astSind(theta); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCEArev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double s; + const double tol = 1.0e-13; + + if (prj->flag != WCS__CEA) { + if (astCEAset(prj)) return 1; + } + + s = y*prj->w[3]; + if (fabs(s) > 1.0) { + if (fabs(s) > 1.0+tol) { + return 2; + } + s = copysign(1.0,s); + } + + *phi = x*prj->w[1]; + *theta = astASind(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->astPRJfwd Pointer to astCARfwd(). +* prj->astPRJrev Pointer to astCARrev(). +*===========================================================================*/ + +int astCARset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "CAR"); + prj->flag = WCS__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->astPRJfwd = astCARfwd; + prj->astPRJrev = astCARrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCARfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + if (prj->flag != WCS__CAR) { + if (astCARset(prj)) return 1; + } + + *x = prj->w[0]*phi; + *y = prj->w[0]*theta; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCARrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + if (prj->flag != WCS__CAR) { + if (astCARset(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->astPRJfwd Pointer to astMERfwd(). +* prj->astPRJrev Pointer to astMERrev(). +*===========================================================================*/ + +int astMERset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "MER"); + prj->flag = WCS__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->astPRJfwd = astMERfwd; + prj->astPRJrev = astMERrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astMERfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + if (prj->flag != WCS__MER) { + if (astMERset(prj)) return 1; + } + + if (theta <= -90.0 || theta >= 90.0) { + return 2; + } + + *x = prj->w[0]*phi; + *y = prj->r0*log(astTand((90.0+theta)/2.0)); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astMERrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + if (prj->flag != WCS__MER) { + if (astMERset(prj)) return 1; + } + + *phi = x*prj->w[1]; + *theta = 2.0*astATand(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->astPRJfwd Pointer to astSFLfwd(). +* prj->astPRJrev Pointer to astSFLrev(). +*===========================================================================*/ + +int astSFLset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "SFL"); + prj->flag = WCS__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->astPRJfwd = astSFLfwd; + prj->astPRJrev = astSFLrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astSFLfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + if (prj->flag != WCS__SFL) { + if (astSFLset(prj)) return 1; + } + + *x = prj->w[0]*phi*astCosd(theta); + *y = prj->w[0]*theta; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astSFLrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double w; + + if (prj->flag != WCS__SFL) { + if (astSFLset(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->astPRJfwd Pointer to astPARfwd(). +* prj->astPRJrev Pointer to astPARrev(). +*===========================================================================*/ + +int astPARset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "PAR"); + prj->flag = WCS__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->astPRJfwd = astPARfwd; + prj->astPRJrev = astPARrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astPARfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double s; + + if (prj->flag != WCS__PAR) { + if (astPARset(prj)) return 1; + } + + s = astSind(theta/3.0); + *x = prj->w[0]*phi*(1.0 - 4.0*s*s); + *y = prj->w[2]*s; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astPARrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double s, t; + + if (prj->flag != WCS__PAR) { + if (astPARset(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*astASind(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->astPRJfwd Pointer to astMOLfwd(). +* prj->astPRJrev Pointer to astMOLrev(). +*===========================================================================*/ + +int astMOLset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "MOL"); + prj->flag = WCS__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->astPRJfwd = astMOLfwd; + prj->astPRJrev = astMOLrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astMOLfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + int j; + double gamma, resid, u, v, v0, v1; + const double tol = 1.0e-13; + + if (prj->flag != WCS__MOL) { + if (astMOLset(prj)) return 1; + } + + if (fabs(theta) == 90.0) { + *x = 0.0; + *y = copysign(prj->w[0],theta); + } else if (theta == 0.0) { + *x = prj->w[1]*phi; + *y = 0.0; + } else { + u = PI*astSind(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 astMOLrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double s, y0, z; + const double tol = 1.0e-12; + + if (prj->flag != WCS__MOL) { + if (astMOLset(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 = copysign(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 = copysign(1.0,z); + } + + *theta = astASind(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->astPRJfwd Pointer to astAITfwd(). +* prj->astPRJrev Pointer to astAITrev(). +*===========================================================================*/ + +int astAITset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "AIT"); + prj->flag = WCS__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->astPRJfwd = astAITfwd; + prj->astPRJrev = astAITrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astAITfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double cthe, w; + + if (prj->flag != WCS__AIT) { + if (astAITset(prj)) return 1; + } + + cthe = astCosd(theta); + w = sqrt(prj->w[0]/(1.0 + cthe*astCosd(phi/2.0))); + *x = 2.0*w*cthe*astSind(phi/2.0); + *y = w*astSind(theta); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astAITrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double s, u, xp, yp, z; + const double tol = 1.0e-13; + + if (prj->flag != WCS__AIT) { + if (astAITset(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 = copysign(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*astATan2d(yp, xp); + } + *theta = astASind(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->astPRJfwd Pointer to astCOPfwd(). +* prj->astPRJrev Pointer to astCOPrev(). +*===========================================================================*/ + +int astCOPset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "COP"); + prj->flag = icopysign(WCS__COP, prj->flag); + prj->phi0 = 0.0; + prj->theta0 = prj->p[1]; + + if (prj->r0 == 0.0) prj->r0 = R2D; + + prj->w[0] = astSind(prj->p[1]); + if (prj->w[0] == 0.0) { + return 1; + } + + prj->w[1] = 1.0/prj->w[0]; + + prj->w[3] = prj->r0*astCosd(prj->p[2]); + if (prj->w[3] == 0.0) { + return 1; + } + + prj->w[4] = 1.0/prj->w[3]; + prj->w[5] = 1.0/astTand(prj->p[1]); + + prj->w[2] = prj->w[3]*prj->w[5]; + + prj->astPRJfwd = astCOPfwd; + prj->astPRJrev = astCOPrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCOPfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, r, s, t; + + if (abs(prj->flag) != WCS__COP) { + if (astCOPset(prj)) return 1; + } + + t = theta - prj->p[1]; + s = astCosd(t); + if (s == 0.0) { + return 2; + } + + a = prj->w[0]*phi; + r = prj->w[2] - prj->w[3]*astSind(t)/s; + + *x = r*astSind(a); + *y = prj->w[2] - r*astCosd(a); + + if (prj->flag > 0 && r*prj->w[0] < 0.0) { + return 2; + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCOPrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double a, dy, r; + + if (abs(prj->flag) != WCS__COP) { + if (astCOPset(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 = astATan2d(x/r, dy/r); + } + + *phi = a*prj->w[1]; + *theta = prj->p[1] + astATand(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*astSind(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->astPRJfwd Pointer to astCOEfwd(). +* prj->astPRJrev Pointer to astCOErev(). +*===========================================================================*/ + +int astCOEset(prj) + +struct AstPrjPrm *prj; + +{ + double theta1, theta2; + + strcpy(prj->code, "COE"); + prj->flag = WCS__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] = (astSind(theta1) + astSind(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 + astSind(theta1)*astSind(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]*astSind(prj->p[1])); + + prj->astPRJfwd = astCOEfwd; + prj->astPRJrev = astCOErev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCOEfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, r; + + if (prj->flag != WCS__COE) { + if (astCOEset(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]*astSind(theta)); + } + + *x = r*astSind(a); + *y = prj->w[2] - r*astCosd(a); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCOErev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double a, dy, r, w; + const double tol = 1.0e-12; + + if (prj->flag != WCS__COE) { + if (astCOEset(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 = astATan2d(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 = astASind(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->astPRJfwd Pointer to astCODfwd(). +* prj->astPRJrev Pointer to astCODrev(). +*===========================================================================*/ + +int astCODset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "COD"); + prj->flag = WCS__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*astSind(prj->p[1])*D2R; + } else { + prj->w[0] = prj->r0*astSind(prj->p[1])*astSind(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*astCosd(prj->p[2])*astCosd(prj->p[1])/prj->w[0]; + prj->w[3] = prj->w[2] + prj->p[1]; + + prj->astPRJfwd = astCODfwd; + prj->astPRJrev = astCODrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCODfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, r; + + if (prj->flag != WCS__COD) { + if (astCODset(prj)) return 1; + } + + a = prj->w[0]*phi; + r = prj->w[3] - theta; + + *x = r*astSind(a); + *y = prj->w[2] - r*astCosd(a); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCODrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double a, dy, r; + + if (prj->flag != WCS__COD) { + if (astCODset(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 = astATan2d(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->astPRJfwd Pointer to astCOOfwd(). +* prj->astPRJrev Pointer to astCOOrev(). +*===========================================================================*/ + +int astCOOset(prj) + +struct AstPrjPrm *prj; + +{ + double cos1, cos2, tan1, tan2, theta1, theta2; + + strcpy(prj->code, "COO"); + prj->flag = WCS__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 = astTand((90.0 - theta1)/2.0); + cos1 = astCosd(theta1); + + if (theta1 == theta2) { + prj->w[0] = astSind(theta1); + } else { + tan2 = astTand((90.0 - theta2)/2.0); + cos2 = astCosd(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(astTand((90.0 - prj->p[1])/2.0),prj->w[0]); + prj->w[4] = 1.0/prj->w[3]; + + prj->astPRJfwd = astCOOfwd; + prj->astPRJrev = astCOOrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCOOfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, r; + + if (prj->flag != WCS__COO) { + if (astCOOset(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(astTand((90.0 - theta)/2.0),prj->w[0]); + } + + *x = r*astSind(a); + *y = prj->w[2] - r*astCosd(a); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCOOrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double a, dy, r; + + if (prj->flag != WCS__COO) { + if (astCOOset(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 = astATan2d(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*astATand(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->astPRJfwd Pointer to astBONfwd(). +* prj->astPRJrev Pointer to astBONrev(). +*===========================================================================*/ + +int astBONset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "BON"); + prj->flag = WCS__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*astCosd(prj->p[1])/astSind(prj->p[1]) + prj->p[1]; + } else { + prj->w[1] = prj->r0*D2R; + prj->w[2] = prj->r0*(astCosd(prj->p[1])/astSind(prj->p[1]) + prj->p[1]*D2R); + } + + prj->astPRJfwd = astBONfwd; + prj->astPRJrev = astBONrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astBONfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, r; + + if (prj->p[1] == 0.0) { + /* Sanson-Flamsteed. */ + return astSFLfwd(phi, theta, prj, x, y); + } + + if (prj->flag != WCS__BON) { + if (astBONset(prj)) return 1; + } + + r = prj->w[2] - theta*prj->w[1]; + a = prj->r0*phi*astCosd(theta)/r; + + *x = r*astSind(a); + *y = prj->w[2] - r*astCosd(a); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astBONrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double a, cthe, dy, r; + + if (prj->p[1] == 0.0) { + /* Sanson-Flamsteed. */ + return astSFLrev(x, y, prj, phi, theta); + } + + if (prj->flag != WCS__BON) { + if (astBONset(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 = astATan2d(x/r, dy/r); + } + + *theta = (prj->w[2] - r)/prj->w[1]; + cthe = astCosd(*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->astPRJfwd Pointer to astPCOfwd(). +* prj->astPRJrev Pointer to astPCOrev(). +*===========================================================================*/ + +int astPCOset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "PCO"); + prj->flag = WCS__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->astPRJfwd = astPCOfwd; + prj->astPRJrev = astPCOrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astPCOfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double a, cthe, cotthe, sthe; + + if (prj->flag != WCS__PCO) { + if (astPCOset(prj)) return 1; + } + + cthe = astCosd(theta); + sthe = astSind(theta); + a = phi*sthe; + + if (sthe == 0.0) { + *x = prj->w[0]*phi; + *y = 0.0; + } else { + cotthe = cthe/sthe; + *x = prj->r0*cotthe*astSind(a); + *y = prj->r0*(cotthe*(1.0 - astCosd(a)) + theta*D2R); + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astPCOrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *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 != WCS__PCO) { + if (astPCOset(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 = copysign(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 = astTand(*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 = astATan2d(yp, xp)/astSind(*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->astPRJfwd Pointer to astTSCfwd(). +* prj->astPRJrev Pointer to astTSCrev(). +*===========================================================================*/ + +int astTSCset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "TSC"); + prj->flag = WCS__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->astPRJfwd = astTSCfwd; + prj->astPRJrev = astTSCrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astTSCfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + int face; + double cthe, l, m, n, rho, x0, xf, y0, yf; + const double tol = 1.0e-12; + + x0 = 0.0; + xf = 0.0; + y0 = 0.0; + yf = 0.0; + + if (prj->flag != WCS__TSC) { + if (astTSCset(prj)) return 1; + } + + cthe = astCosd(theta); + l = cthe*astCosd(phi); + m = cthe*astSind(phi); + n = astSind(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 = copysign(1.0,xf); + } + if (fabs(yf) > 1.0) { + if (fabs(yf) > 1.0+tol) { + return 2; + } + yf = copysign(1.0,yf); + } + + *x = prj->w[0]*(xf + x0); + *y = prj->w[0]*(yf + y0); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astTSCrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double l, m, n, xf, yf; + + if (prj->flag != WCS__TSC) { + if (astTSCset(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 = astATan2d(m, l); + } + *theta = astASind(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->astPRJfwd Pointer to astCSCfwd(). +* prj->astPRJrev Pointer to astCSCrev(). +*===========================================================================*/ + +int astCSCset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "CSC"); + prj->flag = WCS__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->astPRJfwd = astCSCfwd; + prj->astPRJrev = astCSCrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCSCfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + int face; + float cthe, eta, l, m, n, rho, xi; + const float tol = 1.0e-7; + + float a, a2, a2b2, a4, ab, b, b2, b4, ca2, cb2, x0, xf, y0, yf; + 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; + + eta = 0.0; + xi = 0.0; + x0 = 0.0; + y0 = 0.0; + + if (prj->flag != WCS__CSC) { + if (astCSCset(prj)) return 1; + } + + cthe = astCosd(theta); + l = cthe*astCosd(phi); + m = cthe*astSind(phi); + n = astSind(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 = copysign(1.0,xf); + } + if (fabs(yf) > 1.0) { + if (fabs(yf) > 1.0+tol) { + return 2; + } + yf = copysign(1.0,yf); + } + + *x = prj->w[0]*(x0 + xf); + *y = prj->w[0]*(y0 + yf); + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astCSCrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + int face; + float l, m, n; + + 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; + + l = 0.0; + m = 0.0; + n = 0.0; + + if (prj->flag != WCS__CSC) { + if (astCSCset(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 = astATan2d(m, l); + } + *theta = astASind(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->astPRJfwd Pointer to astQSCfwd(). +* prj->astPRJrev Pointer to astQSCrev(). +*===========================================================================*/ + +int astQSCset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "QSC"); + prj->flag = WCS__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->astPRJfwd = astQSCfwd; + prj->astPRJrev = astQSCrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astQSCfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + int face; + double cthe, eta, l, m, n, omega, p, rho, rhu, t, tau, x0, xf, xi, y0, yf; + const double tol = 1.0e-12; + + eta = 0.0; + x0 = 0.0; + xf = 0.0; + xi = 0.0; + y0 = 0.0; + yf = 0.0; + + if (prj->flag != WCS__QSC) { + if (astQSCset(prj)) return 1; + } + + if (fabs(theta) == 90.0) { + *x = 0.0; + *y = copysign(2.0*prj->w[0],theta); + return 0; + } + + cthe = astCosd(theta); + l = cthe*astCosd(phi); + m = cthe*astSind(phi); + n = astSind(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)*(astATand(omega) - astASind(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)*(astATand(omega) - astASind(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)*(astATand(omega) - astASind(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)*(astATand(omega) - astASind(omega/sqrt(tau+tau))); + } + + if (fabs(xf) > 1.0) { + if (fabs(xf) > 1.0+tol) { + return 2; + } + xf = copysign(1.0,xf); + } + if (fabs(yf) > 1.0) { + if (fabs(yf) > 1.0+tol) { + return 2; + } + yf = copysign(1.0,yf); + } + + *x = prj->w[0]*(xf + x0); + *y = prj->w[0]*(yf + y0); + + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astQSCrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + int direct, face; + double l, m, n, omega, rho, rhu, tau, xf, yf, w; + const double tol = 1.0e-12; + + l = 0.0; + m = 0.0; + n = 0.0; + + if (prj->flag != WCS__QSC) { + if (astQSCset(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 = astSind(w)/(astCosd(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 = astSind(w)/(astCosd(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 = astATan2d(m, l); + } + *theta = astASind(n); + + return 0; +} + +/*============================================================================ +* HPX: HEALPix projection. +* +* Given: +* prj->p[1] H - the number of facets in longitude. +* prj->p[2] K - the number of facets in latitude +* +* Given and/or returned: +* prj->r0 Reset to 180/pi if 0. +* prj->phi0 Reset to 0.0 +* prj->theta0 Reset to 0.0 +* +* Returned: +* prj->flag HPX +* prj->code "HPX" +* prj->n True if K is odd. +* prj->w[0] r0*(pi/180) +* prj->w[1] (180/pi)/r0 +* prj->w[2] (K-1)/K +* prj->w[3] 90*K/H +* prj->w[4] (K+1)/2 +* prj->w[5] 90*(K-1)/H +* prj->w[6] 180/H +* prj->w[7] H/360 +* prj->w[8] (90*K/H)*r0*(pi/180) +* prj->w[9] (180/H)*r0*(pi/180) +* prj->astPRJfwd Pointer to astHPXfwd(). +* prj->astPRJrev Pointer to astHPXrev(). + + +*===========================================================================*/ + +int astHPXset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "HPX"); + prj->flag = WCS__HPX; + prj->phi0 = 0.0; + prj->theta0 = 0.0; + + prj->n = ((int)prj->p[2])%2; + + 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] = R2D/prj->r0; + } + + prj->w[2] = (prj->p[2] - 1.0) / prj->p[2]; + prj->w[3] = 90.0 * prj->p[2] / prj->p[1]; + prj->w[4] = (prj->p[2] + 1.0) / 2.0; + prj->w[5] = 90.0 * (prj->p[2] - 1.0) / prj->p[1]; + prj->w[6] = 180.0 / prj->p[1]; + prj->w[7] = prj->p[1] / 360.0; + prj->w[8] = prj->w[3] * prj->w[0]; + prj->w[9] = prj->w[6] * prj->w[0]; + + prj->astPRJfwd = astHPXfwd; + prj->astPRJrev = astHPXrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astHPXfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double abssin, sigma, sinthe, phic; + int hodd; + + if( prj->flag != WCS__HPX ) { + if( astHPXset( prj ) ) return 1; + } + + sinthe = astSind( theta ); + abssin = fabs( sinthe ); + +/* Equatorial zone */ + if( abssin <= prj->w[2] ) { + *x = prj->w[0] * phi; + *y = prj->w[8] * sinthe; + +/* Polar zone */ + } else { + +/* DSB - The expression for phic is conditioned differently to the + WCSLIB code in order to improve accuracy of the floor function for + arguments very slightly below an integer value. */ + hodd = ((int)prj->p[1]) % 2; + if( !prj->n && theta <= 0.0 ) hodd = 1 - hodd; + if( hodd ) { + phic = -180.0 + (2.0*floor( prj->w[7] * phi + 1/2 ) + prj->p[1] ) * prj->w[6]; + } else { + phic = -180.0 + (2.0*floor( prj->w[7] * phi ) + prj->p[1] + 1 ) * prj->w[6]; + } + + sigma = sqrt( prj->p[2]*( 1.0 - abssin )); + + *x = prj->w[0] *( phic + ( phi - phic )*sigma ); + + *y = prj->w[9] * ( prj->w[4] - sigma ); + if( theta < 0 ) *y = -*y; + + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astHPXrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double absy, sigma, t, yr, xc; + int hodd; + + if (prj->flag != WCS__HPX) { + if (astHPXset(prj)) return 1; + } + + yr = prj->w[1]*y; + absy = fabs( yr ); + +/* Equatorial zone */ + if( absy <= prj->w[5] ) { + *phi = prj->w[1] * x; + t = yr/prj->w[3]; + if( t < -1.0 || t > 1.0 ) { + return 2; + } else { + *theta = astASind( t ); + } + +/* Polar zone */ + } else if( absy <= 90 ){ + + +/* DSB - The expression for xc is conditioned differently to the + WCSLIB code in order to improve accuracy of the floor function for + arguments very slightly below an integer value. */ + hodd = ((int)prj->p[1]) % 2; + if( !prj->n && yr <= 0.0 ) hodd = 1 - hodd; + if( hodd ) { + xc = -180.0 + (2.0*floor( prj->w[7] * x + 1/2 ) + prj->p[1] ) * prj->w[6]; + } else { + xc = -180.0 + (2.0*floor( prj->w[7] * x ) + prj->p[1] + 1 ) * prj->w[6]; + } + + sigma = prj->w[4] - absy / prj->w[6]; + + if( sigma == 0.0 ) { + return 2; + } else { + + t = ( x - xc )/sigma; + if( fabs( t ) <= prj->w[6] ) { + *phi = prj->w[1] *( xc + t ); + } else { + return 2; + } + } + + t = 1.0 - sigma*sigma/prj->p[2]; + if( t < -1.0 || t > 1.0 ) { + return 2; + } else { + *theta = astASind( t ); + if( y < 0 ) *theta = -*theta; + } + + } else { + return 2; + } + + return 0; +} + +/*============================================================================ +* XPH: HEALPix polar, aka "butterfly" projection. +* +* Given and/or returned: +* prj->r0 Reset to 180/pi if 0. +* prj->phi0 Reset to 0.0 if undefined. +* prj->theta0 Reset to 0.0 if undefined. +* +* Returned: +* prj->flag XPH +* prj->code "XPH" +* prj->w[0] r0*(pi/180)/sqrt(2) +* prj->w[1] (180/pi)/r0/sqrt(2) +* prj->w[2] 2/3 +* prj->w[3] tol (= 1e-4) +* prj->w[4] sqrt(2/3)*(180/pi) +* prj->w[5] 90 - tol*sqrt(2/3)*(180/pi) +* prj->w[6] sqrt(3/2)*(pi/180) +* prj->astPRJfwd Pointer to astXPHfwd(). +* prj->astPRJrev Pointer to astXPHrev(). +*===========================================================================*/ + +int astXPHset(prj) + +struct AstPrjPrm *prj; + +{ + strcpy(prj->code, "XPH"); + prj->flag = WCS__XPH; + + 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] = R2D/prj->r0; + } + + prj->w[0] /= sqrt(2.0); + prj->w[1] /= sqrt(2.0); + prj->w[2] = 2.0/3.0; + prj->w[3] = 1e-4; + prj->w[4] = sqrt(prj->w[2])*R2D; + prj->w[5] = 90.0 - prj->w[3]*prj->w[4]; + prj->w[6] = sqrt(1.5)*D2R; + + prj->astPRJfwd = astXPHfwd; + prj->astPRJrev = astXPHrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int astXPHfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct AstPrjPrm *prj; +double *x, *y; + +{ + double abssin, chi, eta, psi, sigma, sinthe, xi; + + if (prj->flag != WCS__XPH) { + if (astXPHset(prj)) return 1; + } + + /* Do phi dependence. */ + chi = phi; + if (180.0 <= fabs(chi)) { + chi = fmod(chi, 360.0); + if (chi < -180.0) { + chi += 360.0; + } else if (180.0 <= chi) { + chi -= 360.0; + } + } + + /* phi is also recomputed from chi to avoid rounding problems. */ + chi += 180.0; + psi = fmod(chi, 90.0); + + /* y is used to hold phi (rounded). */ + *x = psi; + *y = chi - 180.0; + + /* Do theta dependence. */ + sinthe = astSind(theta); + abssin = fabs(sinthe); + + if (abssin <= prj->w[2]) { + /* Equatorial regime. */ + xi = *x; + eta = 67.5 * sinthe; + + } else { + /* Polar regime. */ + if (theta < prj->w[5]) { + sigma = sqrt(3.0*(1.0 - abssin)); + } else { + sigma = (90.0 - theta)*prj->w[6]; + } + + xi = 45.0 + (*x - 45.0)*sigma; + eta = 45.0 * (2.0 - sigma); + if (theta < 0.0) eta = -eta; + } + + xi -= 45.0; + eta -= 90.0; + + /* Recall that y holds phi. */ + if (*y < -90.0) { + *x = prj->w[0]*(-xi + eta); + *y = prj->w[0]*(-xi - eta); + + } else if (*y < 0.0) { + *x = prj->w[0]*(+xi + eta); + *y = prj->w[0]*(-xi + eta); + + } else if (*y < 90.0) { + *x = prj->w[0]*( xi - eta); + *y = prj->w[0]*( xi + eta); + + } else { + *x = prj->w[0]*(-xi - eta); + *y = prj->w[0]*( xi - eta); + } + + return 0; + +} + +/*--------------------------------------------------------------------------*/ + +int astXPHrev(x, y, prj, phi, theta) + +const double x, y; +struct AstPrjPrm *prj; +double *phi, *theta; + +{ + double abseta, eta, eta1, sigma, xi, xi1, xr, yr; + const double tol = 1.0e-12; + + if (prj->flag != WCS__XPH) { + if (astXPHset(prj)) return 1; + } + + + xr = x*prj->w[1]; + yr = y*prj->w[1]; + if (xr <= 0.0 && 0.0 < yr) { + xi1 = -xr - yr; + eta1 = xr - yr; + *phi = -180.0; + } else if (xr < 0.0 && yr <= 0.0) { + xi1 = xr - yr; + eta1 = xr + yr; + *phi = -90.0; + } else if (0.0 <= xr && yr < 0.0) { + xi1 = xr + yr; + eta1 = -xr + yr; + *phi = 0.0; + } else { + xi1 = -xr + yr; + eta1 = -xr - yr; + *phi = 90.0; + } + + xi = xi1 + 45.0; + eta = eta1 + 90.0; + abseta = fabs(eta); + + if (abseta <= 90.0) { + if (abseta <= 45.0) { + /* Equatorial regime. */ + *phi += xi; + *theta = astASind(eta/67.5); + + /* Bounds checking. */ + if (45.0+tol < fabs(xi1)) return 2; + + } else { + /* Polar regime. */ + sigma = (90.0 - abseta) / 45.0; + + /* Ensure an exact result for points on the boundary. */ + if (xr == 0.0) { + if (yr <= 0.0) { + *phi = 0.0; + } else { + *phi = 180.0; + } + } else if (yr == 0.0) { + if (xr < 0.0) { + *phi = -90.0; + } else { + *phi = 90.0; + } + } else { + *phi += 45.0 + xi1/sigma; + } + + if (sigma < prj->w[3]) { + *theta = 90.0 - sigma*prj->w[4]; + } else { + *theta = astASind(1.0 - sigma*sigma/3.0); + } + if (eta < 0.0) *theta = -(*theta); + + /* Bounds checking. */ + if (eta < -45.0 && eta+90.0+tol < fabs(xi1)) return 2; + } + + } else { + /* Beyond latitude range. */ + *phi = 0.0; + *theta = 0.0; + return 2; + } + + return 0; +} + + |