Newer
Older
/*
* proj.c
*
* Spherical map projections.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* This file part of: AstrOmatic WCS library
*
* Copyright: (C) 2000-2010 IAP/CNRS/UPMC
* (C) 1995-1999 Mark Calabretta
*
* Authors: Emmanuel Bertin (this version)
* Mark Calabretta (original version)
*
* Licenses: GNU General Public License
*
* AstrOmatic software is free software: you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
* AstrOmatic software 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 General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with AstrOmatic software.
* If not, see <http://www.gnu.org/licenses/>.
*
* Last modified: 10/10/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
/*============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, 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., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* Correspondence concerning WCSLIB may be directed to:
* Internet email: mcalabre@atnf.csiro.au
* Postal address: Dr. Mark Calabretta,
* Australia Telescope National Facility,
* P.O. Box 76,
* Epping, NSW, 2121,
* AUSTRALIA
*
*=============================================================================
*
* C implementation of the spherical map projections recognized by the FITS
* "World Coordinate System" (WCS) convention.
*
* Summary of routines
* -------------------
* Each projection is implemented via separate functions for the forward,
* *fwd(), and reverse, *rev(), transformation.
*
* Initialization routines, *set(), compute intermediate values from the
* projection parameters but need not be called explicitly - see the
* explanation of prj.flag below.
*
* azpset azpfwd azprev AZP: zenithal/azimuthal perspective
* tanset tanfwd tanrev TAN: gnomonic
* sinset sinfwd sinrev SIN: orthographic/synthesis
* stgset stgfwd stgrev STG: stereographic
* arcset arcfwd arcrev ARC: zenithal/azimuthal equidistant
* zpnset zpnfwd zpnrev ZPN: zenithal/azimuthal polynomial
* zeaset zeafwd zearev ZEA: zenithal/azimuthal equal area
* airset airfwd airrev AIR: Airy
* cypset cypfwd cyprev CYP: cylindrical perspective
* carset carfwd carrev CAR: Cartesian
* merset merfwd merrev MER: Mercator
* ceaset ceafwd cearev CEA: cylindrical equal area
* copset copfwd coprev COP: conic perspective
* codset codfwd codrev COD: conic equidistant
* coeset coefwd coerev COE: conic equal area
* cooset coofwd coorev COO: conic orthomorphic
* bonset bonfwd bonrev BON: Bonne
* pcoset pcofwd pcorev PCO: polyconic
* glsset glsfwd glsrev GLS: Sanson-Flamsteed (global sinusoidal)
* parset parfwd parrev PAR: parabolic
* aitset aitfwd aitrev AIT: Hammer-Aitoff
* molset molfwd molrev MOL: Mollweide
* cscset cscfwd cscrev CSC: COBE quadrilateralized spherical cube
* qscset qscfwd qscrev QSC: quadrilateralized spherical cube
* tscset tscfwd tscrev TSC: tangential spherical cube
* tnxset tnxfwd tnxrev TNX: IRAF's gnomonic
*
*
* Initialization routine; *set()
* ------------------------------
* Initializes members of a prjprm data structure which hold intermediate
* values. Note that this routine need not be called directly; it will be
* invoked by prjfwd() and prjrev() if the "flag" structure member is
* anything other than a predefined magic value.
*
* Given and/or returned:
* prj prjprm* Projection parameters (see below).
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid projection parameters.
*
* Forward transformation; *fwd()
* -----------------------------
* Compute (x,y) coordinates in the plane of projection from native spherical
* coordinates (phi,theta).
*
* Given:
* phi, const double
* theta Longitude and latitude of the projected point in
* native spherical coordinates, in degrees.
*
* Given and returned:
* prj prjprm* Projection parameters (see below).
*
* Returned:
* x,y double* Projected coordinates.
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid projection parameters.
* 2: Invalid value of (phi,theta).
*
* Reverse transformation; *rev()
* -----------------------------
* Compute native spherical coordinates (phi,theta) from (x,y) coordinates in
* the plane of projection.
*
* Given:
* x,y const double
* Projected coordinates.
*
* Given and returned:
* prj prjprm* Projection parameters (see below).
*
* Returned:
* phi, double* Longitude and latitude of the projected point in
* theta native spherical coordinates, in degrees.
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid projection parameters.
* 2: Invalid value of (x,y).
*
* Projection parameters
* ---------------------
* The prjprm struct consists of the following:
*
* int flag
* This flag must be set to zero whenever any of p[10] or r0 are set
* or changed. This signals the initialization routine to recompute
* intermediaries. flag may also be set to -1 to disable strict bounds
* checking for the AZP, TAN, SIN, ZPN, and COP projections.
* double r0
* r0; The radius of the generating sphere for the projection, a linear
* scaling parameter. If this is zero, it will be reset to the default
* value of 180/pi (the value for FITS WCS).
* double p[10]
* The first 10 elements contain projection parameters which correspond
* to the PROJPn keywords in FITS, so p[0] is PROJP0, and p[9] is
* PROJP9. Many projections use p[1] (PROJP1) and some also use p[2]
* (PROJP2). ZPN is the only projection which uses any of the others.
*
* The remaining members of the prjprm struct are maintained by the
* initialization routines and should not be modified. This is done for the
* sake of efficiency and to allow an arbitrary number of contexts to be
* maintained simultaneously.
*
* int n
* double w[10]
* Intermediate values derived from the projection parameters.
*
* 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.
*
Loading
Loading full blame…