/*
|
/*
|
* cel.c
|
* cel.c
|
*
|
*
|
* Lower level spherical coordinate transformation and projection routines.
|
* Lower level spherical coordinate transformation and projection routines.
|
*
|
*
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
*
|
*
|
* This file part of: AstrOmatic WCS library
|
* This file part of: AstrOmatic WCS library
|
*
|
*
|
* Copyright: (C) 2000-2010 IAP/CNRS/UPMC
|
* Copyright: (C) 2000-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
|
* (C) 1995-1999 Mark Calabretta
|
* (C) 1995-1999 Mark Calabretta (original version)
|
*
|
|
* Authors: Emmanuel Bertin (this version)
|
|
* Mark Calabretta (original version)
|
|
*
|
*
|
* Licenses: GNU General Public License
|
* Licenses: GNU General Public License
|
*
|
*
|
* AstrOmatic software is free software: you can redistribute it and/or
|
* AstrOmatic software is free software: you can redistribute it and/or
|
* modify it under the terms of the GNU General Public License as
|
* modify it under the terms of the GNU General Public License as
|
* published by the Free Software Foundation, either version 3 of the
|
* published by the Free Software Foundation, either version 3 of the
|
* License, or (at your option) any later version.
|
* License, or (at your option) any later version.
|
* AstrOmatic software is distributed in the hope that it will be useful,
|
* AstrOmatic software is distributed in the hope that it will be useful,
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
* GNU General Public License for more details.
|
* GNU General Public License for more details.
|
* You should have received a copy of the GNU General Public License
|
* You should have received a copy of the GNU General Public License
|
* along with AstrOmatic software.
|
* along with AstrOmatic software.
|
* If not, see <http://www.gnu.org/licenses/>.
|
* If not, see <http://www.gnu.org/licenses/>.
|
*
|
*
|
* Last modified: 10/10/2010
|
* Last modified: 10/10/2010
|
*
|
*
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
|
/*=============================================================================
|
/*=============================================================================
|
*
|
*
|
* WCSLIB - an implementation of the FITS WCS proposal.
|
* WCSLIB - an implementation of the FITS WCS proposal.
|
* Copyright (C) 1995-1999, Mark Calabretta
|
* Copyright (C) 1995-1999, Mark Calabretta
|
*
|
*
|
* This library is free software; you can redistribute it and/or modify it
|
* 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
|
* 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
|
* by the Free Software Foundation; either version 2 of the License, or (at
|
* your option) any later version.
|
* your option) any later version.
|
*
|
*
|
* This library is distributed in the hope that it will be useful, but
|
* This library is distributed in the hope that it will be useful, but
|
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library
|
* General Public License for more details.
|
* General Public License for more details.
|
*
|
*
|
* You should have received a copy of the GNU Library General Public License
|
* 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,
|
* along with this library; if not, write to the Free Software Foundation,
|
* Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
* Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
*
|
*
|
* Correspondence concerning WCSLIB may be directed to:
|
* Correspondence concerning WCSLIB may be directed to:
|
* Internet email: mcalabre@atnf.csiro.au
|
* Internet email: mcalabre@atnf.csiro.au
|
* Postal address: Dr. Mark Calabretta,
|
* Postal address: Dr. Mark Calabretta,
|
* Australia Telescope National Facility,
|
* Australia Telescope National Facility,
|
* P.O. Box 76,
|
* P.O. Box 76,
|
* Epping, NSW, 2121,
|
* Epping, NSW, 2121,
|
* AUSTRALIA
|
* AUSTRALIA
|
*
|
*
|
*=============================================================================
|
*=============================================================================
|
*
|
*
|
* C routines which implement the FITS World Coordinate System (WCS)
|
* C routines which implement the FITS World Coordinate System (WCS)
|
* convention.
|
* convention.
|
*
|
*
|
* Summary of routines
|
* Summary of routines
|
* -------------------
|
* -------------------
|
* These routines are provided as drivers for the lower level spherical
|
* These routines are provided as drivers for the lower level spherical
|
* coordinate transformation and projection routines. There are separate
|
* coordinate transformation and projection routines. There are separate
|
* driver routines for the forward, celfwd(), and reverse, celrev(),
|
* driver routines for the forward, celfwd(), and reverse, celrev(),
|
* transformations.
|
* transformations.
|
*
|
*
|
* An initialization routine, celset(), computes intermediate values from
|
* An initialization routine, celset(), computes intermediate values from
|
* the transformation parameters but need not be called explicitly - see the
|
* the transformation parameters but need not be called explicitly - see the
|
* explanation of cel.flag below.
|
* explanation of cel.flag below.
|
*
|
*
|
*
|
*
|
* Initialization routine; celset()
|
* Initialization routine; celset()
|
* --------------------------------
|
* --------------------------------
|
* Initializes members of a celprm data structure which hold intermediate
|
* Initializes members of a celprm data structure which hold intermediate
|
* values. Note that this routine need not be called directly; it will be
|
* values. Note that this routine need not be called directly; it will be
|
* invoked by celfwd() and celrev() if the "flag" structure member is
|
* invoked by celfwd() and celrev() if the "flag" structure member is
|
* anything other than a predefined magic value.
|
* anything other than a predefined magic value.
|
*
|
*
|
* Given:
|
* Given:
|
* pcode[4] const char
|
* pcode[4] const char
|
* WCS projection code (see below).
|
* WCS projection code (see below).
|
*
|
*
|
* Given and returned:
|
* Given and returned:
|
* cel celprm* Spherical coordinate transformation parameters
|
* cel celprm* Spherical coordinate transformation parameters
|
* (see below).
|
* (see below).
|
* prj prjprm* Projection parameters (usage is described in the
|
* prj prjprm* Projection parameters (usage is described in the
|
* prologue to "proj.c").
|
* prologue to "proj.c").
|
*
|
*
|
* Function return value:
|
* Function return value:
|
* int Error status
|
* int Error status
|
* 0: Success.
|
* 0: Success.
|
* 1: Invalid coordinate transformation parameters.
|
* 1: Invalid coordinate transformation parameters.
|
* 2: Ill-conditioned coordinate transformation
|
* 2: Ill-conditioned coordinate transformation
|
* parameters.
|
* parameters.
|
*
|
*
|
* Forward transformation; celfwd()
|
* Forward transformation; celfwd()
|
* --------------------------------
|
* --------------------------------
|
* Compute (x,y) coordinates in the plane of projection from celestial
|
* Compute (x,y) coordinates in the plane of projection from celestial
|
* coordinates (lng,lat).
|
* coordinates (lng,lat).
|
*
|
*
|
* Given:
|
* Given:
|
* pcode[4] const char
|
* pcode[4] const char
|
* WCS projection code (see below).
|
* WCS projection code (see below).
|
* lng,lat const double
|
* lng,lat const double
|
* Celestial longitude and latitude of the projected
|
* Celestial longitude and latitude of the projected
|
* point, in degrees.
|
* point, in degrees.
|
*
|
*
|
* Given and returned:
|
* Given and returned:
|
* cel celprm* Spherical coordinate transformation parameters
|
* cel celprm* Spherical coordinate transformation parameters
|
* (see below).
|
* (see below).
|
*
|
*
|
* Returned:
|
* Returned:
|
* phi, double* Longitude and latitude in the native coordinate
|
* phi, double* Longitude and latitude in the native coordinate
|
* theta system of the projection, in degrees.
|
* theta system of the projection, in degrees.
|
*
|
*
|
* Given and returned:
|
* Given and returned:
|
* prj prjprm* Projection parameters (usage is described in the
|
* prj prjprm* Projection parameters (usage is described in the
|
* prologue to "proj.c").
|
* prologue to "proj.c").
|
*
|
*
|
* Returned:
|
* Returned:
|
* x,y double* Projected coordinates, "degrees".
|
* x,y double* Projected coordinates, "degrees".
|
*
|
*
|
* Function return value:
|
* Function return value:
|
* int Error status
|
* int Error status
|
* 0: Success.
|
* 0: Success.
|
* 1: Invalid coordinate transformation parameters.
|
* 1: Invalid coordinate transformation parameters.
|
* 2: Invalid projection parameters.
|
* 2: Invalid projection parameters.
|
* 3: Invalid value of (lng,lat).
|
* 3: Invalid value of (lng,lat).
|
*
|
*
|
* Reverse transformation; celrev()
|
* Reverse transformation; celrev()
|
* --------------------------------
|
* --------------------------------
|
* Compute the celestial coordinates (lng,lat) of the point with projected
|
* Compute the celestial coordinates (lng,lat) of the point with projected
|
* coordinates (x,y).
|
* coordinates (x,y).
|
*
|
*
|
* Given:
|
* Given:
|
* pcode[4] const char
|
* pcode[4] const char
|
* WCS projection code (see below).
|
* WCS projection code (see below).
|
* x,y const double
|
* x,y const double
|
* Projected coordinates, "degrees".
|
* Projected coordinates, "degrees".
|
*
|
*
|
* Given and returned:
|
* Given and returned:
|
* prj prjprm* Projection parameters (usage is described in the
|
* prj prjprm* Projection parameters (usage is described in the
|
* prologue to "proj.c").
|
* prologue to "proj.c").
|
*
|
*
|
* Returned:
|
* Returned:
|
* phi, double* Longitude and latitude in the native coordinate
|
* phi, double* Longitude and latitude in the native coordinate
|
* theta system of the projection, in degrees.
|
* theta system of the projection, in degrees.
|
*
|
*
|
* Given and returned:
|
* Given and returned:
|
* cel celprm* Spherical coordinate transformation parameters
|
* cel celprm* Spherical coordinate transformation parameters
|
* (see below).
|
* (see below).
|
*
|
*
|
* Returned:
|
* Returned:
|
* lng,lat double* Celestial longitude and latitude of the projected
|
* lng,lat double* Celestial longitude and latitude of the projected
|
* point, in degrees.
|
* point, in degrees.
|
*
|
*
|
* Function return value:
|
* Function return value:
|
* int Error status
|
* int Error status
|
* 0: Success.
|
* 0: Success.
|
* 1: Invalid coordinate transformation parameters.
|
* 1: Invalid coordinate transformation parameters.
|
* 2: Invalid projection parameters.
|
* 2: Invalid projection parameters.
|
* 3: Invalid value of (x,y).
|
* 3: Invalid value of (x,y).
|
*
|
*
|
* Coordinate transformation parameters
|
* Coordinate transformation parameters
|
* ------------------------------------
|
* ------------------------------------
|
* The celprm struct consists of the following:
|
* The celprm struct consists of the following:
|
*
|
*
|
* int flag
|
* int flag
|
* The celprm struct contains pointers to the forward and reverse
|
* The celprm struct contains pointers to the forward and reverse
|
* projection routines as well as intermediaries computed from the
|
* projection routines as well as intermediaries computed from the
|
* reference coordinates (see below). Whenever the projection code
|
* reference coordinates (see below). Whenever the projection code
|
* (pcode) or any of ref[4] are set or changed then this flag must be
|
* (pcode) or any of ref[4] are set or changed then this flag must be
|
* set to zero to signal the initialization routine, celset(), to
|
* set to zero to signal the initialization routine, celset(), to
|
* redetermine the function pointers and recompute intermediaries.
|
* redetermine the function pointers and recompute intermediaries.
|
* Once this has been done pcode itself is ignored.
|
* Once this has been done pcode itself is ignored.
|
*
|
*
|
* double ref[4]
|
* double ref[4]
|
* The first pair of values should be set to the celestial longitude
|
* The first pair of values should be set to the celestial longitude
|
* and latitude (usually right ascension and declination) of the
|
* and latitude (usually right ascension and declination) of the
|
* reference point of the projection.
|
* reference point of the projection.
|
*
|
*
|
* The second pair of values are the native longitude and latitude of
|
* The second pair of values are the native longitude and latitude of
|
* the pole of the celestial coordinate system and correspond to the
|
* the pole of the celestial coordinate system and correspond to the
|
* FITS keywords LONGPOLE and LATPOLE.
|
* FITS keywords LONGPOLE and LATPOLE.
|
*
|
*
|
* LONGPOLE defaults to 0 degrees if the celestial latitude of the
|
* LONGPOLE defaults to 0 degrees if the celestial latitude of the
|
* reference point of the projection is greater than the native
|
* reference point of the projection is greater than the native
|
* latitude, otherwise 180 degrees. (This is the condition for the
|
* latitude, otherwise 180 degrees. (This is the condition for the
|
* celestial latitude to increase in the same direction as the native
|
* celestial latitude to increase in the same direction as the native
|
* latitude at the reference point.) ref[2] may be set to 999.0 to
|
* latitude at the reference point.) ref[2] may be set to 999.0 to
|
* indicate that the correct default should be substituted.
|
* indicate that the correct default should be substituted.
|
*
|
*
|
* In some circumstances the latitude of the native pole may be
|
* In some circumstances the latitude of the native pole may be
|
* determined by the first three values only to within a sign and
|
* determined by the first three values only to within a sign and
|
* LATPOLE is used to choose between the two solutions. LATPOLE is
|
* LATPOLE is used to choose between the two solutions. LATPOLE is
|
* set in ref[3] and the solution closest to this value is used to
|
* set in ref[3] and the solution closest to this value is used to
|
* reset ref[3]. It is therefore legitimate, for example, to set
|
* reset ref[3]. It is therefore legitimate, for example, to set
|
* ref[3] to 999.0 to choose the more northerly solution - the default
|
* ref[3] to 999.0 to choose the more northerly solution - the default
|
* if the LATPOLE card is omitted from the FITS header. For the
|
* if the LATPOLE card is omitted from the FITS header. For the
|
* special case where the reference point of the projection is at
|
* special case where the reference point of the projection is at
|
* native latitude zero, its celestial latitude is zero, and
|
* native latitude zero, its celestial latitude is zero, and
|
* LONGPOLE = +/- 90 then the native latitude of the pole is not
|
* LONGPOLE = +/- 90 then the native latitude of the pole is not
|
* determined by the first three reference values and LATPOLE
|
* determined by the first three reference values and LATPOLE
|
* specifies it completely.
|
* specifies it completely.
|
*
|
*
|
* The remaining members of the celprm struct are maintained by the
|
* The remaining members of the celprm struct are maintained by the
|
* initialization routines and should not be modified. This is done for 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
|
* sake of efficiency and to allow an arbitrary number of contexts to be
|
* maintained simultaneously.
|
* maintained simultaneously.
|
*
|
*
|
* double euler[5]
|
* double euler[5]
|
* Euler angles and associated intermediaries derived from the
|
* Euler angles and associated intermediaries derived from the
|
* coordinate reference values.
|
* coordinate reference values.
|
* int (*prjfwd)()
|
* int (*prjfwd)()
|
* int (*prjrev)()
|
* int (*prjrev)()
|
* Pointers to the forward and reverse projection routines.
|
* Pointers to the forward and reverse projection routines.
|
*
|
*
|
*
|
*
|
* WCS projection codes
|
* WCS projection codes
|
* --------------------
|
* --------------------
|
* Zenithals/azimuthals:
|
* Zenithals/azimuthals:
|
* AZP: zenithal/azimuthal perspective
|
* AZP: zenithal/azimuthal perspective
|
* TAN: gnomonic
|
* TAN: gnomonic
|
* SIN: synthesis (generalized orthographic)
|
* SIN: synthesis (generalized orthographic)
|
* STG: stereographic
|
* STG: stereographic
|
* ARC: zenithal/azimuthal equidistant
|
* ARC: zenithal/azimuthal equidistant
|
* ZPN: zenithal/azimuthal polynomial
|
* ZPN: zenithal/azimuthal polynomial
|
* ZEA: zenithal/azimuthal equal area
|
* ZEA: zenithal/azimuthal equal area
|
* AIR: Airy
|
* AIR: Airy
|
* TNX: IRAF's polynomial correction to TAN
|
* TNX: IRAF's polynomial correction to TAN
|
*
|
*
|
* Cylindricals:
|
* Cylindricals:
|
* CYP: cylindrical perspective
|
* CYP: cylindrical perspective
|
* CAR: Cartesian
|
* CAR: Cartesian
|
* MER: Mercator
|
* MER: Mercator
|
* CEA: cylindrical equal area
|
* CEA: cylindrical equal area
|
*
|
*
|
* Conics:
|
* Conics:
|
* COP: conic perspective
|
* COP: conic perspective
|
* COD: conic equidistant
|
* COD: conic equidistant
|
* COE: conic equal area
|
* COE: conic equal area
|
* COO: conic orthomorphic
|
* COO: conic orthomorphic
|
*
|
*
|
* Polyconics:
|
* Polyconics:
|
* BON: Bonne
|
* BON: Bonne
|
* PCO: polyconic
|
* PCO: polyconic
|
*
|
*
|
* Pseudo-cylindricals:
|
* Pseudo-cylindricals:
|
* GLS: Sanson-Flamsteed (global sinusoidal)
|
* GLS: Sanson-Flamsteed (global sinusoidal)
|
* PAR: parabolic
|
* PAR: parabolic
|
* MOL: Mollweide
|
* MOL: Mollweide
|
*
|
*
|
* Conventional:
|
* Conventional:
|
* AIT: Hammer-Aitoff
|
* AIT: Hammer-Aitoff
|
*
|
*
|
* Quad-cubes:
|
* Quad-cubes:
|
* CSC: COBE quadrilateralized spherical cube
|
* CSC: COBE quadrilateralized spherical cube
|
* QSC: quadrilateralized spherical cube
|
* QSC: quadrilateralized spherical cube
|
* TSC: tangential spherical cube
|
* TSC: tangential spherical cube
|
*
|
*
|
* Author: Mark Calabretta, Australia Telescope National Facility
|
* Author: Mark Calabretta, Australia Telescope National Facility
|
* IRAF's TNX added by E.Bertin 2000/03/28
|
* IRAF's TNX added by E.Bertin 2000/03/28
|
* Filtering of abs(phi)>180 and abs(theta)>90 added by E.Bertin 2000/11/11
|
* Filtering of abs(phi)>180 and abs(theta)>90 added by E.Bertin 2000/11/11
|
* $Id: cel.c,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
|
* $Id: cel.c,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
|
*===========================================================================*/
|
*===========================================================================*/
|
|
|
#ifdef HAVE_CONFIG_H
|
#ifdef HAVE_CONFIG_H
|
#include "config.h"
|
#include "config.h"
|
#endif
|
#endif
|
|
|
#ifdef HAVE_MATHIMF_H
|
#ifdef HAVE_MATHIMF_H
|
#include <mathimf.h>
|
#include <mathimf.h>
|
#else
|
#else
|
#include <math.h>
|
#include <math.h>
|
#endif
|
#endif
|
#include <string.h>
|
#include <string.h>
|
#include "wcstrig.h"
|
#include "wcstrig.h"
|
#include "cel.h"
|
#include "cel.h"
|
#include "sph.h"
|
#include "sph.h"
|
#include "tnx.h"
|
#include "tnx.h"
|
|
|
int npcode = 26;
|
int npcode = 26;
|
char pcodes[26][4] =
|
char pcodes[26][4] =
|
{"AZP", "TAN", "SIN", "STG", "ARC", "ZPN", "ZEA", "AIR", "CYP", "CAR",
|
{"AZP", "TAN", "SIN", "STG", "ARC", "ZPN", "ZEA", "AIR", "CYP", "CAR",
|
"MER", "CEA", "COP", "COD", "COE", "COO", "BON", "PCO", "GLS", "PAR",
|
"MER", "CEA", "COP", "COD", "COE", "COO", "BON", "PCO", "GLS", "PAR",
|
"AIT", "MOL", "CSC", "QSC", "TSC", "TNX"};
|
"AIT", "MOL", "CSC", "QSC", "TSC", "TNX"};
|
|
|
/* Map error number to error message for each function. */
|
/* Map error number to error message for each function. */
|
const char *celset_errmsg[] = {
|
const char *celset_errmsg[] = {
|
0,
|
0,
|
"Invalid coordinate transformation parameters",
|
"Invalid coordinate transformation parameters",
|
"Ill-conditioned coordinate transformation parameters"};
|
"Ill-conditioned coordinate transformation parameters"};
|
|
|
const char *celfwd_errmsg[] = {
|
const char *celfwd_errmsg[] = {
|
0,
|
0,
|
"Invalid coordinate transformation parameters",
|
"Invalid coordinate transformation parameters",
|
"Invalid projection parameters",
|
"Invalid projection parameters",
|
"Invalid value of (lng,lat)"};
|
"Invalid value of (lng,lat)"};
|
|
|
const char *celrev_errmsg[] = {
|
const char *celrev_errmsg[] = {
|
0,
|
0,
|
"Invalid coordinate transformation parameters",
|
"Invalid coordinate transformation parameters",
|
"Invalid projection parameters",
|
"Invalid projection parameters",
|
"Invalid value of (x,y)"};
|
"Invalid value of (x,y)"};
|
|
|
|
|
int celset(pcode, cel, prj)
|
int celset(pcode, cel, prj)
|
|
|
const char pcode[4];
|
const char pcode[4];
|
struct celprm *cel;
|
struct celprm *cel;
|
struct prjprm *prj;
|
struct prjprm *prj;
|
|
|
{
|
{
|
int dophip;
|
int dophip;
|
const double tol = 1.0e-10;
|
const double tol = 1.0e-10;
|
double clat0, cphip, cthe0, theta0, slat0, sphip, sthe0;
|
double clat0, cphip, cthe0, theta0, slat0, sphip, sthe0;
|
double latp, latp1, latp2;
|
double latp, latp1, latp2;
|
double u, v, x, y, z;
|
double u, v, x, y, z;
|
|
|
/* Set pointers to the forward and reverse projection routines. */
|
/* Set pointers to the forward and reverse projection routines. */
|
if (strcmp(pcode, "AZP") == 0) {
|
if (strcmp(pcode, "AZP") == 0) {
|
cel->prjfwd = azpfwd;
|
cel->prjfwd = azpfwd;
|
cel->prjrev = azprev;
|
cel->prjrev = azprev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "TAN") == 0) {
|
} else if (strcmp(pcode, "TAN") == 0) {
|
cel->prjfwd = tanfwd;
|
cel->prjfwd = tanfwd;
|
cel->prjrev = tanrev;
|
cel->prjrev = tanrev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "SIN") == 0) {
|
} else if (strcmp(pcode, "SIN") == 0) {
|
cel->prjfwd = sinfwd;
|
cel->prjfwd = sinfwd;
|
cel->prjrev = sinrev;
|
cel->prjrev = sinrev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "STG") == 0) {
|
} else if (strcmp(pcode, "STG") == 0) {
|
cel->prjfwd = stgfwd;
|
cel->prjfwd = stgfwd;
|
cel->prjrev = stgrev;
|
cel->prjrev = stgrev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "ARC") == 0) {
|
} else if (strcmp(pcode, "ARC") == 0) {
|
cel->prjfwd = arcfwd;
|
cel->prjfwd = arcfwd;
|
cel->prjrev = arcrev;
|
cel->prjrev = arcrev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "ZPN") == 0) {
|
} else if (strcmp(pcode, "ZPN") == 0) {
|
cel->prjfwd = zpnfwd;
|
cel->prjfwd = zpnfwd;
|
cel->prjrev = zpnrev;
|
cel->prjrev = zpnrev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "ZEA") == 0) {
|
} else if (strcmp(pcode, "ZEA") == 0) {
|
cel->prjfwd = zeafwd;
|
cel->prjfwd = zeafwd;
|
cel->prjrev = zearev;
|
cel->prjrev = zearev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "AIR") == 0) {
|
} else if (strcmp(pcode, "AIR") == 0) {
|
cel->prjfwd = airfwd;
|
cel->prjfwd = airfwd;
|
cel->prjrev = airrev;
|
cel->prjrev = airrev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else if (strcmp(pcode, "CYP") == 0) {
|
} else if (strcmp(pcode, "CYP") == 0) {
|
cel->prjfwd = cypfwd;
|
cel->prjfwd = cypfwd;
|
cel->prjrev = cyprev;
|
cel->prjrev = cyprev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "CAR") == 0) {
|
} else if (strcmp(pcode, "CAR") == 0) {
|
cel->prjfwd = carfwd;
|
cel->prjfwd = carfwd;
|
cel->prjrev = carrev;
|
cel->prjrev = carrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "MER") == 0) {
|
} else if (strcmp(pcode, "MER") == 0) {
|
cel->prjfwd = merfwd;
|
cel->prjfwd = merfwd;
|
cel->prjrev = merrev;
|
cel->prjrev = merrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "CEA") == 0) {
|
} else if (strcmp(pcode, "CEA") == 0) {
|
cel->prjfwd = ceafwd;
|
cel->prjfwd = ceafwd;
|
cel->prjrev = cearev;
|
cel->prjrev = cearev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "COP") == 0) {
|
} else if (strcmp(pcode, "COP") == 0) {
|
cel->prjfwd = copfwd;
|
cel->prjfwd = copfwd;
|
cel->prjrev = coprev;
|
cel->prjrev = coprev;
|
theta0 = prj->p[1];
|
theta0 = prj->p[1];
|
} else if (strcmp(pcode, "COD") == 0) {
|
} else if (strcmp(pcode, "COD") == 0) {
|
cel->prjfwd = codfwd;
|
cel->prjfwd = codfwd;
|
cel->prjrev = codrev;
|
cel->prjrev = codrev;
|
theta0 = prj->p[1];
|
theta0 = prj->p[1];
|
} else if (strcmp(pcode, "COE") == 0) {
|
} else if (strcmp(pcode, "COE") == 0) {
|
cel->prjfwd = coefwd;
|
cel->prjfwd = coefwd;
|
cel->prjrev = coerev;
|
cel->prjrev = coerev;
|
theta0 = prj->p[1];
|
theta0 = prj->p[1];
|
} else if (strcmp(pcode, "COO") == 0) {
|
} else if (strcmp(pcode, "COO") == 0) {
|
cel->prjfwd = coofwd;
|
cel->prjfwd = coofwd;
|
cel->prjrev = coorev;
|
cel->prjrev = coorev;
|
theta0 = prj->p[1];
|
theta0 = prj->p[1];
|
} else if (strcmp(pcode, "BON") == 0) {
|
} else if (strcmp(pcode, "BON") == 0) {
|
cel->prjfwd = bonfwd;
|
cel->prjfwd = bonfwd;
|
cel->prjrev = bonrev;
|
cel->prjrev = bonrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "PCO") == 0) {
|
} else if (strcmp(pcode, "PCO") == 0) {
|
cel->prjfwd = pcofwd;
|
cel->prjfwd = pcofwd;
|
cel->prjrev = pcorev;
|
cel->prjrev = pcorev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "GLS") == 0) {
|
} else if (strcmp(pcode, "GLS") == 0) {
|
cel->prjfwd = glsfwd;
|
cel->prjfwd = glsfwd;
|
cel->prjrev = glsrev;
|
cel->prjrev = glsrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "PAR") == 0) {
|
} else if (strcmp(pcode, "PAR") == 0) {
|
cel->prjfwd = parfwd;
|
cel->prjfwd = parfwd;
|
cel->prjrev = parrev;
|
cel->prjrev = parrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "AIT") == 0) {
|
} else if (strcmp(pcode, "AIT") == 0) {
|
cel->prjfwd = aitfwd;
|
cel->prjfwd = aitfwd;
|
cel->prjrev = aitrev;
|
cel->prjrev = aitrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "MOL") == 0) {
|
} else if (strcmp(pcode, "MOL") == 0) {
|
cel->prjfwd = molfwd;
|
cel->prjfwd = molfwd;
|
cel->prjrev = molrev;
|
cel->prjrev = molrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "CSC") == 0) {
|
} else if (strcmp(pcode, "CSC") == 0) {
|
cel->prjfwd = cscfwd;
|
cel->prjfwd = cscfwd;
|
cel->prjrev = cscrev;
|
cel->prjrev = cscrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "QSC") == 0) {
|
} else if (strcmp(pcode, "QSC") == 0) {
|
cel->prjfwd = qscfwd;
|
cel->prjfwd = qscfwd;
|
cel->prjrev = qscrev;
|
cel->prjrev = qscrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "TSC") == 0) {
|
} else if (strcmp(pcode, "TSC") == 0) {
|
cel->prjfwd = tscfwd;
|
cel->prjfwd = tscfwd;
|
cel->prjrev = tscrev;
|
cel->prjrev = tscrev;
|
theta0 = 0.0;
|
theta0 = 0.0;
|
} else if (strcmp(pcode, "TNX") == 0) {
|
} else if (strcmp(pcode, "TNX") == 0) {
|
cel->prjfwd = tnxfwd;
|
cel->prjfwd = tnxfwd;
|
cel->prjrev = tnxrev;
|
cel->prjrev = tnxrev;
|
theta0 = 90.0;
|
theta0 = 90.0;
|
} else {
|
} else {
|
/* Unrecognized projection code. */
|
/* Unrecognized projection code. */
|
return 1;
|
return 1;
|
}
|
}
|
|
|
/* Set default for native longitude of the celestial pole? */
|
/* Set default for native longitude of the celestial pole? */
|
dophip = (cel->ref[2] == 999.0);
|
dophip = (cel->ref[2] == 999.0);
|
|
|
/* Compute celestial coordinates of the native pole. */
|
/* Compute celestial coordinates of the native pole. */
|
if (theta0 == 90.0) {
|
if (theta0 == 90.0) {
|
/* Reference point is at the native pole. */
|
/* Reference point is at the native pole. */
|
|
|
if (dophip) {
|
if (dophip) {
|
/* Set default for longitude of the celestial pole. */
|
/* Set default for longitude of the celestial pole. */
|
cel->ref[2] = 180.0;
|
cel->ref[2] = 180.0;
|
}
|
}
|
|
|
latp = cel->ref[1];
|
latp = cel->ref[1];
|
cel->ref[3] = latp;
|
cel->ref[3] = latp;
|
|
|
cel->euler[0] = cel->ref[0];
|
cel->euler[0] = cel->ref[0];
|
cel->euler[1] = 90.0 - latp;
|
cel->euler[1] = 90.0 - latp;
|
} else {
|
} else {
|
/* Reference point away from the native pole. */
|
/* Reference point away from the native pole. */
|
|
|
/* Set default for longitude of the celestial pole. */
|
/* Set default for longitude of the celestial pole. */
|
if (dophip) {
|
if (dophip) {
|
cel->ref[2] = (cel->ref[1] < theta0) ? 180.0 : 0.0;
|
cel->ref[2] = (cel->ref[1] < theta0) ? 180.0 : 0.0;
|
}
|
}
|
|
|
clat0 = wcs_cosd(cel->ref[1]);
|
clat0 = wcs_cosd(cel->ref[1]);
|
slat0 = wcs_sind(cel->ref[1]);
|
slat0 = wcs_sind(cel->ref[1]);
|
cphip = wcs_cosd(cel->ref[2]);
|
cphip = wcs_cosd(cel->ref[2]);
|
sphip = wcs_sind(cel->ref[2]);
|
sphip = wcs_sind(cel->ref[2]);
|
cthe0 = wcs_cosd(theta0);
|
cthe0 = wcs_cosd(theta0);
|
sthe0 = wcs_sind(theta0);
|
sthe0 = wcs_sind(theta0);
|
|
|
x = cthe0*cphip;
|
x = cthe0*cphip;
|
y = sthe0;
|
y = sthe0;
|
z = sqrt(x*x + y*y);
|
z = sqrt(x*x + y*y);
|
if (z == 0.0) {
|
if (z == 0.0) {
|
if (slat0 != 0.0) {
|
if (slat0 != 0.0) {
|
return 1;
|
return 1;
|
}
|
}
|
|
|
/* latp determined by LATPOLE in this case. */
|
/* latp determined by LATPOLE in this case. */
|
latp = cel->ref[3];
|
latp = cel->ref[3];
|
} else {
|
} else {
|
if (fabs(slat0/z) > 1.0) {
|
if (fabs(slat0/z) > 1.0) {
|
return 1;
|
return 1;
|
}
|
}
|
|
|
u = wcs_atan2d(y,x);
|
u = wcs_atan2d(y,x);
|
v = wcs_acosd(slat0/z);
|
v = wcs_acosd(slat0/z);
|
|
|
latp1 = u + v;
|
latp1 = u + v;
|
if (latp1 > 180.0) {
|
if (latp1 > 180.0) {
|
latp1 -= 360.0;
|
latp1 -= 360.0;
|
} else if (latp1 < -180.0) {
|
} else if (latp1 < -180.0) {
|
latp1 += 360.0;
|
latp1 += 360.0;
|
}
|
}
|
|
|
latp2 = u - v;
|
latp2 = u - v;
|
if (latp2 > 180.0) {
|
if (latp2 > 180.0) {
|
latp2 -= 360.0;
|
latp2 -= 360.0;
|
} else if (latp2 < -180.0) {
|
} else if (latp2 < -180.0) {
|
latp2 += 360.0;
|
latp2 += 360.0;
|
}
|
}
|
|
|
if (fabs(cel->ref[3]-latp1) < fabs(cel->ref[3]-latp2)) {
|
if (fabs(cel->ref[3]-latp1) < fabs(cel->ref[3]-latp2)) {
|
if (fabs(latp1) < 90.0+tol) {
|
if (fabs(latp1) < 90.0+tol) {
|
latp = latp1;
|
latp = latp1;
|
} else {
|
} else {
|
latp = latp2;
|
latp = latp2;
|
}
|
}
|
} else {
|
} else {
|
if (fabs(latp2) < 90.0+tol) {
|
if (fabs(latp2) < 90.0+tol) {
|
latp = latp2;
|
latp = latp2;
|
} else {
|
} else {
|
latp = latp1;
|
latp = latp1;
|
}
|
}
|
}
|
}
|
|
|
cel->ref[3] = latp;
|
cel->ref[3] = latp;
|
}
|
}
|
|
|
cel->euler[1] = 90.0 - latp;
|
cel->euler[1] = 90.0 - latp;
|
|
|
z = wcs_cosd(latp)*clat0;
|
z = wcs_cosd(latp)*clat0;
|
if (fabs(z) < tol) {
|
if (fabs(z) < tol) {
|
if (fabs(clat0) < tol) {
|
if (fabs(clat0) < tol) {
|
/* Celestial pole at the reference point. */
|
/* Celestial pole at the reference point. */
|
cel->euler[0] = cel->ref[0];
|
cel->euler[0] = cel->ref[0];
|
cel->euler[1] = 90.0 - theta0;
|
cel->euler[1] = 90.0 - theta0;
|
} else if (latp > 0.0) {
|
} else if (latp > 0.0) {
|
/* Celestial pole at the native north pole.*/
|
/* Celestial pole at the native north pole.*/
|
cel->euler[0] = cel->ref[0] + cel->ref[2] - 180.0;
|
cel->euler[0] = cel->ref[0] + cel->ref[2] - 180.0;
|
cel->euler[1] = 0.0;
|
cel->euler[1] = 0.0;
|
} else if (latp < 0.0) {
|
} else if (latp < 0.0) {
|
/* Celestial pole at the native south pole. */
|
/* Celestial pole at the native south pole. */
|
cel->euler[0] = cel->ref[0] - cel->ref[2];
|
cel->euler[0] = cel->ref[0] - cel->ref[2];
|
cel->euler[1] = 180.0;
|
cel->euler[1] = 180.0;
|
}
|
}
|
} else {
|
} else {
|
x = (sthe0 - wcs_sind(latp)*slat0)/z;
|
x = (sthe0 - wcs_sind(latp)*slat0)/z;
|
y = sphip*cthe0/clat0;
|
y = sphip*cthe0/clat0;
|
if (x == 0.0 && y == 0.0) {
|
if (x == 0.0 && y == 0.0) {
|
return 1;
|
return 1;
|
}
|
}
|
cel->euler[0] = cel->ref[0] - wcs_atan2d(y,x);
|
cel->euler[0] = cel->ref[0] - wcs_atan2d(y,x);
|
}
|
}
|
|
|
/* Make euler[0] the same sign as ref[0]. */
|
/* Make euler[0] the same sign as ref[0]. */
|
if (cel->ref[0] >= 0.0) {
|
if (cel->ref[0] >= 0.0) {
|
if (cel->euler[0] < 0.0) cel->euler[0] += 360.0;
|
if (cel->euler[0] < 0.0) cel->euler[0] += 360.0;
|
} else {
|
} else {
|
if (cel->euler[0] > 0.0) cel->euler[0] -= 360.0;
|
if (cel->euler[0] > 0.0) cel->euler[0] -= 360.0;
|
}
|
}
|
}
|
}
|
|
|
cel->euler[2] = cel->ref[2];
|
cel->euler[2] = cel->ref[2];
|
cel->euler[3] = wcs_cosd(cel->euler[1]);
|
cel->euler[3] = wcs_cosd(cel->euler[1]);
|
cel->euler[4] = wcs_sind(cel->euler[1]);
|
cel->euler[4] = wcs_sind(cel->euler[1]);
|
cel->flag = CELSET;
|
cel->flag = CELSET;
|
|
|
/* Check for ill-conditioned parameters. */
|
/* Check for ill-conditioned parameters. */
|
if (fabs(latp) > 90.0+tol) {
|
if (fabs(latp) > 90.0+tol) {
|
return 2;
|
return 2;
|
}
|
}
|
|
|
return 0;
|
return 0;
|
}
|
}
|
|
|
/*--------------------------------------------------------------------------*/
|
/*--------------------------------------------------------------------------*/
|
|
|
int celfwd(pcode, lng, lat, cel, phi, theta, prj, x, y)
|
int celfwd(pcode, lng, lat, cel, phi, theta, prj, x, y)
|
|
|
const char pcode[4];
|
const char pcode[4];
|
const double lng, lat;
|
const double lng, lat;
|
struct celprm *cel;
|
struct celprm *cel;
|
double *phi, *theta;
|
double *phi, *theta;
|
struct prjprm *prj;
|
struct prjprm *prj;
|
double *x, *y;
|
double *x, *y;
|
|
|
{
|
{
|
int err;
|
int err;
|
|
|
if (cel->flag != CELSET) {
|
if (cel->flag != CELSET) {
|
if (celset(pcode, cel, prj)) return 1;
|
if (celset(pcode, cel, prj)) return 1;
|
}
|
}
|
|
|
/* Compute native coordinates. */
|
/* Compute native coordinates. */
|
sphfwd(lng, lat, cel->euler, phi, theta);
|
sphfwd(lng, lat, cel->euler, phi, theta);
|
|
|
/* Apply forward projection. */
|
/* Apply forward projection. */
|
if ((err = cel->prjfwd(*phi, *theta, prj, x, y))) {
|
if ((err = cel->prjfwd(*phi, *theta, prj, x, y))) {
|
return err == 1 ? 2 : 3;
|
return err == 1 ? 2 : 3;
|
}
|
}
|
|
|
return 0;
|
return 0;
|
}
|
}
|
|
|
/*--------------------------------------------------------------------------*/
|
/*--------------------------------------------------------------------------*/
|
|
|
int celrev(pcode, x, y, prj, phi, theta, cel, lng, lat)
|
int celrev(pcode, x, y, prj, phi, theta, cel, lng, lat)
|
|
|
const char pcode[4];
|
const char pcode[4];
|
const double x, y;
|
const double x, y;
|
struct prjprm *prj;
|
struct prjprm *prj;
|
double *phi, *theta;
|
double *phi, *theta;
|
struct celprm *cel;
|
struct celprm *cel;
|
double *lng, *lat;
|
double *lng, *lat;
|
|
|
{
|
{
|
int err;
|
int err;
|
|
|
if (cel->flag != CELSET) {
|
if (cel->flag != CELSET) {
|
if(celset(pcode, cel, prj)) return 1;
|
if(celset(pcode, cel, prj)) return 1;
|
}
|
}
|
|
|
/* Apply reverse projection. */
|
/* Apply reverse projection. */
|
if ((err = cel->prjrev(x, y, prj, phi, theta))) {
|
if ((err = cel->prjrev(x, y, prj, phi, theta))) {
|
return err == 1 ? 2 : 3;
|
return err == 1 ? 2 : 3;
|
}
|
}
|
if (fabs(*phi)>180.0 || fabs(*theta)>90.0)
|
if (fabs(*phi)>180.0 || fabs(*theta)>90.0)
|
return 2;
|
return 2;
|
|
|
/* Compute native coordinates. */
|
/* Compute native coordinates. */
|
sphrev(*phi, *theta, cel->euler, lng, lat);
|
sphrev(*phi, *theta, cel->euler, lng, lat);
|
|
|
return 0;
|
return 0;
|
}
|
}
|
|
|
