1 /*============================================================================
2   WCSLIB 7.3 - an implementation of the FITS WCS standard.
3   Copyright (C) 1995-2020, Mark Calabretta
4 
5   This file is part of WCSLIB.
6 
7   WCSLIB is free software: you can redistribute it and/or modify it under the
8   terms of the GNU Lesser General Public License as published by the Free
9   Software Foundation, either version 3 of the License, or (at your option)
10   any later version.
11 
12   WCSLIB is distributed in the hope that it will be useful, but WITHOUT ANY
13   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14   FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
15   more details.
16 
17   You should have received a copy of the GNU Lesser General Public License
18   along with WCSLIB.  If not, see http://www.gnu.org/licenses.
19 
20   Direct correspondence concerning WCSLIB to mark@calabretta.id.au
21 
22   Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
23   http://www.atnf.csiro.au/people/Mark.Calabretta
24   $Id: wcs.h,v 7.3.1.2 2020/08/17 12:10:44 mcalabre Exp mcalabre $
25 *=============================================================================
26 *
27 * WCSLIB 7.3 - C routines that implement the FITS World Coordinate System
28 * (WCS) standard.  Refer to the README file provided with WCSLIB for an
29 * overview of the library.
30 *
31 *
32 * Summary of the wcs routines
33 * ---------------------------
34 * Routines in this suite implement the FITS World Coordinate System (WCS)
35 * standard which defines methods to be used for computing world coordinates
36 * from image pixel coordinates, and vice versa.  The standard, and proposed
37 * extensions for handling distortions, are described in
38 *
39 =   "Representations of world coordinates in FITS",
40 =   Greisen, E.W., & Calabretta, M.R. 2002, A&A, 395, 1061 (WCS Paper I)
41 =
42 =   "Representations of celestial coordinates in FITS",
43 =   Calabretta, M.R., & Greisen, E.W. 2002, A&A, 395, 1077 (WCS Paper II)
44 =
45 =   "Representations of spectral coordinates in FITS",
46 =   Greisen, E.W., Calabretta, M.R., Valdes, F.G., & Allen, S.L.
47 =   2006, A&A, 446, 747 (WCS Paper III)
48 =
49 =   "Representations of distortions in FITS world coordinate systems",
50 =   Calabretta, M.R. et al. (WCS Paper IV, draft dated 2004/04/22),
51 =   available from http://www.atnf.csiro.au/people/Mark.Calabretta
52 =
53 =   "Mapping on the HEALPix grid",
54 =   Calabretta, M.R., & Roukema, B.F. 2007, MNRAS, 381, 865 (WCS Paper V)
55 =
56 =   "Representing the 'Butterfly' Projection in FITS -- Projection Code XPH",
57 =   Calabretta, M.R., & Lowe, S.R. 2013, PASA, 30, e050 (WCS Paper VI)
58 =
59 =   "Representations of time coordinates in FITS -
60 =    Time and relative dimension in space",
61 =   Rots, A.H., Bunclark, P.S., Calabretta, M.R., Allen, S.L.,
62 =   Manchester, R.N., & Thompson, W.T. 2015, A&A, 574, A36 (WCS Paper VII)
63 *
64 * These routines are based on the wcsprm struct which contains all information
65 * needed for the computations.  The struct contains some members that must be
66 * set by the user, and others that are maintained by these routines, somewhat
67 * like a C++ class but with no encapsulation.
68 *
69 * wcsnpv(), wcsnps(), wcsini(), wcsinit(), wcssub(), and wcsfree() are
70 * provided to manage the wcsprm struct and another, wcsprt(), prints its
71 * contents.  Refer to the description of the wcsprm struct for an explanation
72 * of the anticipated usage of these routines.  wcscopy(), which does a deep
73 * copy of one wcsprm struct to another, is defined as a preprocessor macro
74 * function that invokes wcssub().
75 *
76 * wcsperr() prints the error message(s) (if any) stored in a wcsprm struct,
77 * and the linprm, celprm, prjprm, spcprm, and tabprm structs that it contains.
78 *
79 * A setup routine, wcsset(), computes intermediate values in the wcsprm struct
80 * from parameters in it that were supplied by the user.  The struct always
81 * needs to be set up by wcsset() but this need not be called explicitly -
82 * refer to the explanation of wcsprm::flag.
83 *
84 * wcsp2s() and wcss2p() implement the WCS world coordinate transformations.
85 * In fact, they are high level driver routines for the WCS linear,
86 * logarithmic, celestial, spectral and tabular transformation routines
87 * described in lin.h, log.h, cel.h, spc.h and tab.h.
88 *
89 * Given either the celestial longitude or latitude plus an element of the
90 * pixel coordinate a hybrid routine, wcsmix(), iteratively solves for the
91 * unknown elements.
92 *
93 * wcssptr() translates the spectral axis in a wcsprm struct.  For example, a
94 * 'FREQ' axis may be translated into 'ZOPT-F2W' and vice versa.
95 *
96 * wcslib_version() returns the WCSLIB version number.
97 *
98 * Quadcube projections:
99 * ---------------------
100 *   The quadcube projections (TSC, CSC, QSC) may be represented in FITS in
101 *   either of two ways:
102 *
103 *     a: The six faces may be laid out in one plane and numbered as follows:
104 *
105 =                                 0
106 =
107 =                        4  3  2  1  4  3  2
108 =
109 =                                 5
110 *
111 *        Faces 2, 3 and 4 may appear on one side or the other (or both).  The
112 *        world-to-pixel routines map faces 2, 3 and 4 to the left but the
113 *        pixel-to-world routines accept them on either side.
114 *
115 *     b: The "COBE" convention in which the six faces are stored in a
116 *        three-dimensional structure using a CUBEFACE axis indexed from
117 *        0 to 5 as above.
118 *
119 *   These routines support both methods; wcsset() determines which is being
120 *   used by the presence or absence of a CUBEFACE axis in ctype[].  wcsp2s()
121 *   and wcss2p() translate the CUBEFACE axis representation to the single
122 *   plane representation understood by the lower-level WCSLIB projection
123 *   routines.
124 *
125 *
126 * wcsnpv() - Memory allocation for PVi_ma
127 * ---------------------------------------
128 * wcsnpv() sets or gets the value of NPVMAX (default 64).  This global
129 * variable controls the number of pvcard structs, for holding PVi_ma
130 * keyvalues, that wcsini() should allocate space for.  It is also used by
131 * wcsinit() as the default value of npvmax.
132 *
133 * PLEASE NOTE: This function is not thread-safe.
134 *
135 * Given:
136 *   n         int       Value of NPVMAX; ignored if < 0.  Use a value less
137 *                       than zero to get the current value.
138 *
139 * Function return value:
140 *             int       Current value of NPVMAX.
141 *
142 *
143 * wcsnps() - Memory allocation for PSi_ma
144 * ---------------------------------------
145 * wcsnps() sets or gets the value of NPSMAX (default 8).  This global variable
146 * controls the number of pscard structs, for holding PSi_ma keyvalues, that
147 * wcsini() should allocate space for.  It is also used by wcsinit() as the
148 * default value of npsmax.
149 *
150 * PLEASE NOTE: This function is not thread-safe.
151 *
152 * Given:
153 *   n         int       Value of NPSMAX; ignored if < 0.  Use a value less
154 *                       than zero to get the current value.
155 *
156 * Function return value:
157 *             int       Current value of NPSMAX.
158 *
159 *
160 * wcsini() - Default constructor for the wcsprm struct
161 * ----------------------------------------------------
162 * wcsini() is a thin wrapper on wcsinit().  It invokes it with npvmax,
163 * npsmax, and ndpmax set to -1 which causes it to use the values of the
164 * global variables NDPMAX, NPSMAX, and NDPMAX.  It is thereby potentially
165 * thread-unsafe if these variables are altered dynamically via wcsnpv(),
166 * wcsnps(), and disndp().  Use wcsinit() for a thread-safe alternative in
167 * this case.
168 *
169 *
170 * wcsinit() - Default constructor for the wcsprm struct
171 * -----------------------------------------------------
172 * wcsinit() optionally allocates memory for arrays in a wcsprm struct and sets
173 * all members of the struct to default values.
174 *
175 * PLEASE NOTE: every wcsprm struct should be initialized by wcsinit(),
176 * possibly repeatedly.  On the first invokation, and only the first
177 * invokation, wcsprm::flag must be set to -1 to initialize memory management,
178 * regardless of whether wcsinit() will actually be used to allocate memory.
179 *
180 * Given:
181 *   alloc     int       If true, allocate memory unconditionally for the
182 *                       crpix, etc. arrays.  Please note that memory is never
183 *                       allocated by wcsinit() for the auxprm, tabprm, nor
184 *                       wtbarr structs.
185 *
186 *                       If false, it is assumed that pointers to these arrays
187 *                       have been set by the user except if they are null
188 *                       pointers in which case memory will be allocated for
189 *                       them regardless.  (In other words, setting alloc true
190 *                       saves having to initalize these pointers to zero.)
191 *
192 *   naxis     int       The number of world coordinate axes.  This is used to
193 *                       determine the length of the various wcsprm vectors and
194 *                       matrices and therefore the amount of memory to
195 *                       allocate for them.
196 *
197 * Given and returned:
198 *   wcs       struct wcsprm*
199 *                       Coordinate transformation parameters.
200 *
201 *                       Note that, in order to initialize memory management,
202 *                       wcsprm::flag should be set to -1 when wcs is
203 *                       initialized for the first time (memory leaks may
204 *                       result if it had already been initialized).
205 *
206 * Given:
207 *   npvmax    int       The number of PVi_ma keywords to allocate space for.
208 *                       If set to -1, the value of the global variable NPVMAX
209 *                       will be used.  This is potentially thread-unsafe if
210 *                       wcsnpv() is being used dynamically to alter its value.
211 *
212 *   npsmax    int       The number of PSi_ma keywords to allocate space for.
213 *                       If set to -1, the value of the global variable NPSMAX
214 *                       will be used.  This is potentially thread-unsafe if
215 *                       wcsnps() is being used dynamically to alter its value.
216 *
217 *   ndpmax    int       The number of DPja or DQia keywords to allocate space
218 *                       for.  If set to -1, the value of the global variable
219 *                       NDPMAX will be used.  This is potentially
220 *                       thread-unsafe if disndp() is being used dynamically to
221 *                       alter its value.
222 *
223 * Function return value:
224 *             int       Status return value:
225 *                         0: Success.
226 *                         1: Null wcsprm pointer passed.
227 *                         2: Memory allocation failed.
228 *
229 *                       For returns > 1, a detailed error message is set in
230 *                       wcsprm::err if enabled, see wcserr_enable().
231 *
232 *
233 * wcsauxi() - Default constructor for the auxprm struct
234 * -----------------------------------------------------
235 * wcsauxi() optionally allocates memory for an auxprm struct, attaches it to
236 * wcsprm, and sets all members of the struct to default values.
237 *
238 * Given:
239 *   alloc     int       If true, allocate memory unconditionally for the
240 *                       auxprm struct.
241 *
242 *                       If false, it is assumed that wcsprm::aux has already
243 *                       been set to point to an auxprm struct, in which case
244 *                       the user is responsible for managing that memory.
245 *                       However, if wcsprm::aux is a null pointer, memory will
246 *                       be allocated regardless.  (In other words, setting
247 *                       alloc true saves having to initalize the pointer to
248 *                       zero.)
249 *
250 * Given and returned:
251 *   wcs       struct wcsprm*
252 *                       Coordinate transformation parameters.
253 *
254 * Function return value:
255 *             int       Status return value:
256 *                         0: Success.
257 *                         1: Null wcsprm pointer passed.
258 *                         2: Memory allocation failed.
259 *
260 *
261 * wcssub() - Subimage extraction routine for the wcsprm struct
262 * ------------------------------------------------------------
263 * wcssub() extracts the coordinate description for a subimage from a wcsprm
264 * struct.  It does a deep copy, using wcsinit() to allocate memory for its
265 * arrays if required.  Only the "information to be provided" part of the
266 * struct is extracted.  Consequently, wcsset() need not have been, and won't
267 * be invoked on the struct from which the subimage is extracted.  A call to
268 * wcsset() is required to set up the subimage struct.
269 *
270 * The world coordinate system of the subimage must be separable in the sense
271 * that the world coordinates at any point in the subimage must depend only on
272 * the pixel coordinates of the axes extracted.  In practice, this means that
273 * the linear transformation matrix of the original image must not contain
274 * non-zero off-diagonal terms that associate any of the subimage axes with any
275 * of the non-subimage axes.  Likewise, if any distortions are associated with
276 * the subimage axes, they must not depend on any of the axes that are not
277 * being extracted.
278 *
279 * Note that while the required elements of the tabprm array are extracted, the
280 * wtbarr array is not.  (Thus it is not appropriate to call wcssub() after
281 * wcstab() but before filling the tabprm structs - refer to wcshdr.h.)
282 *
283 * wcssub() can also add axes to a wcsprm struct.  The new axes will be created
284 * using the defaults set by wcsinit() which produce a simple, unnamed, linear
285 * axis with world coordinate equal to the pixel coordinate.  These default
286 * values can be changed afterwards, before invoking wcsset().
287 *
288 * Given:
289 *   alloc     int       If true, allocate memory for the crpix, etc. arrays in
290 *                       the destination.  Otherwise, it is assumed that
291 *                       pointers to these arrays have been set by the user
292 *                       except if they are null pointers in which case memory
293 *                       will be allocated for them regardless.
294 *
295 *   wcssrc    const struct wcsprm*
296 *                       Struct to extract from.
297 *
298 * Given and returned:
299 *   nsub      int*
300 *   axes      int[]     Vector of length *nsub containing the image axis
301 *                       numbers (1-relative) to extract.  Order is
302 *                       significant; axes[0] is the axis number of the input
303 *                       image that corresponds to the first axis in the
304 *                       subimage, etc.
305 *
306 *                       Use an axis number of 0 to create a new axis using
307 *                       the defaults set by wcsinit().  They can be changed
308 *                       later.
309 *
310 *                       nsub (the pointer) may be set to zero, and so also may
311 *                       *nsub, which is interpreted to mean all axes in the
312 *                       input image; the number of axes will be returned if
313 *                       nsub != 0x0.  axes itself (the pointer) may be set to
314 *                       zero to indicate the first *nsub axes in their
315 *                       original order.
316 *
317 *                       Set both nsub (or *nsub) and axes to zero to do a deep
318 *                       copy of one wcsprm struct to another.
319 *
320 *                       Subimage extraction by coordinate axis type may be
321 *                       done by setting the elements of axes[] to the
322 *                       following special preprocessor macro values:
323 *
324 *                         WCSSUB_LONGITUDE: Celestial longitude.
325 *                         WCSSUB_LATITUDE:  Celestial latitude.
326 *                         WCSSUB_CUBEFACE:  Quadcube CUBEFACE axis.
327 *                         WCSSUB_SPECTRAL:  Spectral axis.
328 *                         WCSSUB_STOKES:    Stokes axis.
329 *
330 *                       Refer to the notes (below) for further usage examples.
331 *
332 *                       On return, *nsub will be set to the number of axes in
333 *                       the subimage; this may be zero if there were no axes
334 *                       of the required type(s) (in which case no memory will
335 *                       be allocated).  axes[] will contain the axis numbers
336 *                       that were extracted, or 0 for newly created axes.  The
337 *                       vector length must be sufficient to contain all axis
338 *                       numbers.  No checks are performed to verify that the
339 *                       coordinate axes are consistent, this is done by
340 *                       wcsset().
341 *
342 *   wcsdst    struct wcsprm*
343 *                       Struct describing the subimage.  wcsprm::flag should
344 *                       be set to -1 if wcsdst was not previously initialized
345 *                       (memory leaks may result if it was previously
346 *                       initialized).
347 *
348 * Function return value:
349 *             int       Status return value:
350 *                         0: Success.
351 *                         1: Null wcsprm pointer passed.
352 *                         2: Memory allocation failed.
353 *                        12: Invalid subimage specification.
354 *                        13: Non-separable subimage coordinate system.
355 *
356 *                       For returns > 1, a detailed error message is set in
357 *                       wcsprm::err if enabled, see wcserr_enable().
358 *
359 * Notes:
360 *   Combinations of subimage axes of particular types may be extracted in the
361 *   same order as they occur in the input image by combining preprocessor
362 *   codes, for example
363 *
364 =     *nsub = 1;
365 =     axes[0] = WCSSUB_LONGITUDE | WCSSUB_LATITUDE | WCSSUB_SPECTRAL;
366 *
367 *   would extract the longitude, latitude, and spectral axes in the same order
368 *   as the input image.  If one of each were present, *nsub = 3 would be
369 *   returned.
370 *
371 *   For convenience, WCSSUB_CELESTIAL is defined as the combination
372 *   WCSSUB_LONGITUDE | WCSSUB_LATITUDE | WCSSUB_CUBEFACE.
373 *
374 *   The codes may also be negated to extract all but the types specified, for
375 *   example
376 *
377 =     *nsub = 4;
378 =     axes[0] = WCSSUB_LONGITUDE;
379 =     axes[1] = WCSSUB_LATITUDE;
380 =     axes[2] = WCSSUB_CUBEFACE;
381 =     axes[3] = -(WCSSUB_SPECTRAL | WCSSUB_STOKES);
382 *
383 *   The last of these specifies all axis types other than spectral or Stokes.
384 *   Extraction is done in the order specified by axes[] a longitude axis (if
385 *   present) would be extracted first (via axes[0]) and not subsequently (via
386 *   axes[3]).  Likewise for the latitude and cubeface axes in this example.
387 *
388 *   From the foregoing, it is apparent that the value of *nsub returned may be
389 *   less than or greater than that given.  However, it will never exceed the
390 *   number of axes in the input image (plus the number of newly-created axes
391 *   if any were specified on input).
392 *
393 *
394 * wcscompare() - Compare two wcsprm structs for equality
395 * ------------------------------------------------------
396 * wcscompare() compares two wcsprm structs for equality.
397 *
398 * Given:
399 *   cmp       int       A bit field controlling the strictness of the
400 *                       comparison.  When 0, all fields must be identical.
401 *
402 *                       The following constants may be or'ed together to
403 *                       relax the comparison:
404 *                         WCSCOMPARE_ANCILLARY: Ignore ancillary keywords
405 *                           that don't change the WCS transformation, such
406 *                           as DATE-OBS or EQUINOX.
407 *                         WCSCOMPARE_TILING: Ignore integral differences in
408 *                           CRPIXja.  This is the 'tiling' condition, where
409 *                           two WCSes cover different regions of the same
410 *                           map projection and align on the same map grid.
411 *                         WCSCOMPARE_CRPIX: Ignore any differences at all in
412 *                           CRPIXja.  The two WCSes cover different regions
413 *                           of the same map projection but may not align on
414 *                           the same map grid.  Overrides WCSCOMPARE_TILING.
415 *
416 *   tol       double    Tolerance for comparison of floating-point values.
417 *                       For example, for tol == 1e-6, all floating-point
418 *                       values in the structs must be equal to the first 6
419 *                       decimal places.  A value of 0 implies exact equality.
420 *
421 *   wcs1      const struct wcsprm*
422 *                       The first wcsprm struct to compare.
423 *
424 *   wcs2      const struct wcsprm*
425 *                       The second wcsprm struct to compare.
426 *
427 * Returned:
428 *   equal     int*      Non-zero when the given structs are equal.
429 *
430 * Function return value:
431 *             int       Status return value:
432 *                         0: Success.
433 *                         1: Null pointer passed.
434 *
435 *
436 * wcscopy() macro - Copy routine for the wcsprm struct
437 * ----------------------------------------------------
438 * wcscopy() does a deep copy of one wcsprm struct to another.  As of
439 * WCSLIB 3.6, it is implemented as a preprocessor macro that invokes
440 * wcssub() with the nsub and axes pointers both set to zero.
441 *
442 *
443 * wcsfree() - Destructor for the wcsprm struct
444 * --------------------------------------------
445 * wcsfree() frees memory allocated for the wcsprm arrays by wcsinit() and/or
446 * wcsset().  wcsinit() records the memory it allocates and wcsfree() will only
447 * attempt to free this.
448 *
449 * PLEASE NOTE: wcsfree() must not be invoked on a wcsprm struct that was not
450 * initialized by wcsinit().
451 *
452 * Returned:
453 *   wcs       struct wcsprm*
454 *                       Coordinate transformation parameters.
455 *
456 * Function return value:
457 *             int       Status return value:
458 *                         0: Success.
459 *                         1: Null wcsprm pointer passed.
460 *
461 *
462 * wcsprt() - Print routine for the wcsprm struct
463 * ----------------------------------------------
464 * wcsprt() prints the contents of a wcsprm struct using wcsprintf().  Mainly
465 * intended for diagnostic purposes.
466 *
467 * Given:
468 *   wcs       const struct wcsprm*
469 *                       Coordinate transformation parameters.
470 *
471 * Function return value:
472 *             int       Status return value:
473 *                         0: Success.
474 *                         1: Null wcsprm pointer passed.
475 *
476 *
477 * wcsperr() - Print error messages from a wcsprm struct
478 * -----------------------------------------------------
479 * wcsperr() prints the error message(s), if any, stored in a wcsprm struct,
480 * and the linprm, celprm, prjprm, spcprm, and tabprm structs that it contains.
481 * If there are no errors then nothing is printed.  It uses wcserr_prt(), q.v.
482 *
483 * Given:
484 *   wcs       const struct wcsprm*
485 *                       Coordinate transformation parameters.
486 *
487 *   prefix    const char *
488 *                       If non-NULL, each output line will be prefixed with
489 *                       this string.
490 *
491 * Function return value:
492 *             int       Status return value:
493 *                         0: Success.
494 *                         1: Null wcsprm pointer passed.
495 *
496 *
497 * wcsbchk() - Enable/disable bounds checking
498 * ------------------------------------------
499 * wcsbchk() is used to control bounds checking in the projection routines.
500 * Note that wcsset() always enables bounds checking.  wcsbchk() will invoke
501 * wcsset() on the wcsprm struct beforehand if necessary.
502 *
503 * Given and returned:
504 *   wcs       struct wcsprm*
505 *                       Coordinate transformation parameters.
506 *
507 * Given:
508 *   bounds    int       If bounds&1 then enable strict bounds checking for the
509 *                       spherical-to-Cartesian (s2x) transformation for the
510 *                       AZP, SZP, TAN, SIN, ZPN, and COP projections.
511 *
512 *                       If bounds&2 then enable strict bounds checking for the
513 *                       Cartesian-to-spherical (x2s) transformation for the
514 *                       HPX and XPH projections.
515 *
516 *                       If bounds&4 then enable bounds checking on the native
517 *                       coordinates returned by the Cartesian-to-spherical
518 *                       (x2s) transformations using prjchk().
519 *
520 *                       Zero it to disable all checking.
521 *
522 * Function return value:
523 *             int       Status return value:
524 *                         0: Success.
525 *                         1: Null wcsprm pointer passed.
526 *
527 *
528 * wcsset() - Setup routine for the wcsprm struct
529 * ----------------------------------------------
530 * wcsset() sets up a wcsprm struct according to information supplied within
531 * it (refer to the description of the wcsprm struct).
532 *
533 * wcsset() recognizes the NCP projection and converts it to the equivalent SIN
534 * projection and likewise translates GLS into SFL.  It also translates the
535 * AIPS spectral types ('FREQ-LSR', 'FELO-HEL', etc.), possibly changing the
536 * input header keywords wcsprm::ctype and/or wcsprm::specsys if necessary.
537 *
538 * Note that this routine need not be called directly; it will be invoked by
539 * wcsp2s() and wcss2p() if the wcsprm::flag is anything other than a
540 * predefined magic value.
541 *
542 * Given and returned:
543 *   wcs       struct wcsprm*
544 *                       Coordinate transformation parameters.
545 *
546 * Function return value:
547 *             int       Status return value:
548 *                         0: Success.
549 *                         1: Null wcsprm pointer passed.
550 *                         2: Memory allocation failed.
551 *                         3: Linear transformation matrix is singular.
552 *                         4: Inconsistent or unrecognized coordinate axis
553 *                            types.
554 *                         5: Invalid parameter value.
555 *                         6: Invalid coordinate transformation parameters.
556 *                         7: Ill-conditioned coordinate transformation
557 *                            parameters.
558 *
559 *                       For returns > 1, a detailed error message is set in
560 *                       wcsprm::err if enabled, see wcserr_enable().
561 *
562 * Notes:
563 *   wcsset() always enables strict bounds checking in the projection routines
564 *   (via a call to prjini()).  Use wcsbchk() to modify bounds-checking after
565 *   wcsset() is invoked.
566 *
567 *
568 * wcsp2s() - Pixel-to-world transformation
569 * ----------------------------------------
570 * wcsp2s() transforms pixel coordinates to world coordinates.
571 *
572 * Given and returned:
573 *   wcs       struct wcsprm*
574 *                       Coordinate transformation parameters.
575 *
576 * Given:
577 *   ncoord,
578 *   nelem     int       The number of coordinates, each of vector length
579 *                       nelem but containing wcs.naxis coordinate elements.
580 *                       Thus nelem must equal or exceed the value of the
581 *                       NAXIS keyword unless ncoord == 1, in which case nelem
582 *                       is not used.
583 *
584 *   pixcrd    const double[ncoord][nelem]
585 *                       Array of pixel coordinates.
586 *
587 * Returned:
588 *   imgcrd    double[ncoord][nelem]
589 *                       Array of intermediate world coordinates.  For
590 *                       celestial axes, imgcrd[][wcs.lng] and
591 *                       imgcrd[][wcs.lat] are the projected x-, and
592 *                       y-coordinates in pseudo "degrees".  For spectral
593 *                       axes, imgcrd[][wcs.spec] is the intermediate spectral
594 *                       coordinate, in SI units.
595 *
596 *   phi,theta double[ncoord]
597 *                       Longitude and latitude in the native coordinate system
598 *                       of the projection [deg].
599 *
600 *   world     double[ncoord][nelem]
601 *                       Array of world coordinates.  For celestial axes,
602 *                       world[][wcs.lng] and world[][wcs.lat] are the
603 *                       celestial longitude and latitude [deg].  For
604 *                       spectral axes, imgcrd[][wcs.spec] is the intermediate
605 *                       spectral coordinate, in SI units.
606 *
607 *   stat      int[ncoord]
608 *                       Status return value for each coordinate:
609 *                         0: Success.
610 *                        1+: A bit mask indicating invalid pixel coordinate
611 *                            element(s).
612 *
613 * Function return value:
614 *             int       Status return value:
615 *                         0: Success.
616 *                         1: Null wcsprm pointer passed.
617 *                         2: Memory allocation failed.
618 *                         3: Linear transformation matrix is singular.
619 *                         4: Inconsistent or unrecognized coordinate axis
620 *                            types.
621 *                         5: Invalid parameter value.
622 *                         6: Invalid coordinate transformation parameters.
623 *                         7: Ill-conditioned coordinate transformation
624 *                            parameters.
625 *                         8: One or more of the pixel coordinates were
626 *                            invalid, as indicated by the stat vector.
627 *
628 *                       For returns > 1, a detailed error message is set in
629 *                       wcsprm::err if enabled, see wcserr_enable().
630 *
631 *
632 * wcss2p() - World-to-pixel transformation
633 * ----------------------------------------
634 * wcss2p() transforms world coordinates to pixel coordinates.
635 *
636 * Given and returned:
637 *   wcs       struct wcsprm*
638 *                       Coordinate transformation parameters.
639 *
640 * Given:
641 *   ncoord,
642 *   nelem     int       The number of coordinates, each of vector length nelem
643 *                       but containing wcs.naxis coordinate elements.  Thus
644 *                       nelem must equal or exceed the value of the NAXIS
645 *                       keyword unless ncoord == 1, in which case nelem is not
646 *                       used.
647 *
648 *   world     const double[ncoord][nelem]
649 *                       Array of world coordinates.  For celestial axes,
650 *                       world[][wcs.lng] and world[][wcs.lat] are the
651 *                       celestial longitude and latitude [deg]. For spectral
652 *                       axes, world[][wcs.spec] is the spectral coordinate, in
653 *                       SI units.
654 *
655 * Returned:
656 *   phi,theta double[ncoord]
657 *                       Longitude and latitude in the native coordinate
658 *                       system of the projection [deg].
659 *
660 *   imgcrd    double[ncoord][nelem]
661 *                       Array of intermediate world coordinates.  For
662 *                       celestial axes, imgcrd[][wcs.lng] and
663 *                       imgcrd[][wcs.lat] are the projected x-, and
664 *                       y-coordinates in pseudo "degrees".  For quadcube
665 *                       projections with a CUBEFACE axis the face number is
666 *                       also returned in imgcrd[][wcs.cubeface].  For
667 *                       spectral axes, imgcrd[][wcs.spec] is the intermediate
668 *                       spectral coordinate, in SI units.
669 *
670 *   pixcrd    double[ncoord][nelem]
671 *                       Array of pixel coordinates.
672 *
673 *   stat      int[ncoord]
674 *                       Status return value for each coordinate:
675 *                         0: Success.
676 *                        1+: A bit mask indicating invalid world coordinate
677 *                            element(s).
678 *
679 * Function return value:
680 *             int       Status return value:
681 *                         0: Success.
682 *                         1: Null wcsprm pointer passed.
683 *                         2: Memory allocation failed.
684 *                         3: Linear transformation matrix is singular.
685 *                         4: Inconsistent or unrecognized coordinate axis
686 *                            types.
687 *                         5: Invalid parameter value.
688 *                         6: Invalid coordinate transformation parameters.
689 *                         7: Ill-conditioned coordinate transformation
690 *                            parameters.
691 *                         9: One or more of the world coordinates were
692 *                            invalid, as indicated by the stat vector.
693 *
694 *                       For returns > 1, a detailed error message is set in
695 *                       wcsprm::err if enabled, see wcserr_enable().
696 *
697 *
698 * wcsmix() - Hybrid coordinate transformation
699 * -------------------------------------------
700 * wcsmix(), given either the celestial longitude or latitude plus an element
701 * of the pixel coordinate, solves for the remaining elements by iterating on
702 * the unknown celestial coordinate element using wcss2p().  Refer also to the
703 * notes below.
704 *
705 * Given and returned:
706 *   wcs       struct wcsprm*
707 *                       Indices for the celestial coordinates obtained
708 *                       by parsing the wcsprm::ctype[].
709 *
710 * Given:
711 *   mixpix    int       Which element of the pixel coordinate is given.
712 *
713 *   mixcel    int       Which element of the celestial coordinate is given:
714 *                         1: Celestial longitude is given in
715 *                            world[wcs.lng], latitude returned in
716 *                            world[wcs.lat].
717 *                         2: Celestial latitude is given in
718 *                            world[wcs.lat], longitude returned in
719 *                            world[wcs.lng].
720 *
721 *   vspan     const double[2]
722 *                       Solution interval for the celestial coordinate [deg].
723 *                       The ordering of the two limits is irrelevant.
724 *                       Longitude ranges may be specified with any convenient
725 *                       normalization, for example [-120,+120] is the same as
726 *                       [240,480], except that the solution will be returned
727 *                       with the same normalization, i.e. lie within the
728 *                       interval specified.
729 *
730 *   vstep     const double
731 *                       Step size for solution search [deg].  If zero, a
732 *                       sensible, although perhaps non-optimal default will be
733 *                       used.
734 *
735 *   viter     int       If a solution is not found then the step size will be
736 *                       halved and the search recommenced.  viter controls how
737 *                       many times the step size is halved.  The allowed range
738 *                       is 5 - 10.
739 *
740 * Given and returned:
741 *   world     double[naxis]
742 *                       World coordinate elements.  world[wcs.lng] and
743 *                       world[wcs.lat] are the celestial longitude and
744 *                       latitude [deg].  Which is given and which returned
745 *                       depends on the value of mixcel.  All other elements
746 *                       are given.
747 *
748 * Returned:
749 *   phi,theta double[naxis]
750 *                       Longitude and latitude in the native coordinate
751 *                       system of the projection [deg].
752 *
753 *   imgcrd    double[naxis]
754 *                       Image coordinate elements.  imgcrd[wcs.lng] and
755 *                       imgcrd[wcs.lat] are the projected x-, and
756 *                       y-coordinates in pseudo "degrees".
757 *
758 * Given and returned:
759 *   pixcrd    double[naxis]
760 *                       Pixel coordinate.  The element indicated by mixpix is
761 *                       given and the remaining elements are returned.
762 *
763 * Function return value:
764 *             int       Status return value:
765 *                         0: Success.
766 *                         1: Null wcsprm pointer passed.
767 *                         2: Memory allocation failed.
768 *                         3: Linear transformation matrix is singular.
769 *                         4: Inconsistent or unrecognized coordinate axis
770 *                            types.
771 *                         5: Invalid parameter value.
772 *                         6: Invalid coordinate transformation parameters.
773 *                         7: Ill-conditioned coordinate transformation
774 *                            parameters.
775 *                        10: Invalid world coordinate.
776 *                        11: No solution found in the specified interval.
777 *
778 *                       For returns > 1, a detailed error message is set in
779 *                       wcsprm::err if enabled, see wcserr_enable().
780 *
781 * Notes:
782 *   Initially the specified solution interval is checked to see if it's a
783 *   "crossing" interval.  If it isn't, a search is made for a crossing
784 *   solution by iterating on the unknown celestial coordinate starting at the
785 *   upper limit of the solution interval and decrementing by the specified
786 *   step size.  A crossing is indicated if the trial value of the pixel
787 *   coordinate steps through the value specified.  If a crossing interval is
788 *   found then the solution is determined by a modified form of "regula falsi"
789 *   division of the crossing interval.  If no crossing interval was found
790 *   within the specified solution interval then a search is made for a
791 *   "non-crossing" solution as may arise from a point of tangency.  The
792 *   process is complicated by having to make allowance for the discontinuities
793 *   that occur in all map projections.
794 *
795 *   Once one solution has been determined others may be found by subsequent
796 *   invokations of wcsmix() with suitably restricted solution intervals.
797 *
798 *   Note the circumstance that arises when the solution point lies at a native
799 *   pole of a projection in which the pole is represented as a finite curve,
800 *   for example the zenithals and conics.  In such cases two or more valid
801 *   solutions may exist but wcsmix() only ever returns one.
802 *
803 *   Because of its generality wcsmix() is very compute-intensive.  For
804 *   compute-limited applications more efficient special-case solvers could be
805 *   written for simple projections, for example non-oblique cylindrical
806 *   projections.
807 *
808 *
809 * wcssptr() - Spectral axis translation
810 * -------------------------------------
811 * wcssptr() translates the spectral axis in a wcsprm struct.  For example, a
812 * 'FREQ' axis may be translated into 'ZOPT-F2W' and vice versa.
813 *
814 * Given and returned:
815 *   wcs       struct wcsprm*
816 *                       Coordinate transformation parameters.
817 *
818 *   i         int*      Index of the spectral axis (0-relative).  If given < 0
819 *                       it will be set to the first spectral axis identified
820 *                       from the ctype[] keyvalues in the wcsprm struct.
821 *
822 *   ctype     char[9]   Desired spectral CTYPEia.  Wildcarding may be used as
823 *                       for the ctypeS2 argument to spctrn() as described in
824 *                       the prologue of spc.h, i.e. if the final three
825 *                       characters are specified as "???", or if just the
826 *                       eighth character is specified as '?', the correct
827 *                       algorithm code will be substituted and returned.
828 *
829 * Function return value:
830 *             int       Status return value:
831 *                         0: Success.
832 *                         1: Null wcsprm pointer passed.
833 *                         2: Memory allocation failed.
834 *                         3: Linear transformation matrix is singular.
835 *                         4: Inconsistent or unrecognized coordinate axis
836 *                            types.
837 *                         5: Invalid parameter value.
838 *                         6: Invalid coordinate transformation parameters.
839 *                         7: Ill-conditioned coordinate transformation
840 *                            parameters.
841 *                        12: Invalid subimage specification (no spectral
842 *                            axis).
843 *
844 *                       For returns > 1, a detailed error message is set in
845 *                       wcsprm::err if enabled, see wcserr_enable().
846 *
847 *
848 * wcslib_version() - WCSLIB version number
849 * ----------------------------------------
850 * wcslib_version() returns the WCSLIB version number.
851 *
852 * The major version number changes when the ABI changes or when the license
853 * conditions change.  ABI changes typically result from a change to the
854 * contents of one of the structs.  The major version number is used to
855 * distinguish between incompatible versions of the sharable library.
856 *
857 * The minor version number changes with new functionality or bug fixes that do
858 * not involve a change in the ABI.
859 *
860 * The auxiliary version number (which is often absent) signals changes to the
861 * documentation, test suite, build procedures, or any other change that does
862 * not affect the compiled library.
863 *
864 * Returned:
865 *   vers[3]   int[3]    The broken-down version number:
866 *                         0: Major version number.
867 *                         1: Minor version number.
868 *                         2: Auxiliary version number (zero if absent).
869 *                       May be given as a null pointer if not required.
870 *
871 * Function return value:
872 *             char*     A null-terminated, statically allocated string
873 *                       containing the version number in the usual form, i.e.
874 *                       "<major>.<minor>.<auxiliary>".
875 *
876 *
877 * wcsprm struct - Coordinate transformation parameters
878 * ----------------------------------------------------
879 * The wcsprm struct contains information required to transform world
880 * coordinates.  It consists of certain members that must be set by the user
881 * ("given") and others that are set by the WCSLIB routines ("returned").
882 * While the addresses of the arrays themselves may be set by wcsinit() if it
883 * (optionally) allocates memory, their contents must be set by the user.
884 *
885 * Some parameters that are given are not actually required for transforming
886 * coordinates.  These are described as "auxiliary"; the struct simply provides
887 * a place to store them, though they may be used by wcshdo() in constructing a
888 * FITS header from a wcsprm struct.  Some of the returned values are supplied
889 * for informational purposes and others are for internal use only as
890 * indicated.
891 *
892 * In practice, it is expected that a WCS parser would scan the FITS header to
893 * determine the number of coordinate axes.  It would then use wcsinit() to
894 * allocate memory for arrays in the wcsprm struct and set default values.
895 * Then as it reread the header and identified each WCS keyrecord it would load
896 * the value into the relevant wcsprm array element.  This is essentially what
897 * wcspih() does - refer to the prologue of wcshdr.h.  As the final step,
898 * wcsset() is invoked, either directly or indirectly, to set the derived
899 * members of the wcsprm struct.  wcsset() strips off trailing blanks in all
900 * string members and null-fills the character array.
901 *
902 *   int flag
903 *     (Given and returned) This flag must be set to zero whenever any of the
904 *     following wcsprm struct members are set or changed:
905 *
906 *       - wcsprm::naxis (q.v., not normally set by the user),
907 *       - wcsprm::crpix,
908 *       - wcsprm::pc,
909 *       - wcsprm::cdelt,
910 *       - wcsprm::crval,
911 *       - wcsprm::cunit,
912 *       - wcsprm::ctype,
913 *       - wcsprm::lonpole,
914 *       - wcsprm::latpole,
915 *       - wcsprm::restfrq,
916 *       - wcsprm::restwav,
917 *       - wcsprm::npv,
918 *       - wcsprm::pv,
919 *       - wcsprm::nps,
920 *       - wcsprm::ps,
921 *       - wcsprm::cd,
922 *       - wcsprm::crota,
923 *       - wcsprm::altlin,
924 *       - wcsprm::ntab,
925 *       - wcsprm::nwtb,
926 *       - wcsprm::tab,
927 *       - wcsprm::wtb.
928 *
929 *     This signals the initialization routine, wcsset(), to recompute the
930 *     returned members of the celprm struct.  celset() will reset flag to
931 *     indicate that this has been done.
932 *
933 *     PLEASE NOTE: flag should be set to -1 when wcsinit() is called for the
934 *     first time for a particular wcsprm struct in order to initialize memory
935 *     management.  It must ONLY be used on the first initialization otherwise
936 *     memory leaks may result.
937 *
938 *   int naxis
939 *     (Given or returned) Number of pixel and world coordinate elements.
940 *
941 *     If wcsinit() is used to initialize the linprm struct (as would normally
942 *     be the case) then it will set naxis from the value passed to it as a
943 *     function argument.  The user should not subsequently modify it.
944 *
945 *   double *crpix
946 *     (Given) Address of the first element of an array of double containing
947 *     the coordinate reference pixel, CRPIXja.
948 *
949 *   double *pc
950 *     (Given) Address of the first element of the PCi_ja (pixel coordinate)
951 *     transformation matrix.  The expected order is
952 *
953 =       struct wcsprm wcs;
954 =       wcs.pc = {PC1_1, PC1_2, PC2_1, PC2_2};
955 *
956 *     This may be constructed conveniently from a 2-D array via
957 *
958 =       double m[2][2] = {{PC1_1, PC1_2},
959 =                         {PC2_1, PC2_2}};
960 *
961 *     which is equivalent to
962 *
963 =       double m[2][2];
964 =       m[0][0] = PC1_1;
965 =       m[0][1] = PC1_2;
966 =       m[1][0] = PC2_1;
967 =       m[1][1] = PC2_2;
968 *
969 *     The storage order for this 2-D array is the same as for the 1-D array,
970 *     whence
971 *
972 =       wcs.pc = *m;
973 *
974 *     would be legitimate.
975 *
976 *   double *cdelt
977 *     (Given) Address of the first element of an array of double containing
978 *     the coordinate increments, CDELTia.
979 *
980 *   double *crval
981 *     (Given) Address of the first element of an array of double containing
982 *     the coordinate reference values, CRVALia.
983 *
984 *   char (*cunit)[72]
985 *     (Given) Address of the first element of an array of char[72] containing
986 *     the CUNITia keyvalues which define the units of measurement of the
987 *     CRVALia, CDELTia, and CDi_ja keywords.
988 *
989 *     As CUNITia is an optional header keyword, cunit[][72] may be left blank
990 *     but otherwise is expected to contain a standard units specification as
991 *     defined by WCS Paper I.  Utility function wcsutrn(), described in
992 *     wcsunits.h, is available to translate commonly used non-standard units
993 *     specifications but this must be done as a separate step before invoking
994 *     wcsset().
995 *
996 *     For celestial axes, if cunit[][72] is not blank, wcsset() uses
997 *     wcsunits() to parse it and scale cdelt[], crval[], and cd[][*] to
998 *     degrees.  It then resets cunit[][72] to "deg".
999 *
1000 *     For spectral axes, if cunit[][72] is not blank, wcsset() uses wcsunits()
1001 *     to parse it and scale cdelt[], crval[], and cd[][*] to SI units.  It
1002 *     then resets cunit[][72] accordingly.
1003 *
1004 *     wcsset() ignores cunit[][72] for other coordinate types; cunit[][72] may
1005 *     be used to label coordinate values.
1006 *
1007 *     These variables accomodate the longest allowed string-valued FITS
1008 *     keyword, being limited to 68 characters, plus the null-terminating
1009 *     character.
1010 *
1011 *   char (*ctype)[72]
1012 *     (Given) Address of the first element of an array of char[72] containing
1013 *     the coordinate axis types, CTYPEia.
1014 *
1015 *     The ctype[][72] keyword values must be in upper case and there must be
1016 *     zero or one pair of matched celestial axis types, and zero or one
1017 *     spectral axis.  The ctype[][72] strings should be padded with blanks on
1018 *     the right and null-terminated so that they are at least eight characters
1019 *     in length.
1020 *
1021 *     These variables accomodate the longest allowed string-valued FITS
1022 *     keyword, being limited to 68 characters, plus the null-terminating
1023 *     character.
1024 *
1025 *   double lonpole
1026 *     (Given and returned) The native longitude of the celestial pole, phi_p,
1027 *     given by LONPOLEa [deg] or by PVi_2a [deg] attached to the longitude
1028 *     axis which takes precedence if defined, and ...
1029 *   double latpole
1030 *     (Given and returned) ... the native latitude of the celestial pole,
1031 *     theta_p, given by LATPOLEa [deg] or by PVi_3a [deg] attached to the
1032 *     longitude axis which takes precedence if defined.
1033 *
1034 *     lonpole and latpole may be left to default to values set by wcsinit()
1035 *     (see celprm::ref), but in any case they will be reset by wcsset() to
1036 *     the values actually used.  Note therefore that if the wcsprm struct is
1037 *     reused without resetting them, whether directly or via wcsinit(), they
1038 *     will no longer have their default values.
1039 *
1040 *   double restfrq
1041 *     (Given) The rest frequency [Hz], and/or ...
1042 *   double restwav
1043 *     (Given) ... the rest wavelength in vacuo [m], only one of which need be
1044 *     given, the other should be set to zero.
1045 *
1046 *   int npv
1047 *     (Given) The number of entries in the wcsprm::pv[] array.
1048 *
1049 *   int npvmax
1050 *     (Given or returned) The length of the wcsprm::pv[] array.
1051 *
1052 *     npvmax will be set by wcsinit() if it allocates memory for wcsprm::pv[],
1053 *     otherwise it must be set by the user.  See also wcsnpv().
1054 *
1055 *   struct pvcard *pv
1056 *     (Given) Address of the first element of an array of length npvmax of
1057 *     pvcard structs.
1058 *
1059 *     As a FITS header parser encounters each PVi_ma keyword it should load it
1060 *     into a pvcard struct in the array and increment npv.  wcsset()
1061 *     interprets these as required.
1062 *
1063 *     Note that, if they were not given, wcsset() resets the entries for
1064 *     PVi_1a, PVi_2a, PVi_3a, and PVi_4a for longitude axis i to match
1065 *     phi_0 and theta_0 (the native longitude and latitude of the reference
1066 *     point), LONPOLEa and LATPOLEa respectively.
1067 *
1068 *   int nps
1069 *     (Given) The number of entries in the wcsprm::ps[] array.
1070 *
1071 *   int npsmax
1072 *     (Given or returned) The length of the wcsprm::ps[] array.
1073 *
1074 *     npsmax will be set by wcsinit() if it allocates memory for wcsprm::ps[],
1075 *     otherwise it must be set by the user.  See also wcsnps().
1076 *
1077 *   struct pscard *ps
1078 *     (Given) Address of the first element of an array of length npsmax of
1079 *     pscard structs.
1080 *
1081 *     As a FITS header parser encounters each PSi_ma keyword it should load it
1082 *     into a pscard struct in the array and increment nps.  wcsset()
1083 *     interprets these as required (currently no PSi_ma keyvalues are
1084 *     recognized).
1085 *
1086 *   double *cd
1087 *     (Given) For historical compatibility, the wcsprm struct supports two
1088 *     alternate specifications of the linear transformation matrix, those
1089 *     associated with the CDi_ja keywords, and ...
1090 *   double *crota
1091 *     (Given) ... those associated with the CROTAi keywords.  Although these
1092 *     may not formally co-exist with PCi_ja, the approach taken here is simply
1093 *     to ignore them if given in conjunction with PCi_ja.
1094 *
1095 *   int altlin
1096 *     (Given) altlin is a bit flag that denotes which of the PCi_ja, CDi_ja
1097 *     and CROTAi keywords are present in the header:
1098 *
1099 *     - Bit 0: PCi_ja is present.
1100 *
1101 *     - Bit 1: CDi_ja is present.
1102 *
1103 *       Matrix elements in the IRAF convention are
1104 *       equivalent to the product CDi_ja = CDELTia * PCi_ja, but the
1105 *       defaults differ from that of the PCi_ja matrix.  If one or more
1106 *       CDi_ja keywords are present then all unspecified CDi_ja default to
1107 *       zero.  If no CDi_ja (or CROTAi) keywords are present, then the
1108 *       header is assumed to be in PCi_ja form whether or not any PCi_ja
1109 *       keywords are present since this results in an interpretation of
1110 *       CDELTia consistent with the original FITS specification.
1111 *
1112 *       While CDi_ja may not formally co-exist with PCi_ja, it may co-exist
1113 *       with CDELTia and CROTAi which are to be ignored.
1114 *
1115 *     - Bit 2: CROTAi is present.
1116 *
1117 *       In the AIPS convention, CROTAi may only be
1118 *       associated with the latitude axis of a celestial axis pair.  It
1119 *       specifies a rotation in the image plane that is applied AFTER the
1120 *       CDELTia; any other CROTAi keywords are ignored.
1121 *
1122 *       CROTAi may not formally co-exist with PCi_ja.
1123 *
1124 *       CROTAi and CDELTia may formally co-exist with CDi_ja but if so are to
1125 *       be ignored.
1126 *
1127 *     CDi_ja and CROTAi keywords, if found, are to be stored in the
1128 *     wcsprm::cd and wcsprm::crota arrays which are dimensioned similarly to
1129 *     wcsprm::pc and wcsprm::cdelt.  FITS
1130 *     header parsers should use the following procedure:
1131 *
1132 *     - Whenever a PCi_ja  keyword is encountered: altlin |= 1;
1133 *
1134 *     - Whenever a CDi_ja  keyword is encountered: altlin |= 2;
1135 *
1136 *     - Whenever a CROTAi keyword is encountered: altlin |= 4;
1137 *
1138 *     If none of these bits are set the PCi_ja representation results, i.e.
1139 *     wcsprm::pc and wcsprm::cdelt will be used as given.
1140 *
1141 *     These alternate specifications of the linear transformation matrix are
1142 *     translated immediately to PCi_ja by wcsset() and are invisible to the
1143 *     lower-level WCSLIB routines.  In particular, wcsset() resets
1144 *     wcsprm::cdelt to unity if CDi_ja is present (and no PCi_ja).
1145 *
1146 *     If CROTAi are present but none is associated with the latitude axis
1147 *     (and no PCi_ja or CDi_ja), then wcsset() reverts to a unity PCi_ja
1148 *     matrix.
1149 *
1150 *   int velref
1151 *     (Given) AIPS velocity code VELREF, refer to spcaips().
1152 *
1153 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1154 *     wcsprm::velref is changed.
1155 *
1156 *   char alt[4]
1157 *     (Given, auxiliary) Character code for alternate coordinate descriptions
1158 *     (i.e. the 'a' in keyword names such as CTYPEia).  This is blank for the
1159 *     primary coordinate description, or one of the 26 upper-case letters,
1160 *     A-Z.
1161 *
1162 *     An array of four characters is provided for alignment purposes, only the
1163 *     first is used.
1164 *
1165 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1166 *     wcsprm::alt is changed.
1167 *
1168 *   int colnum
1169 *     (Given, auxiliary) Where the coordinate representation is associated
1170 *     with an image-array column in a FITS binary table, this variable may be
1171 *     used to record the relevant column number.
1172 *
1173 *     It should be set to zero for an image header or pixel list.
1174 *
1175 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1176 *     wcsprm::colnum is changed.
1177 *
1178 *   int *colax
1179 *     (Given, auxiliary) Address of the first element of an array of int
1180 *     recording the column numbers for each axis in a pixel list.
1181 *
1182 *     The array elements should be set to zero for an image header or image
1183 *     array in a binary table.
1184 *
1185 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1186 *     wcsprm::colax is changed.
1187 *
1188 *   char (*cname)[72]
1189 *     (Given, auxiliary) The address of the first element of an array of
1190 *     char[72] containing the coordinate axis names, CNAMEia.
1191 *
1192 *     These variables accomodate the longest allowed string-valued FITS
1193 *     keyword, being limited to 68 characters, plus the null-terminating
1194 *     character.
1195 *
1196 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1197 *     wcsprm::cname is changed.
1198 *
1199 *   double *crder
1200 *     (Given, auxiliary) Address of the first element of an array of double
1201 *     recording the random error in the coordinate value, CRDERia.
1202 *
1203 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1204 *     wcsprm::crder is changed.
1205 *
1206 *   double *csyer
1207 *     (Given, auxiliary) Address of the first element of an array of double
1208 *     recording the systematic error in the coordinate value, CSYERia.
1209 *
1210 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1211 *     wcsprm::csyer is changed.
1212 *
1213 *   double *czphs
1214 *     (Given, auxiliary) Address of the first element of an array of double
1215 *     recording the time at the zero point of a phase axis, CZPHSia.
1216 *
1217 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1218 *     wcsprm::czphs is changed.
1219 *
1220 *   double *cperi
1221 *     (Given, auxiliary) Address of the first element of an array of double
1222 *     recording the period of a phase axis, CPERIia.
1223 *
1224 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1225 *     wcsprm::cperi is changed.
1226 *
1227 *   char wcsname[72]
1228 *     (Given, auxiliary) The name given to the coordinate representation,
1229 *     WCSNAMEa.  This variable accomodates the longest allowed string-valued
1230 *     FITS keyword, being limited to 68 characters, plus the null-terminating
1231 *     character.
1232 *
1233 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1234 *     wcsprm::wcsname is changed.
1235 *
1236 *   char timesys[72]
1237 *     (Given, auxiliary) TIMESYS keyvalue, being the time scale (UTC, TAI,
1238 *     etc.) in which all other time-related auxiliary header values are
1239 *     recorded.  Also defines the time scale for an image axis with CTYPEia
1240 *     set to 'TIME'.
1241 *
1242 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1243 *     wcsprm::timesys is changed.
1244 *
1245 *   char trefpos[72]
1246 *     (Given, auxiliary) TREFPOS keyvalue, being the location in space where
1247 *     the recorded time is valid.
1248 *
1249 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1250 *     wcsprm::trefpos is changed.
1251 *
1252 *   char trefdir[72]
1253 *     (Given, auxiliary) TREFDIR keyvalue, being the reference direction used
1254 *     in calculating a pathlength delay.
1255 *
1256 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1257 *     wcsprm::trefdir is changed.
1258 *
1259 *   char plephem[72]
1260 *     (Given, auxiliary) PLEPHEM keyvalue, being the Solar System ephemeris
1261 *     used for calculating a pathlength delay.
1262 *
1263 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1264 *     wcsprm::plephem is changed.
1265 *
1266 *   char timeunit[72]
1267 *     (Given, auxiliary) TIMEUNIT keyvalue, being the time units in which
1268 *     the following header values are expressed: TSTART, TSTOP, TIMEOFFS,
1269 *     TIMSYER, TIMRDER, TIMEDEL.  It also provides the default value for
1270 *     CUNITia for time axes.
1271 *
1272 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1273 *     wcsprm::timeunit is changed.
1274 *
1275 *   char dateref[72]
1276 *     (Given, auxiliary) DATEREF keyvalue, being the date of a reference epoch
1277 *     relative to which other time measurements refer.
1278 *
1279 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1280 *     wcsprm::dateref is changed.
1281 *
1282 *   double mjdref[2]
1283 *     (Given, auxiliary) MJDREF keyvalue, equivalent to DATEREF expressed as
1284 *     a Modified Julian Date (MJD = JD - 2400000.5).  The value is given as
1285 *     the sum of the two-element vector, allowing increased precision.
1286 *
1287 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1288 *     wcsprm::mjdref is changed.
1289 *
1290 *   double timeoffs
1291 *     (Given, auxiliary) TIMEOFFS keyvalue, being a time offset, which may be
1292 *     used, for example, to provide a uniform clock correction for times
1293 *     referenced to DATEREF.
1294 *
1295 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1296 *     wcsprm::timeoffs is changed.
1297 *
1298 *   char dateobs[72]
1299 *     (Given, auxiliary) DATE-OBS keyvalue, being the date at the start of the
1300 *     observation unless otherwise explained in the DATE-OBS keycomment, in
1301 *     ISO format, yyyy-mm-ddThh:mm:ss.
1302 *
1303 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1304 *     wcsprm::dateobs is changed.
1305 *
1306 *   char datebeg[72]
1307 *     (Given, auxiliary) DATE-BEG keyvalue, being the date at the start of the
1308 *     observation in ISO format, yyyy-mm-ddThh:mm:ss.
1309 *
1310 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1311 *     wcsprm::datebeg is changed.
1312 *
1313 *   char dateavg[72]
1314 *     (Given, auxiliary) DATE-AVG keyvalue, being the date at a representative
1315 *     mid-point of the observation in ISO format, yyyy-mm-ddThh:mm:ss.
1316 *
1317 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1318 *     wcsprm::dateavg is changed.
1319 *
1320 *   char dateend[72]
1321 *     (Given, auxiliary) DATE-END keyvalue, baing the date at the end of the
1322 *     observation in ISO format, yyyy-mm-ddThh:mm:ss.
1323 *
1324 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1325 *     wcsprm::dateend is changed.
1326 *
1327 *   double mjdobs
1328 *     (Given, auxiliary) MJD-OBS keyvalue, equivalent to DATE-OBS expressed
1329 *     as a Modified Julian Date (MJD = JD - 2400000.5).
1330 *
1331 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1332 *     wcsprm::mjdobs is changed.
1333 *
1334 *   double mjdbeg
1335 *     (Given, auxiliary) MJD-BEG keyvalue, equivalent to DATE-BEG expressed
1336 *     as a Modified Julian Date (MJD = JD - 2400000.5).
1337 *
1338 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1339 *     wcsprm::mjdbeg is changed.
1340 *
1341 *   double mjdavg
1342 *     (Given, auxiliary) MJD-AVG keyvalue, equivalent to DATE-AVG expressed
1343 *     as a Modified Julian Date (MJD = JD - 2400000.5).
1344 *
1345 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1346 *     wcsprm::mjdavg is changed.
1347 *
1348 *   double mjdend
1349 *     (Given, auxiliary) MJD-END keyvalue, equivalent to DATE-END expressed
1350 *     as a Modified Julian Date (MJD = JD - 2400000.5).
1351 *
1352 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1353 *     wcsprm::mjdend is changed.
1354 *
1355 *   double jepoch
1356 *     (Given, auxiliary) JEPOCH keyvalue, equivalent to DATE-OBS expressed
1357 *     as a Julian epoch.
1358 *
1359 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1360 *     wcsprm::jepoch is changed.
1361 *
1362 *   double bepoch
1363 *     (Given, auxiliary) BEPOCH keyvalue, equivalent to DATE-OBS expressed
1364 *     as a Besselian epoch
1365 *
1366 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1367 *     wcsprm::bepoch is changed.
1368 *
1369 *   double tstart
1370 *     (Given, auxiliary) TSTART keyvalue, equivalent to DATE-BEG expressed
1371 *     as a time in units of TIMEUNIT relative to DATEREF+TIMEOFFS.
1372 *
1373 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1374 *     wcsprm::tstart is changed.
1375 *
1376 *   double tstop
1377 *     (Given, auxiliary) TSTOP keyvalue, equivalent to DATE-END expressed
1378 *     as a time in units of TIMEUNIT relative to DATEREF+TIMEOFFS.
1379 *
1380 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1381 *     wcsprm::tstop is changed.
1382 *
1383 *   double xposure
1384 *     (Given, auxiliary) XPOSURE keyvalue, being the effective exposure time
1385 *     in units of TIMEUNIT.
1386 *
1387 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1388 *     wcsprm::xposure is changed.
1389 *
1390 *   double telapse
1391 *     (Given, auxiliary) TELAPSE keyvalue, equivalent to the elapsed time
1392 *     between DATE-BEG and DATE-END, in units of TIMEUNIT.
1393 *
1394 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1395 *     wcsprm::telapse is changed.
1396 *
1397 *   double timsyer
1398 *     (Given, auxiliary) TIMSYER keyvalue, being the absolute error of the
1399 *     time values, in units of TIMEUNIT.
1400 *
1401 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1402 *     wcsprm::timsyer is changed.
1403 *
1404 *   double timrder
1405 *     (Given, auxiliary) TIMRDER keyvalue, being the accuracy of time stamps
1406 *     relative to each other, in units of TIMEUNIT.
1407 *
1408 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1409 *     wcsprm::timrder is changed.
1410 *
1411 *   double timedel
1412 *     (Given, auxiliary) TIMEDEL keyvalue, being the resolution of the time
1413 *     stamps.
1414 *
1415 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1416 *     wcsprm::timedel is changed.
1417 *
1418 *   double timepixr
1419 *     (Given, auxiliary) TIMEPIXR keyvalue, being the relative position of the
1420 *     time stamps in binned time intervals, a value between 0.0 and 1.0.
1421 *
1422 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1423 *     wcsprm::timepixr is changed.
1424 *
1425 *   double obsgeo[6]
1426 *     (Given, auxiliary) Location of the observer in a standard terrestrial
1427 *     reference frame.  The first three give ITRS Cartesian coordinates
1428 *     OBSGEO-X [m],   OBSGEO-Y [m],   OBSGEO-Z [m], and the second three give
1429 *     OBSGEO-L [deg], OBSGEO-B [deg], OBSGEO-H [m], which are related through
1430 *     a standard transformation.
1431 *
1432 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1433 *     wcsprm::obsgeo is changed.
1434 *
1435 *   char obsorbit[72]
1436 *     (Given, auxiliary) OBSORBIT keyvalue, being the URI, URL, or name of an
1437 *     orbit ephemeris file giving spacecraft coordinates relating to TREFPOS.
1438 *
1439 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1440 *     wcsprm::obsorbit is changed.
1441 *
1442 *   char radesys[72]
1443 *     (Given, auxiliary) The equatorial or ecliptic coordinate system type,
1444 *     RADESYSa.
1445 *
1446 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1447 *     wcsprm::radesys is changed.
1448 *
1449 *   double equinox
1450 *     (Given, auxiliary) The equinox associated with dynamical equatorial or
1451 *     ecliptic coordinate systems, EQUINOXa (or EPOCH in older headers).  Not
1452 *     applicable to ICRS equatorial or ecliptic coordinates.
1453 *
1454 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1455 *     wcsprm::equinox is changed.
1456 *
1457 *   char specsys[72]
1458 *     (Given, auxiliary) Spectral reference frame (standard of rest),
1459 *     SPECSYSa.
1460 *
1461 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1462 *     wcsprm::specsys is changed.
1463 *
1464 *   char ssysobs[72]
1465 *     (Given, auxiliary) The spectral reference frame in which there is no
1466 *     differential variation in the spectral coordinate across the
1467 *     field-of-view, SSYSOBSa.
1468 *
1469 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1470 *     wcsprm::ssysobs is changed.
1471 *
1472 *   double velosys
1473 *     (Given, auxiliary) The relative radial velocity [m/s] between the
1474 *     observer and the selected standard of rest in the direction of the
1475 *     celestial reference coordinate, VELOSYSa.
1476 *
1477 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1478 *     wcsprm::velosys is changed.
1479 *
1480 *   double zsource
1481 *     (Given, auxiliary) The redshift, ZSOURCEa, of the source.
1482 *
1483 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1484 *     wcsprm::zsource is changed.
1485 *
1486 *   char ssyssrc[72]
1487 *     (Given, auxiliary) The spectral reference frame (standard of rest),
1488 *     SSYSSRCa, in which wcsprm::zsource was measured.
1489 *
1490 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1491 *     wcsprm::ssyssrc is changed.
1492 *
1493 *   double velangl
1494 *     (Given, auxiliary) The angle [deg] that should be used to decompose an
1495 *     observed velocity into radial and transverse components.
1496 *
1497 *     It is not necessary to reset the wcsprm struct (via wcsset()) when
1498 *     wcsprm::velangl is changed.
1499 *
1500 *   struct auxprm *aux
1501 *     (Given, auxiliary) This struct holds auxiliary coordinate system
1502 *     information of a specialist nature.  While these parameters may be
1503 *     widely recognized within particular fields of astronomy, they differ
1504 *     from the above auxiliary parameters in not being defined by any of the
1505 *     FITS WCS standards.  Collecting them together in a separate struct that
1506 *     is allocated only when required helps to control bloat in the size of
1507 *     the wcsprm struct.
1508 *
1509 *   int ntab
1510 *     (Given) See wcsprm::tab.
1511 *
1512 *   int nwtb
1513 *     (Given) See wcsprm::wtb.
1514 *
1515 *   struct tabprm *tab
1516 *     (Given) Address of the first element of an array of ntab tabprm structs
1517 *     for which memory has been allocated.  These are used to store tabular
1518 *     transformation parameters.
1519 *
1520 *     Although technically wcsprm::ntab and tab are "given", they will
1521 *     normally be set by invoking wcstab(), whether directly or indirectly.
1522 *
1523 *     The tabprm structs contain some members that must be supplied and others
1524 *     that are derived.  The information to be supplied comes primarily from
1525 *     arrays stored in one or more FITS binary table extensions.  These
1526 *     arrays, referred to here as "wcstab arrays", are themselves located by
1527 *     parameters stored in the FITS image header.
1528 *
1529 *   struct wtbarr *wtb
1530 *     (Given) Address of the first element of an array of nwtb wtbarr structs
1531 *     for which memory has been allocated.  These are used in extracting
1532 *     wcstab arrays from a FITS binary table.
1533 *
1534 *     Although technically wcsprm::nwtb and wtb are "given", they will
1535 *     normally be set by invoking wcstab(), whether directly or indirectly.
1536 *
1537 *   char lngtyp[8]
1538 *     (Returned) Four-character WCS celestial longitude and ...
1539 *   char lattyp[8]
1540 *     (Returned) ... latitude axis types. e.g. "RA", "DEC", "GLON", "GLAT",
1541 *     etc. extracted from 'RA--', 'DEC-', 'GLON', 'GLAT', etc. in the first
1542 *     four characters of CTYPEia but with trailing dashes removed.  (Declared
1543 *     as char[8] for alignment reasons.)
1544 *
1545 *   int lng
1546 *     (Returned) Index for the longitude coordinate, and ...
1547 *   int lat
1548 *     (Returned) ... index for the latitude coordinate, and ...
1549 *   int spec
1550 *     (Returned) ... index for the spectral coordinate in the imgcrd[][] and
1551 *     world[][] arrays in the API of wcsp2s(), wcss2p() and wcsmix().
1552 *
1553 *     These may also serve as indices into the pixcrd[][] array provided that
1554 *     the PCi_ja matrix does not transpose axes.
1555 *
1556 *   int cubeface
1557 *     (Returned) Index into the pixcrd[][] array for the CUBEFACE axis.  This
1558 *     is used for quadcube projections where the cube faces are stored on a
1559 *     separate axis (see wcs.h).
1560 *
1561 *   int *types
1562 *     (Returned) Address of the first element of an array of int containing a
1563 *     four-digit type code for each axis.
1564 *
1565 *     - First digit (i.e. 1000s):
1566 *       - 0: Non-specific coordinate type.
1567 *       - 1: Stokes coordinate.
1568 *       - 2: Celestial coordinate (including CUBEFACE).
1569 *       - 3: Spectral coordinate.
1570 *
1571 *     - Second digit (i.e. 100s):
1572 *       - 0: Linear axis.
1573 *       - 1: Quantized axis (STOKES, CUBEFACE).
1574 *       - 2: Non-linear celestial axis.
1575 *       - 3: Non-linear spectral axis.
1576 *       - 4: Logarithmic axis.
1577 *       - 5: Tabular axis.
1578 *
1579 *     - Third digit (i.e. 10s):
1580 *       - 0: Group number, e.g. lookup table number, being an index into the
1581 *            tabprm array (see above).
1582 *
1583 *     - The fourth digit is used as a qualifier depending on the axis type.
1584 *
1585 *       - For celestial axes:
1586 *         - 0: Longitude coordinate.
1587 *         - 1: Latitude coordinate.
1588 *         - 2: CUBEFACE number.
1589 *
1590 *       - For lookup tables: the axis number in a multidimensional table.
1591 *
1592 *     CTYPEia in "4-3" form with unrecognized algorithm code will have its
1593 *     type set to -1 and generate an error.
1594 *
1595 *   struct linprm lin
1596 *     (Returned) Linear transformation parameters (usage is described in the
1597 *     prologue to lin.h).
1598 *
1599 *   struct celprm cel
1600 *     (Returned) Celestial transformation parameters (usage is described in
1601 *     the prologue to cel.h).
1602 *
1603 *   struct spcprm spc
1604 *     (Returned) Spectral transformation parameters (usage is described in the
1605 *     prologue to spc.h).
1606 *
1607 *   struct wcserr *err
1608 *     (Returned) If enabled, when an error status is returned, this struct
1609 *     contains detailed information about the error, see wcserr_enable().
1610 *
1611 *   int m_flag
1612 *     (For internal use only.)
1613 *   int m_naxis
1614 *     (For internal use only.)
1615 *   double *m_crpix
1616 *     (For internal use only.)
1617 *   double *m_pc
1618 *     (For internal use only.)
1619 *   double *m_cdelt
1620 *     (For internal use only.)
1621 *   double *m_crval
1622 *     (For internal use only.)
1623 *   char (*m_cunit)[72]
1624 *     (For internal use only.)
1625 *   char (*m_ctype)[72]
1626 *     (For internal use only.)
1627 *   struct pvcard *m_pv
1628 *     (For internal use only.)
1629 *   struct pscard *m_ps
1630 *     (For internal use only.)
1631 *   double *m_cd
1632 *     (For internal use only.)
1633 *   double *m_crota
1634 *     (For internal use only.)
1635 *   int *m_colax
1636 *     (For internal use only.)
1637 *   char (*m_cname)[72]
1638 *     (For internal use only.)
1639 *   double *m_crder
1640 *     (For internal use only.)
1641 *   double *m_csyer
1642 *     (For internal use only.)
1643 *   double *m_czphs
1644 *     (For internal use only.)
1645 *   double *m_cperi
1646 *     (For internal use only.)
1647 *   struct tabprm *m_tab
1648 *     (For internal use only.)
1649 *   struct wtbarr *m_wtb
1650 *     (For internal use only.)
1651 *
1652 *
1653 * pvcard struct - Store for PVi_ma keyrecords
1654 * -------------------------------------------
1655 * The pvcard struct is used to pass the parsed contents of PVi_ma keyrecords
1656 * to wcsset() via the wcsprm struct.
1657 *
1658 * All members of this struct are to be set by the user.
1659 *
1660 *   int i
1661 *     (Given) Axis number (1-relative), as in the FITS PVi_ma keyword.  If
1662 *     i == 0, wcsset() will replace it with the latitude axis number.
1663 *
1664 *   int m
1665 *     (Given) Parameter number (non-negative), as in the FITS PVi_ma keyword.
1666 *
1667 *   double value
1668 *     (Given) Parameter value.
1669 *
1670 *
1671 * pscard struct - Store for PSi_ma keyrecords
1672 * -------------------------------------------
1673 * The pscard struct is used to pass the parsed contents of PSi_ma keyrecords
1674 * to wcsset() via the wcsprm struct.
1675 *
1676 * All members of this struct are to be set by the user.
1677 *
1678 *   int i
1679 *     (Given) Axis number (1-relative), as in the FITS PSi_ma keyword.
1680 *
1681 *   int m
1682 *     (Given) Parameter number (non-negative), as in the FITS PSi_ma keyword.
1683 *
1684 *   char value[72]
1685 *     (Given) Parameter value.
1686 *
1687 *
1688 * auxprm struct - Additional auxiliary parameters
1689 * -----------------------------------------------
1690 * The auxprm struct holds auxiliary coordinate system information of a
1691 * specialist nature.  It is anticipated that this struct will expand in future
1692 * to accomodate additional parameters.
1693 *
1694 * All members of this struct are to be set by the user.
1695 *
1696 *   double rsun_ref
1697 *     (Given, auxiliary) Reference radius of the Sun used in coordinate
1698 *     calculations (m).
1699 *
1700 *   double dsun_obs
1701 *     (Given, auxiliary) Distance between the centre of the Sun and the
1702 *     observer (m).
1703 *
1704 *   double crln_obs
1705 *     (Given, auxiliary) Carrington heliographic longitude of the observer
1706 *     (deg).
1707 *
1708 *   double hgln_obs
1709 *     (Given, auxiliary) Stonyhurst heliographic longitude of the observer
1710 *     (deg).
1711 *
1712 *   double hglt_obs
1713 *     (Given, auxiliary) Heliographic latitude (Carrington or Stonyhurst) of
1714 *     the observer (deg).
1715 *
1716 *
1717 * Global variable: const char *wcs_errmsg[] - Status return messages
1718 * ------------------------------------------------------------------
1719 * Error messages to match the status value returned from each function.
1720 *
1721 *===========================================================================*/
1722 
1723 #ifndef WCSLIB_WCS
1724 #define WCSLIB_WCS
1725 
1726 #define wcsset wcs_set
1727 
1728 #include "lin.h"
1729 #include "cel.h"
1730 #include "spc.h"
1731 
1732 #ifdef __cplusplus
1733 extern "C" {
1734 #define wtbarr wtbarr_s		// See prologue of wtbarr.h.
1735 #endif
1736 
1737 #define WCSSUB_LONGITUDE 0x1001
1738 #define WCSSUB_LATITUDE  0x1002
1739 #define WCSSUB_CUBEFACE  0x1004
1740 #define WCSSUB_CELESTIAL 0x1007
1741 #define WCSSUB_SPECTRAL  0x1008
1742 #define WCSSUB_STOKES    0x1010
1743 
1744 
1745 #define WCSCOMPARE_ANCILLARY 0x0001
1746 #define WCSCOMPARE_TILING    0x0002
1747 #define WCSCOMPARE_CRPIX     0x0004
1748 
1749 
1750 extern const char *wcs_errmsg[];
1751 
1752 enum wcs_errmsg_enum {
1753   WCSERR_SUCCESS         =  0,	// Success.
1754   WCSERR_NULL_POINTER    =  1,	// Null wcsprm pointer passed.
1755   WCSERR_MEMORY          =  2,	// Memory allocation failed.
1756   WCSERR_SINGULAR_MTX    =  3,	// Linear transformation matrix is singular.
1757   WCSERR_BAD_CTYPE       =  4,	// Inconsistent or unrecognized coordinate
1758 				// axis type.
1759   WCSERR_BAD_PARAM       =  5,	// Invalid parameter value.
1760   WCSERR_BAD_COORD_TRANS =  6,	// Unrecognized coordinate transformation
1761 				// parameter.
1762   WCSERR_ILL_COORD_TRANS =  7,	// Ill-conditioned coordinate transformation
1763 				// parameter.
1764   WCSERR_BAD_PIX         =  8,	// One or more of the pixel coordinates were
1765 				// invalid.
1766   WCSERR_BAD_WORLD       =  9,	// One or more of the world coordinates were
1767 				// invalid.
1768   WCSERR_BAD_WORLD_COORD = 10,	// Invalid world coordinate.
1769   WCSERR_NO_SOLUTION     = 11,	// No solution found in the specified
1770 				// interval.
1771   WCSERR_BAD_SUBIMAGE    = 12,	// Invalid subimage specification.
1772   WCSERR_NON_SEPARABLE   = 13 	// Non-separable subimage coordinate system.
1773 };
1774 
1775 
1776 // Struct used for storing PVi_ma keywords.
1777 struct pvcard {
1778   int i;			// Axis number, as in PVi_ma (1-relative).
1779   int m;			// Parameter number, ditto  (0-relative).
1780   double value;			// Parameter value.
1781 };
1782 
1783 // Size of the pvcard struct in int units, used by the Fortran wrappers.
1784 #define PVLEN (sizeof(struct pvcard)/sizeof(int))
1785 
1786 // Struct used for storing PSi_ma keywords.
1787 struct pscard {
1788   int i;			// Axis number, as in PSi_ma (1-relative).
1789   int m;			// Parameter number, ditto  (0-relative).
1790   char value[72];		// Parameter value.
1791 };
1792 
1793 // Size of the pscard struct in int units, used by the Fortran wrappers.
1794 #define PSLEN (sizeof(struct pscard)/sizeof(int))
1795 
1796 // Struct used to hold additional auxiliary parameters.
1797 struct auxprm {
1798   double rsun_ref;              // Solar radius.
1799   double dsun_obs;              // Distance from Sun centre to observer.
1800   double crln_obs;              // Carrington heliographic lng of observer.
1801   double hgln_obs;              // Stonyhurst heliographic lng of observer.
1802   double hglt_obs;              // Heliographic latitude of observer.
1803 };
1804 
1805 // Size of the auxprm struct in int units, used by the Fortran wrappers.
1806 #define AUXLEN (sizeof(struct auxprm)/sizeof(int))
1807 
1808 
1809 struct wcsprm {
1810   // Initialization flag (see the prologue above).
1811   //--------------------------------------------------------------------------
1812   int    flag;			// Set to zero to force initialization.
1813 
1814   // FITS header keyvalues to be provided (see the prologue above).
1815   //--------------------------------------------------------------------------
1816   int    naxis;			// Number of axes (pixel and coordinate).
1817   double *crpix;		// CRPIXja keyvalues for each pixel axis.
1818   double *pc;			// PCi_ja  linear transformation matrix.
1819   double *cdelt;		// CDELTia keyvalues for each coord axis.
1820   double *crval;		// CRVALia keyvalues for each coord axis.
1821 
1822   char   (*cunit)[72];		// CUNITia keyvalues for each coord axis.
1823   char   (*ctype)[72];		// CTYPEia keyvalues for each coord axis.
1824 
1825   double lonpole;		// LONPOLEa keyvalue.
1826   double latpole;		// LATPOLEa keyvalue.
1827 
1828   double restfrq;		// RESTFRQa keyvalue.
1829   double restwav;		// RESTWAVa keyvalue.
1830 
1831   int    npv;			// Number of PVi_ma keywords, and the
1832   int    npvmax;		// number for which space was allocated.
1833   struct pvcard *pv;		// PVi_ma keywords for each i and m.
1834 
1835   int    nps;			// Number of PSi_ma keywords, and the
1836   int    npsmax;		// number for which space was allocated.
1837   struct pscard *ps;		// PSi_ma keywords for each i and m.
1838 
1839   // Alternative header keyvalues (see the prologue above).
1840   //--------------------------------------------------------------------------
1841   double *cd;			// CDi_ja linear transformation matrix.
1842   double *crota;		// CROTAi keyvalues for each coord axis.
1843   int    altlin;		// Alternative representations
1844 				//   Bit 0: PCi_ja  is present,
1845 				//   Bit 1: CDi_ja  is present,
1846 				//   Bit 2: CROTAi is present.
1847   int    velref;		// AIPS velocity code, VELREF.
1848 
1849   // Auxiliary coordinate system information of a general nature.  Not
1850   // used by WCSLIB.  Refer to the prologue comments above for a brief
1851   // explanation of these values.
1852   char   alt[4];
1853   int    colnum;
1854   int    *colax;
1855 				// Auxiliary coordinate axis information.
1856   char   (*cname)[72];
1857   double *crder;
1858   double *csyer;
1859   double *czphs;
1860   double *cperi;
1861 
1862   char   wcsname[72];
1863 				// Time reference system and measurement.
1864   char   timesys[72], trefpos[72], trefdir[72], plephem[72];
1865   char   timeunit[72];
1866   char   dateref[72];
1867   double mjdref[2];
1868   double timeoffs;
1869 				// Data timestamps and durations.
1870   char   dateobs[72], datebeg[72], dateavg[72], dateend[72];
1871   double mjdobs, mjdbeg, mjdavg, mjdend;
1872   double jepoch, bepoch;
1873   double tstart, tstop;
1874   double xposure, telapse;
1875 				// Timing accuracy.
1876   double timsyer, timrder;
1877   double timedel, timepixr;
1878 				// Spatial & celestial reference frame.
1879   double obsgeo[6];
1880   char   obsorbit[72];
1881   char   radesys[72];
1882   double equinox;
1883   char   specsys[72];
1884   char   ssysobs[72];
1885   double velosys;
1886   double zsource;
1887   char   ssyssrc[72];
1888   double velangl;
1889 
1890   // Additional auxiliary coordinate system information of a specialist
1891   // nature.  Not used by WCSLIB.  Refer to the prologue comments above.
1892   struct auxprm *aux;
1893 
1894   // Coordinate lookup tables (see the prologue above).
1895   //--------------------------------------------------------------------------
1896   int    ntab;			// Number of separate tables.
1897   int    nwtb;			// Number of wtbarr structs.
1898   struct tabprm *tab;		// Tabular transformation parameters.
1899   struct wtbarr *wtb;		// Array of wtbarr structs.
1900 
1901   //--------------------------------------------------------------------------
1902   // Information derived from the FITS header keyvalues by wcsset().
1903   //--------------------------------------------------------------------------
1904   char   lngtyp[8], lattyp[8];	// Celestial axis types, e.g. RA, DEC.
1905   int    lng, lat, spec;	// Longitude, latitude and spectral axis
1906 				// indices (0-relative).
1907   int    cubeface;		// True if there is a CUBEFACE axis.
1908   int    *types;		// Coordinate type codes for each axis.
1909 
1910   struct linprm lin;		//    Linear transformation parameters.
1911   struct celprm cel;		// Celestial transformation parameters.
1912   struct spcprm spc;		//  Spectral transformation parameters.
1913 
1914   //--------------------------------------------------------------------------
1915   //             THE REMAINDER OF THE WCSPRM STRUCT IS PRIVATE.
1916   //--------------------------------------------------------------------------
1917 
1918   // Error handling, if enabled.
1919   //--------------------------------------------------------------------------
1920   struct wcserr *err;
1921 
1922   // Memory management.
1923   //--------------------------------------------------------------------------
1924   int    m_flag, m_naxis;
1925   double *m_crpix, *m_pc, *m_cdelt, *m_crval;
1926   char  (*m_cunit)[72], (*m_ctype)[72];
1927   struct pvcard *m_pv;
1928   struct pscard *m_ps;
1929   double *m_cd, *m_crota;
1930   int    *m_colax;
1931   char  (*m_cname)[72];
1932   double *m_crder, *m_csyer, *m_czphs, *m_cperi;
1933   struct auxprm *m_aux;
1934   struct tabprm *m_tab;
1935   struct wtbarr *m_wtb;
1936 };
1937 
1938 // Size of the wcsprm struct in int units, used by the Fortran wrappers.
1939 #define WCSLEN (sizeof(struct wcsprm)/sizeof(int))
1940 
1941 
1942 int wcsnpv(int n);
1943 
1944 int wcsnps(int n);
1945 
1946 int wcsini(int alloc, int naxis, struct wcsprm *wcs);
1947 
1948 int wcsinit(int alloc, int naxis, struct wcsprm *wcs, int npvmax, int npsmax,
1949             int ndpmax);
1950 
1951 int wcsauxi(int alloc, struct wcsprm *wcs);
1952 
1953 int wcssub(int alloc, const struct wcsprm *wcssrc, int *nsub, int axes[],
1954            struct wcsprm *wcsdst);
1955 
1956 int wcscompare(int cmp, double tol, const struct wcsprm *wcs1,
1957                const struct wcsprm *wcs2, int *equal);
1958 
1959 int wcsfree(struct wcsprm *wcs);
1960 
1961 int wcsprt(const struct wcsprm *wcs);
1962 
1963 int wcsperr(const struct wcsprm *wcs, const char *prefix);
1964 
1965 int wcsbchk(struct wcsprm *wcs, int bounds);
1966 
1967 int wcsset(struct wcsprm *wcs);
1968 
1969 int wcsp2s(struct wcsprm *wcs, int ncoord, int nelem, const double pixcrd[],
1970            double imgcrd[], double phi[], double theta[], double world[],
1971            int stat[]);
1972 
1973 int wcss2p(struct wcsprm *wcs, int ncoord, int nelem, const double world[],
1974            double phi[], double theta[], double imgcrd[], double pixcrd[],
1975            int stat[]);
1976 
1977 int wcsmix(struct wcsprm *wcs, int mixpix, int mixcel, const double vspan[],
1978            double vstep, int viter, double world[], double phi[],
1979            double theta[], double imgcrd[], double pixcrd[]);
1980 
1981 int wcssptr(struct wcsprm *wcs, int *i, char ctype[9]);
1982 
1983 const char* wcslib_version(int vers[3]);
1984 
1985 // Defined mainly for backwards compatibility, use wcssub() instead.
1986 #define wcscopy(alloc, wcssrc, wcsdst) wcssub(alloc, wcssrc, 0x0, 0x0, wcsdst)
1987 
1988 
1989 // Deprecated.
1990 #define wcsini_errmsg wcs_errmsg
1991 #define wcssub_errmsg wcs_errmsg
1992 #define wcscopy_errmsg wcs_errmsg
1993 #define wcsfree_errmsg wcs_errmsg
1994 #define wcsprt_errmsg wcs_errmsg
1995 #define wcsset_errmsg wcs_errmsg
1996 #define wcsp2s_errmsg wcs_errmsg
1997 #define wcss2p_errmsg wcs_errmsg
1998 #define wcsmix_errmsg wcs_errmsg
1999 
2000 #ifdef __cplusplus
2001 #undef wtbarr
2002 }
2003 #endif
2004 
2005 #endif // WCSLIB_WCS
2006