1 /*** File libwcs/wcs.c
2 *** June 24, 2016
3 *** By Jessica Mink, jmink@cfa.harvard.edu
4 *** Harvard-Smithsonian Center for Astrophysics
5 *** Copyright (C) 1994-2016
6 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
7
8 This library is free software; you can redistribute it and/or
9 modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2 of the License, or (at your option) any later version.
12
13 This library is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with this library; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22 Correspondence concerning WCSTools should be addressed as follows:
23 Internet email: jmink@cfa.harvard.edu
24 Postal address: Jessica Mink
25 Smithsonian Astrophysical Observatory
26 60 Garden St.
27 Cambridge, MA 02138 USA
28
29 * Module: wcs.c (World Coordinate Systems)
30 * Purpose: Convert FITS WCS to pixels and vice versa:
31 * Subroutine: wcsxinit (cra,cdec,secpix,xrpix,yrpix,nxpix,nypix,rotate,equinox,epoch,proj)
32 * sets a WCS structure from arguments
33 * Subroutine: wcskinit (nxpix,nypix,ctype1,ctype2,crpix1,crpix2,crval1,crval2,
34 cd,cdelt1,cdelt2,crota,equinox,epoch)
35 * sets a WCS structure from keyword-based arguments
36 * Subroutine: wcsreset (wcs,crpix1,crpix2,crval1,crval2,cdelt1,cdelt2,crota,cd, equinox)
37 * resets an existing WCS structure from arguments
38 * Subroutine: wcsdeltset (wcs,cdelt1,cdelt2,crota) sets rotation and scaling
39 * Subroutine: wcscdset (wcs, cd) sets rotation and scaling from a CD matrix
40 * Subroutine: wcspcset (wcs,cdelt1,cdelt2,pc) sets rotation and scaling
41 * Subroutine: wcseqset (wcs, equinox) resets an existing WCS structure to new equinox
42 * Subroutine: iswcs(wcs) returns 1 if WCS structure is filled, else 0
43 * Subroutine: nowcs(wcs) returns 0 if WCS structure is filled, else 1
44 * Subroutine: wcscent (wcs) prints the image center and size in WCS units
45 * Subroutine: wcssize (wcs, cra, cdec, dra, ddec) returns image center and size
46 * Subroutine: wcsfull (wcs, cra, cdec, width, height) returns image center and size
47 * Subroutine: wcsrange (wcs, ra1, ra2, dec1, dec2) returns image coordinate limits
48
49 * Subroutine: wcsshift (wcs,cra,cdec) resets the center of a WCS structure
50 * Subroutine: wcsdist (x1,y1,x2,y2) compute angular distance between ra/dec or lat/long
51 * Subroutine: wcsdiff (x1,y1,x2,y2) compute angular distance between ra/dec or lat/long
52 * Subroutine: wcscominit (wcs,command) sets up a command format for execution by wcscom
53 * Subroutine: wcsoutinit (wcs,coor) sets up the coordinate system used by pix2wcs
54 * Subroutine: getwcsout (wcs) returns current output coordinate system used by pix2wcs
55 * Subroutine: wcsininit (wcs,coor) sets up the coordinate system used by wcs2pix
56 * Subroutine: getwcsin (wcs) returns current input coordinate system used by wcs2pix
57 * Subroutine: setwcsdeg(wcs, new) sets WCS output in degrees or hh:mm:ss
58 * Subroutine: getradecsys(wcs) returns current coordinate system type
59 * Subroutine: wcscom (wcs,file,x,y,wcstr) executes a command using the current world coordinates
60 * Subroutine: setwcslin (wcs, mode) sets output string mode for LINEAR
61 * Subroutine: pix2wcst (wcs,xpix,ypix,wcstring,lstr) pixels -> sky coordinate string
62 * Subroutine: pix2wcs (wcs,xpix,ypix,xpos,ypos) pixel coordinates -> sky coordinates
63 * Subroutine: wcsc2pix (wcs,xpos,ypos,coorsys,xpix,ypix,offscl) sky coordinates -> pixel coordinates
64 * Subroutine: wcs2pix (wcs,xpos,ypos,xpix,ypix,offscl) sky coordinates -> pixel coordinates
65 * Subroutine: wcszin (izpix) sets third dimension for pix2wcs() and pix2wcst()
66 * Subroutine: wcszout (wcs) returns third dimension from wcs2pix()
67 * Subroutine: setwcsfile (filename) Set file name for error messages
68 * Subroutine: setwcserr (errmsg) Set error message
69 * Subroutine: wcserr() Print error message
70 * Subroutine: setdefwcs (wcsproj) Set flag to choose AIPS or WCSLIB WCS subroutines
71 * Subroutine: getdefwcs() Get flag to switch between AIPS and WCSLIB subroutines
72 * Subroutine: savewcscoor (wcscoor)
73 * Subroutine: getwcscoor() Return preset output default coordinate system
74 * Subroutine: savewcscom (i, wcscom) Save specified WCS command
75 * Subroutine: setwcscom (wcs) Initialize WCS commands
76 * Subroutine: getwcscom (i) Return specified WCS command
77 * Subroutine: wcsfree (wcs) Free storage used by WCS structure
78 * Subroutine: freewcscom (wcs) Free storage used by WCS commands
79 * Subroutine: cpwcs (&header, cwcs)
80 */
81
82 #include <string.h> /* strstr, NULL */
83 #include <stdio.h> /* stderr */
84 #include <math.h>
85 #include "wcs.h"
86 #ifndef VMS
87 #include <stdlib.h>
88 #endif
89
90 static char wcserrmsg[80];
91 static char wcsfile[256]={""};
92 static void wcslibrot();
93 void wcsrotset();
94 static int wcsproj0 = 0;
95 static int izpix = 0;
96 static double zpix = 0.0;
97
98 void
wcsfree(wcs)99 wcsfree (wcs)
100 struct WorldCoor *wcs; /* WCS structure */
101 {
102 if (nowcs (wcs)) {
103
104 /* Free WCS structure if allocated but not filled */
105 if (wcs)
106 free (wcs);
107
108 return;
109 }
110
111 /* Free WCS on which this WCS depends */
112 if (wcs->wcs) {
113 wcsfree (wcs->wcs);
114 wcs->wcs = NULL;
115 }
116
117 freewcscom (wcs);
118 if (wcs->wcsname != NULL)
119 free (wcs->wcsname);
120 if (wcs->lin.imgpix != NULL)
121 free (wcs->lin.imgpix);
122 if (wcs->lin.piximg != NULL)
123 free (wcs->lin.piximg);
124 if (wcs->inv_x != NULL)
125 poly_end (wcs->inv_x);
126 if (wcs->inv_y != NULL)
127 poly_end (wcs->inv_y);
128 free (wcs);
129 return;
130 }
131
132 /* Set up a WCS structure from subroutine arguments */
133
134 struct WorldCoor *
wcsxinit(cra,cdec,secpix,xrpix,yrpix,nxpix,nypix,rotate,equinox,epoch,proj)135 wcsxinit (cra,cdec,secpix,xrpix,yrpix,nxpix,nypix,rotate,equinox,epoch,proj)
136
137 double cra; /* Center right ascension in degrees */
138 double cdec; /* Center declination in degrees */
139 double secpix; /* Number of arcseconds per pixel */
140 double xrpix; /* Reference pixel X coordinate */
141 double yrpix; /* Reference pixel X coordinate */
142 int nxpix; /* Number of pixels along x-axis */
143 int nypix; /* Number of pixels along y-axis */
144 double rotate; /* Rotation angle (clockwise positive) in degrees */
145 int equinox; /* Equinox of coordinates, 1950 and 2000 supported */
146 double epoch; /* Epoch of coordinates, used for FK4/FK5 conversion
147 * no effect if 0 */
148 char *proj; /* Projection */
149
150 {
151 struct WorldCoor *wcs;
152 double cdelt1, cdelt2;
153
154 wcs = (struct WorldCoor *) calloc (1, sizeof(struct WorldCoor));
155
156 /* Set WCSLIB flags so that structures will be reinitialized */
157 wcs->cel.flag = 0;
158 wcs->lin.flag = 0;
159 wcs->wcsl.flag = 0;
160
161 /* Image dimensions */
162 wcs->naxis = 2;
163 wcs->naxes = 2;
164 wcs->lin.naxis = 2;
165 wcs->nxpix = nxpix;
166 wcs->nypix = nypix;
167
168 wcs->wcsproj = wcsproj0;
169
170 wcs->crpix[0] = xrpix;
171 wcs->crpix[1] = yrpix;
172 wcs->xrefpix = wcs->crpix[0];
173 wcs->yrefpix = wcs->crpix[1];
174 wcs->lin.crpix = wcs->crpix;
175
176 wcs->crval[0] = cra;
177 wcs->crval[1] = cdec;
178 wcs->xref = wcs->crval[0];
179 wcs->yref = wcs->crval[1];
180 wcs->cel.ref[0] = wcs->crval[0];
181 wcs->cel.ref[1] = wcs->crval[1];
182 wcs->cel.ref[2] = 999.0;
183
184 strcpy (wcs->c1type,"RA");
185 strcpy (wcs->c2type,"DEC");
186
187 /* Allan Brighton: 28.4.98: for backward compat., remove leading "--" */
188 while (proj && *proj == '-')
189 proj++;
190 strncpy (wcs->ptype,proj, 3);
191 wcs->ptype[3] = (char) 0;
192 strcpy (wcs->ctype[0],"RA---");
193 strcpy (wcs->ctype[1],"DEC--");
194 strcat (wcs->ctype[0],proj);
195 strcat (wcs->ctype[1],proj);
196
197 if (wcstype (wcs, wcs->ctype[0], wcs->ctype[1])) {
198 wcsfree (wcs);
199 return (NULL);
200 }
201
202 /* Approximate world coordinate system from a known plate scale */
203 cdelt1 = -secpix / 3600.0;
204 cdelt2 = secpix / 3600.0;
205 wcsdeltset (wcs, cdelt1, cdelt2, rotate);
206 wcs->lin.cdelt = wcs->cdelt;
207 wcs->lin.pc = wcs->pc;
208
209 /* Coordinate reference frame and equinox */
210 wcs->equinox = (double) equinox;
211 if (equinox > 1980)
212 strcpy (wcs->radecsys,"FK5");
213 else
214 strcpy (wcs->radecsys,"FK4");
215 if (epoch > 0)
216 wcs->epoch = epoch;
217 else
218 wcs->epoch = 0.0;
219 wcs->wcson = 1;
220
221 wcs->syswcs = wcscsys (wcs->radecsys);
222 wcsoutinit (wcs, wcs->radecsys);
223 wcsininit (wcs, wcs->radecsys);
224 wcs->eqout = 0.0;
225 wcs->printsys = 1;
226 wcs->tabsys = 0;
227
228 /* Initialize special WCS commands */
229 setwcscom (wcs);
230
231 return (wcs);
232 }
233
234
235 /* Set up a WCS structure from subroutine arguments based on FITS keywords */
236
237 struct WorldCoor *
wcskinit(naxis1,naxis2,ctype1,ctype2,crpix1,crpix2,crval1,crval2,cd,cdelt1,cdelt2,crota,equinox,epoch)238 wcskinit (naxis1, naxis2, ctype1, ctype2, crpix1, crpix2, crval1, crval2,
239 cd, cdelt1, cdelt2, crota, equinox, epoch)
240
241 int naxis1; /* Number of pixels along x-axis */
242 int naxis2; /* Number of pixels along y-axis */
243 char *ctype1; /* FITS WCS projection for axis 1 */
244 char *ctype2; /* FITS WCS projection for axis 2 */
245 double crpix1, crpix2; /* Reference pixel coordinates */
246 double crval1, crval2; /* Coordinates at reference pixel in degrees */
247 double *cd; /* Rotation matrix, used if not NULL */
248 double cdelt1, cdelt2; /* scale in degrees/pixel, ignored if cd is not NULL */
249 double crota; /* Rotation angle in degrees, ignored if cd is not NULL */
250 int equinox; /* Equinox of coordinates, 1950 and 2000 supported */
251 double epoch; /* Epoch of coordinates, used for FK4/FK5 conversion
252 * no effect if 0 */
253 {
254 struct WorldCoor *wcs;
255
256 wcs = (struct WorldCoor *) calloc (1, sizeof(struct WorldCoor));
257
258 /* Set WCSLIB flags so that structures will be reinitialized */
259 wcs->cel.flag = 0;
260 wcs->lin.flag = 0;
261 wcs->wcsl.flag = 0;
262
263 /* Image dimensions */
264 wcs->naxis = 2;
265 wcs->naxes = 2;
266 wcs->lin.naxis = 2;
267 wcs->nxpix = naxis1;
268 wcs->nypix = naxis2;
269
270 wcs->wcsproj = wcsproj0;
271
272 wcs->crpix[0] = crpix1;
273 wcs->crpix[1] = crpix2;
274 wcs->xrefpix = wcs->crpix[0];
275 wcs->yrefpix = wcs->crpix[1];
276 wcs->lin.crpix = wcs->crpix;
277
278 if (wcstype (wcs, ctype1, ctype2)) {
279 wcsfree (wcs);
280 return (NULL);
281 }
282 if (wcs->latbase == 90)
283 crval2 = 90.0 - crval2;
284 else if (wcs->latbase == -90)
285 crval2 = crval2 - 90.0;
286
287 wcs->crval[0] = crval1;
288 wcs->crval[1] = crval2;
289 wcs->xref = wcs->crval[0];
290 wcs->yref = wcs->crval[1];
291 wcs->cel.ref[0] = wcs->crval[0];
292 wcs->cel.ref[1] = wcs->crval[1];
293 wcs->cel.ref[2] = 999.0;
294
295 if (cd != NULL)
296 wcscdset (wcs, cd);
297
298 else if (cdelt1 != 0.0)
299 wcsdeltset (wcs, cdelt1, cdelt2, crota);
300
301 else {
302 wcsdeltset (wcs, 1.0, 1.0, crota);
303 setwcserr ("WCSRESET: setting CDELT to 1");
304 }
305 wcs->lin.cdelt = wcs->cdelt;
306 wcs->lin.pc = wcs->pc;
307
308 /* Coordinate reference frame and equinox */
309 wcs->equinox = (double) equinox;
310 if (equinox > 1980)
311 strcpy (wcs->radecsys,"FK5");
312 else
313 strcpy (wcs->radecsys,"FK4");
314 if (epoch > 0)
315 wcs->epoch = epoch;
316 else
317 wcs->epoch = 0.0;
318 wcs->wcson = 1;
319
320 strcpy (wcs->radecout, wcs->radecsys);
321 wcs->syswcs = wcscsys (wcs->radecsys);
322 wcsoutinit (wcs, wcs->radecsys);
323 wcsininit (wcs, wcs->radecsys);
324 wcs->eqout = 0.0;
325 wcs->printsys = 1;
326 wcs->tabsys = 0;
327
328 /* Initialize special WCS commands */
329 setwcscom (wcs);
330
331 return (wcs);
332 }
333
334
335 /* Set projection in WCS structure from FITS keyword values */
336
337 int
wcstype(wcs,ctype1,ctype2)338 wcstype (wcs, ctype1, ctype2)
339
340 struct WorldCoor *wcs; /* World coordinate system structure */
341 char *ctype1; /* FITS WCS projection for axis 1 */
342 char *ctype2; /* FITS WCS projection for axis 2 */
343
344 {
345 int i, iproj;
346 int nctype = NWCSTYPE;
347 char ctypes[NWCSTYPE][4];
348 char dtypes[10][4];
349
350 /* Initialize projection types */
351 strcpy (ctypes[0], "LIN");
352 strcpy (ctypes[1], "AZP");
353 strcpy (ctypes[2], "SZP");
354 strcpy (ctypes[3], "TAN");
355 strcpy (ctypes[4], "SIN");
356 strcpy (ctypes[5], "STG");
357 strcpy (ctypes[6], "ARC");
358 strcpy (ctypes[7], "ZPN");
359 strcpy (ctypes[8], "ZEA");
360 strcpy (ctypes[9], "AIR");
361 strcpy (ctypes[10], "CYP");
362 strcpy (ctypes[11], "CAR");
363 strcpy (ctypes[12], "MER");
364 strcpy (ctypes[13], "CEA");
365 strcpy (ctypes[14], "COP");
366 strcpy (ctypes[15], "COD");
367 strcpy (ctypes[16], "COE");
368 strcpy (ctypes[17], "COO");
369 strcpy (ctypes[18], "BON");
370 strcpy (ctypes[19], "PCO");
371 strcpy (ctypes[20], "SFL");
372 strcpy (ctypes[21], "PAR");
373 strcpy (ctypes[22], "AIT");
374 strcpy (ctypes[23], "MOL");
375 strcpy (ctypes[24], "CSC");
376 strcpy (ctypes[25], "QSC");
377 strcpy (ctypes[26], "TSC");
378 strcpy (ctypes[27], "NCP");
379 strcpy (ctypes[28], "GLS");
380 strcpy (ctypes[29], "DSS");
381 strcpy (ctypes[30], "PLT");
382 strcpy (ctypes[31], "TNX");
383 strcpy (ctypes[32], "ZPX");
384 strcpy (ctypes[33], "TPV");
385
386 /* Initialize distortion types */
387 strcpy (dtypes[1], "SIP");
388
389 if (!strncmp (ctype1, "LONG",4))
390 strncpy (ctype1, "XLON",4);
391
392 strcpy (wcs->ctype[0], ctype1);
393
394 /* This is only to catch special non-standard projections */
395 strncpy (wcs->ptype, ctype1, 3);
396 wcs->ptype[3] = 0;
397
398 /* Linear coordinates */
399 if (!strncmp (ctype1,"LINEAR",6)) {
400 wcs->prjcode = WCS_LIN;
401 strcpy (wcs->c1type, "LIN");
402 strcpy (wcs->ptype, "LIN");
403 }
404
405 /* Pixel coordinates */
406 else if (!strncmp (ctype1,"PIXEL",6)) {
407 wcs->prjcode = WCS_PIX;
408 strcpy (wcs->c1type, "PIX");
409 strcpy (wcs->ptype, "PIX");
410 }
411
412 /*Detector pixel coordinates */
413 else if (strsrch (ctype1,"DET")) {
414 wcs->prjcode = WCS_PIX;
415 strcpy (wcs->c1type, "PIX");
416 strcpy (wcs->ptype, "PIX");
417 }
418
419 /* Set up right ascension, declination, latitude, or longitude */
420 else if (ctype1[0] == 'R' || ctype1[0] == 'D' ||
421 ctype1[0] == 'A' || ctype1[1] == 'L') {
422 wcs->c1type[0] = ctype1[0];
423 wcs->c1type[1] = ctype1[1];
424 if (ctype1[2] == '-') {
425 wcs->c1type[2] = 0;
426 iproj = 3;
427 }
428 else {
429 wcs->c1type[2] = ctype1[2];
430 iproj = 4;
431 if (ctype1[3] == '-') {
432 wcs->c1type[3] = 0;
433 }
434 else {
435 wcs->c1type[3] = ctype1[3];
436 wcs->c1type[4] = 0;
437 }
438 }
439 if (ctype1[iproj] == '-') iproj = iproj + 1;
440 if (ctype1[iproj] == '-') iproj = iproj + 1;
441 if (ctype1[iproj] == '-') iproj = iproj + 1;
442 if (ctype1[iproj] == '-') iproj = iproj + 1;
443 wcs->ptype[0] = ctype1[iproj];
444 wcs->ptype[1] = ctype1[iproj+1];
445 wcs->ptype[2] = ctype1[iproj+2];
446 wcs->ptype[3] = 0;
447 sprintf (wcs->ctype[0],"%-4s %3s",wcs->c1type,wcs->ptype);
448 for (i = 0; i < 8; i++)
449 if (wcs->ctype[0][i] == ' ') wcs->ctype[0][i] = '-';
450
451 /* Find projection type */
452 wcs->prjcode = 0; /* default type is linear */
453 for (i = 1; i < nctype; i++) {
454 if (!strncmp(wcs->ptype, ctypes[i], 3))
455 wcs->prjcode = i;
456 }
457
458 /* Handle "obsolete" NCP projection (now WCSLIB should be OK)
459 if (wcs->prjcode == WCS_NCP) {
460 if (wcs->wcsproj == WCS_BEST)
461 wcs->wcsproj = WCS_OLD;
462 else if (wcs->wcsproj == WCS_ALT)
463 wcs->wcsproj = WCS_NEW;
464 } */
465
466 /* Work around bug in WCSLIB handling of CAR projection
467 else if (wcs->prjcode == WCS_CAR) {
468 if (wcs->wcsproj == WCS_BEST)
469 wcs->wcsproj = WCS_OLD;
470 else if (wcs->wcsproj == WCS_ALT)
471 wcs->wcsproj = WCS_NEW;
472 } */
473
474 /* Work around bug in WCSLIB handling of COE projection
475 else if (wcs->prjcode == WCS_COE) {
476 if (wcs->wcsproj == WCS_BEST)
477 wcs->wcsproj = WCS_OLD;
478 else if (wcs->wcsproj == WCS_ALT)
479 wcs->wcsproj = WCS_NEW;
480 }
481
482 else if (wcs->wcsproj == WCS_BEST) */
483 if (wcs->wcsproj == WCS_BEST)
484 wcs->wcsproj = WCS_NEW;
485
486 else if (wcs->wcsproj == WCS_ALT)
487 wcs->wcsproj = WCS_OLD;
488
489 /* if (wcs->wcsproj == WCS_OLD && (
490 wcs->prjcode != WCS_STG && wcs->prjcode != WCS_AIT &&
491 wcs->prjcode != WCS_MER && wcs->prjcode != WCS_GLS &&
492 wcs->prjcode != WCS_ARC && wcs->prjcode != WCS_TAN &&
493 wcs->prjcode != WCS_TNX && wcs->prjcode != WCS_SIN &&
494 wcs->prjcode != WCS_PIX && wcs->prjcode != WCS_LIN &&
495 wcs->prjcode != WCS_CAR && wcs->prjcode != WCS_COE &&
496 wcs->prjcode != WCS_NCP && wcs->prjcode != WCS_ZPX))
497 wcs->wcsproj = WCS_NEW; */
498
499 /* Handle NOAO corrected TNX as uncorrected TAN if oldwcs is set */
500 if (wcs->wcsproj == WCS_OLD && wcs->prjcode == WCS_TNX) {
501 wcs->ctype[0][6] = 'A';
502 wcs->ctype[0][7] = 'N';
503 wcs->prjcode = WCS_TAN;
504 }
505
506 /* Handle NOAO corrected ZPX as uncorrected ZPN if oldwcs is set */
507 if (wcs->wcsproj == WCS_OLD && wcs->prjcode == WCS_ZPX) {
508 wcs->ctype[0][6] = 'P';
509 wcs->ctype[0][7] = 'N';
510 wcs->prjcode = WCS_ZPN;
511 }
512 }
513
514 /* If not sky coordinates, assume linear */
515 else {
516 wcs->prjcode = WCS_LIN;
517 strcpy (wcs->c1type, "LIN");
518 strcpy (wcs->ptype, "LIN");
519 return (0);
520 }
521
522 /* Second coordinate type */
523 if (!strncmp (ctype2, "NPOL",4)) {
524 ctype2[0] = ctype1[0];
525 strncpy (ctype2+1, "LAT",3);
526 wcs->latbase = 90;
527 strcpy (wcs->radecsys,"NPOLE");
528 wcs->syswcs = WCS_NPOLE;
529 }
530 else if (!strncmp (ctype2, "SPA-",4)) {
531 ctype2[0] = ctype1[0];
532 strncpy (ctype2+1, "LAT",3);
533 wcs->latbase = -90;
534 strcpy (wcs->radecsys,"SPA");
535 wcs->syswcs = WCS_SPA;
536 }
537 else
538 wcs->latbase = 0;
539 strcpy (wcs->ctype[1], ctype2);
540
541 /* Linear coordinates */
542 if (!strncmp (ctype2,"LINEAR",6)) {
543 wcs->prjcode = WCS_LIN;
544 strcpy (wcs->c2type, "LIN");
545 }
546
547 /* Pixel coordinates */
548 else if (!strncmp (ctype2,"PIXEL",6)) {
549 wcs->prjcode = WCS_PIX;
550 strcpy (wcs->c2type, "PIX");
551 }
552
553 /* Detector coordinates */
554 else if (!strncmp (ctype2,"DET",3)) {
555 wcs->prjcode = WCS_PIX;
556 strcpy (wcs->c2type, "PIX");
557 }
558
559 /* Set up right ascension, declination, latitude, or longitude */
560 else if (ctype2[0] == 'R' || ctype2[0] == 'D' ||
561 ctype2[0] == 'A' || ctype2[1] == 'L') {
562 wcs->c2type[0] = ctype2[0];
563 wcs->c2type[1] = ctype2[1];
564 if (ctype2[2] == '-') {
565 wcs->c2type[2] = 0;
566 iproj = 3;
567 }
568 else {
569 wcs->c2type[2] = ctype2[2];
570 iproj = 4;
571 if (ctype2[3] == '-') {
572 wcs->c2type[3] = 0;
573 }
574 else {
575 wcs->c2type[3] = ctype2[3];
576 wcs->c2type[4] = 0;
577 }
578 }
579 if (ctype2[iproj] == '-') iproj = iproj + 1;
580 if (ctype2[iproj] == '-') iproj = iproj + 1;
581 if (ctype2[iproj] == '-') iproj = iproj + 1;
582 if (ctype2[iproj] == '-') iproj = iproj + 1;
583
584 if (!strncmp (ctype1, "DEC", 3) ||
585 !strncmp (ctype1+1, "LAT", 3))
586 wcs->coorflip = 1;
587 else
588 wcs->coorflip = 0;
589 if (ctype2[1] == 'L' || ctype2[0] == 'A') {
590 wcs->degout = 1;
591 wcs->ndec = 5;
592 }
593 else {
594 wcs->degout = 0;
595 wcs->ndec = 3;
596 }
597 sprintf (wcs->ctype[1],"%-4s %3s",wcs->c2type,wcs->ptype);
598 for (i = 0; i < 8; i++)
599 if (wcs->ctype[1][i] == ' ') wcs->ctype[1][i] = '-';
600 }
601
602 /* If not sky coordinates, assume linear */
603 else {
604 strcpy (wcs->c2type, "LIN");
605 wcs->prjcode = WCS_LIN;
606 }
607
608 /* Set distortion code from CTYPE1 extension */
609 setdistcode (wcs, ctype1);
610
611 return (0);
612 }
613
614
615 int
wcsreset(wcs,crpix1,crpix2,crval1,crval2,cdelt1,cdelt2,crota,cd)616 wcsreset (wcs, crpix1, crpix2, crval1, crval2, cdelt1, cdelt2, crota, cd)
617
618 struct WorldCoor *wcs; /* World coordinate system data structure */
619 double crpix1, crpix2; /* Reference pixel coordinates */
620 double crval1, crval2; /* Coordinates at reference pixel in degrees */
621 double cdelt1, cdelt2; /* scale in degrees/pixel, ignored if cd is not NULL */
622 double crota; /* Rotation angle in degrees, ignored if cd is not NULL */
623 double *cd; /* Rotation matrix, used if not NULL */
624 {
625
626 if (nowcs (wcs))
627 return (-1);
628
629 /* Set WCSLIB flags so that structures will be reinitialized */
630 wcs->cel.flag = 0;
631 wcs->lin.flag = 0;
632 wcs->wcsl.flag = 0;
633
634 /* Reference pixel coordinates and WCS value */
635 wcs->crpix[0] = crpix1;
636 wcs->crpix[1] = crpix2;
637 wcs->xrefpix = wcs->crpix[0];
638 wcs->yrefpix = wcs->crpix[1];
639 wcs->lin.crpix = wcs->crpix;
640
641 wcs->crval[0] = crval1;
642 wcs->crval[1] = crval2;
643 wcs->xref = wcs->crval[0];
644 wcs->yref = wcs->crval[1];
645 if (wcs->coorflip) {
646 wcs->cel.ref[1] = wcs->crval[0];
647 wcs->cel.ref[0] = wcs->crval[1];
648 }
649 else {
650 wcs->cel.ref[0] = wcs->crval[0];
651 wcs->cel.ref[1] = wcs->crval[1];
652 }
653 /* Keep ref[2] and ref[3] from input */
654
655 /* Initialize to no plate fit */
656 wcs->ncoeff1 = 0;
657 wcs->ncoeff2 = 0;
658
659 if (cd != NULL)
660 wcscdset (wcs, cd);
661
662 else if (cdelt1 != 0.0)
663 wcsdeltset (wcs, cdelt1, cdelt2, crota);
664
665 else {
666 wcs->xinc = 1.0;
667 wcs->yinc = 1.0;
668 setwcserr ("WCSRESET: setting CDELT to 1");
669 }
670
671 /* Coordinate reference frame, equinox, and epoch */
672 if (!strncmp (wcs->ptype,"LIN",3) ||
673 !strncmp (wcs->ptype,"PIX",3))
674 wcs->degout = -1;
675
676 wcs->wcson = 1;
677 return (0);
678 }
679
680 void
wcseqset(wcs,equinox)681 wcseqset (wcs, equinox)
682
683 struct WorldCoor *wcs; /* World coordinate system data structure */
684 double equinox; /* Desired equinox as fractional year */
685 {
686
687 if (nowcs (wcs))
688 return;
689
690 /* Leave WCS alone if already at desired equinox */
691 if (wcs->equinox == equinox)
692 return;
693
694 /* Convert center from B1950 (FK4) to J2000 (FK5) */
695 if (equinox == 2000.0 && wcs->equinox == 1950.0) {
696 if (wcs->coorflip) {
697 fk425e (&wcs->crval[1], &wcs->crval[0], wcs->epoch);
698 wcs->cel.ref[1] = wcs->crval[0];
699 wcs->cel.ref[0] = wcs->crval[1];
700 }
701 else {
702 fk425e (&wcs->crval[0], &wcs->crval[1], wcs->epoch);
703 wcs->cel.ref[0] = wcs->crval[0];
704 wcs->cel.ref[1] = wcs->crval[1];
705 }
706 wcs->xref = wcs->crval[0];
707 wcs->yref = wcs->crval[1];
708 wcs->equinox = 2000.0;
709 strcpy (wcs->radecsys, "FK5");
710 wcs->syswcs = WCS_J2000;
711 wcs->cel.flag = 0;
712 wcs->wcsl.flag = 0;
713 }
714
715 /* Convert center from J2000 (FK5) to B1950 (FK4) */
716 else if (equinox == 1950.0 && wcs->equinox == 2000.0) {
717 if (wcs->coorflip) {
718 fk524e (&wcs->crval[1], &wcs->crval[0], wcs->epoch);
719 wcs->cel.ref[1] = wcs->crval[0];
720 wcs->cel.ref[0] = wcs->crval[1];
721 }
722 else {
723 fk524e (&wcs->crval[0], &wcs->crval[1], wcs->epoch);
724 wcs->cel.ref[0] = wcs->crval[0];
725 wcs->cel.ref[1] = wcs->crval[1];
726 }
727 wcs->xref = wcs->crval[0];
728 wcs->yref = wcs->crval[1];
729 wcs->equinox = 1950.0;
730 strcpy (wcs->radecsys, "FK4");
731 wcs->syswcs = WCS_B1950;
732 wcs->cel.flag = 0;
733 wcs->wcsl.flag = 0;
734 }
735 wcsoutinit (wcs, wcs->radecsys);
736 wcsininit (wcs, wcs->radecsys);
737 return;
738 }
739
740
741 /* Set scale and rotation in WCS structure */
742
743 void
wcscdset(wcs,cd)744 wcscdset (wcs, cd)
745
746 struct WorldCoor *wcs; /* World coordinate system structure */
747 double *cd; /* CD matrix, ignored if NULL */
748 {
749 double tcd;
750
751 if (cd == NULL)
752 return;
753
754 wcs->rotmat = 1;
755 wcs->cd[0] = cd[0];
756 wcs->cd[1] = cd[1];
757 wcs->cd[2] = cd[2];
758 wcs->cd[3] = cd[3];
759 (void) matinv (2, wcs->cd, wcs->dc);
760
761 /* Compute scale */
762 wcs->xinc = sqrt (cd[0]*cd[0] + cd[2]*cd[2]);
763 wcs->yinc = sqrt (cd[1]*cd[1] + cd[3]*cd[3]);
764
765 /* Deal with x=Dec/y=RA case */
766 if (wcs->coorflip) {
767 tcd = cd[1];
768 cd[1] = -cd[2];
769 cd[2] = -tcd;
770 }
771 wcslibrot (wcs);
772 wcs->wcson = 1;
773
774 /* Compute image rotation */
775 wcsrotset (wcs);
776
777 wcs->cdelt[0] = wcs->xinc;
778 wcs->cdelt[1] = wcs->yinc;
779
780 return;
781 }
782
783
784 /* Set scale and rotation in WCS structure from axis scale and rotation */
785
786 void
wcsdeltset(wcs,cdelt1,cdelt2,crota)787 wcsdeltset (wcs, cdelt1, cdelt2, crota)
788
789 struct WorldCoor *wcs; /* World coordinate system structure */
790 double cdelt1; /* degrees/pixel in first axis (or both axes) */
791 double cdelt2; /* degrees/pixel in second axis if nonzero */
792 double crota; /* Rotation counterclockwise in degrees */
793 {
794 double *pci;
795 double crot, srot;
796 int i, j, naxes;
797
798 naxes = wcs->naxis;
799 if (naxes > 2)
800 naxes = 2;
801 wcs->cdelt[0] = cdelt1;
802 if (cdelt2 != 0.0)
803 wcs->cdelt[1] = cdelt2;
804 else
805 wcs->cdelt[1] = cdelt1;
806 wcs->xinc = wcs->cdelt[0];
807 wcs->yinc = wcs->cdelt[1];
808 pci = wcs->pc;
809 for (i = 0; i < naxes; i++) {
810 for (j = 0; j < naxes; j++) {
811 if (i ==j)
812 *pci = 1.0;
813 else
814 *pci = 0.0;
815 pci++;
816 }
817 }
818 wcs->rotmat = 0;
819
820 /* If image is reversed, value of CROTA is flipped, too */
821 wcs->rot = crota;
822 if (wcs->rot < 0.0)
823 wcs->rot = wcs->rot + 360.0;
824 if (wcs->rot >= 360.0)
825 wcs->rot = wcs->rot - 360.0;
826 crot = cos (degrad(wcs->rot));
827 if (cdelt1 * cdelt2 > 0)
828 srot = sin (-degrad(wcs->rot));
829 else
830 srot = sin (degrad(wcs->rot));
831
832 /* Set CD matrix */
833 wcs->cd[0] = wcs->cdelt[0] * crot;
834 if (wcs->cdelt[0] < 0)
835 wcs->cd[1] = -fabs (wcs->cdelt[1]) * srot;
836 else
837 wcs->cd[1] = fabs (wcs->cdelt[1]) * srot;
838 if (wcs->cdelt[1] < 0)
839 wcs->cd[2] = fabs (wcs->cdelt[0]) * srot;
840 else
841 wcs->cd[2] = -fabs (wcs->cdelt[0]) * srot;
842 wcs->cd[3] = wcs->cdelt[1] * crot;
843 (void) matinv (2, wcs->cd, wcs->dc);
844
845 /* Set rotation matrix */
846 wcslibrot (wcs);
847
848 /* Set image rotation and mirroring */
849 if (wcs->coorflip) {
850 if (wcs->cdelt[0] < 0 && wcs->cdelt[1] > 0) {
851 wcs->imflip = 1;
852 wcs->imrot = wcs->rot - 90.0;
853 if (wcs->imrot < -180.0) wcs->imrot = wcs->imrot + 360.0;
854 wcs->pa_north = wcs->rot;
855 wcs->pa_east = wcs->rot - 90.0;
856 if (wcs->pa_east < -180.0) wcs->pa_east = wcs->pa_east + 360.0;
857 }
858 else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] < 0) {
859 wcs->imflip = 1;
860 wcs->imrot = wcs->rot + 90.0;
861 if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
862 wcs->pa_north = wcs->rot;
863 wcs->pa_east = wcs->rot - 90.0;
864 if (wcs->pa_east < -180.0) wcs->pa_east = wcs->pa_east + 360.0;
865 }
866 else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] > 0) {
867 wcs->imflip = 0;
868 wcs->imrot = wcs->rot + 90.0;
869 if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
870 wcs->pa_north = wcs->imrot;
871 wcs->pa_east = wcs->rot + 90.0;
872 if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
873 }
874 else if (wcs->cdelt[0] < 0 && wcs->cdelt[1] < 0) {
875 wcs->imflip = 0;
876 wcs->imrot = wcs->rot - 90.0;
877 if (wcs->imrot < -180.0) wcs->imrot = wcs->imrot + 360.0;
878 wcs->pa_north = wcs->imrot;
879 wcs->pa_east = wcs->rot + 90.0;
880 if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
881 }
882 }
883 else {
884 if (wcs->cdelt[0] < 0 && wcs->cdelt[1] > 0) {
885 wcs->imflip = 0;
886 wcs->imrot = wcs->rot;
887 wcs->pa_north = wcs->rot + 90.0;
888 if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
889 wcs->pa_east = wcs->rot + 180.0;
890 if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
891 }
892 else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] < 0) {
893 wcs->imflip = 0;
894 wcs->imrot = wcs->rot + 180.0;
895 if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
896 wcs->pa_north = wcs->imrot + 90.0;
897 if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
898 wcs->pa_east = wcs->imrot + 180.0;
899 if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
900 }
901 else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] > 0) {
902 wcs->imflip = 1;
903 wcs->imrot = -wcs->rot;
904 wcs->pa_north = wcs->imrot + 90.0;
905 if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
906 wcs->pa_east = wcs->rot;
907 }
908 else if (wcs->cdelt[0] < 0 && wcs->cdelt[1] < 0) {
909 wcs->imflip = 1;
910 wcs->imrot = wcs->rot + 180.0;
911 if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
912 wcs->pa_north = wcs->imrot + 90.0;
913 if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
914 wcs->pa_east = wcs->rot + 90.0;
915 if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
916 }
917 }
918
919 return;
920 }
921
922
923 /* Set scale and rotation in WCS structure */
924
925 void
wcspcset(wcs,cdelt1,cdelt2,pc)926 wcspcset (wcs, cdelt1, cdelt2, pc)
927
928 struct WorldCoor *wcs; /* World coordinate system structure */
929 double cdelt1; /* degrees/pixel in first axis (or both axes) */
930 double cdelt2; /* degrees/pixel in second axis if nonzero */
931 double *pc; /* Rotation matrix, ignored if NULL */
932 {
933 double *pci, *pc0i;
934 int i, j, naxes;
935
936 if (pc == NULL)
937 return;
938
939 naxes = wcs->naxis;
940 /* if (naxes > 2)
941 naxes = 2; */
942 if (naxes < 1 || naxes > 9) {
943 naxes = wcs->naxes;
944 wcs->naxis = naxes;
945 }
946 wcs->cdelt[0] = cdelt1;
947 if (cdelt2 != 0.0)
948 wcs->cdelt[1] = cdelt2;
949 else
950 wcs->cdelt[1] = cdelt1;
951 wcs->xinc = wcs->cdelt[0];
952 wcs->yinc = wcs->cdelt[1];
953
954 /* Set rotation matrix */
955 pci = wcs->pc;
956 pc0i = pc;
957 for (i = 0; i < naxes; i++) {
958 for (j = 0; j < naxes; j++) {
959 *pci = *pc0i;
960 pci++;
961 pc0i++;
962 }
963 }
964
965 /* Set CD matrix */
966 if (naxes > 1) {
967 wcs->cd[0] = pc[0] * wcs->cdelt[0];
968 wcs->cd[1] = pc[1] * wcs->cdelt[0];
969 wcs->cd[2] = pc[naxes] * wcs->cdelt[1];
970 wcs->cd[3] = pc[naxes+1] * wcs->cdelt[1];
971 }
972 else {
973 wcs->cd[0] = pc[0] * wcs->cdelt[0];
974 wcs->cd[1] = 0.0;
975 wcs->cd[2] = 0.0;
976 wcs->cd[3] = 1.0;
977 }
978 (void) matinv (2, wcs->cd, wcs->dc);
979 wcs->rotmat = 1;
980
981 (void)linset (&wcs->lin);
982 wcs->wcson = 1;
983
984 wcsrotset (wcs);
985
986 return;
987 }
988
989
990 /* Set up rotation matrix for WCSLIB projection subroutines */
991
992 static void
wcslibrot(wcs)993 wcslibrot (wcs)
994
995 struct WorldCoor *wcs; /* World coordinate system structure */
996
997 {
998 int i, mem, naxes;
999
1000 naxes = wcs->naxis;
1001 if (naxes > 2)
1002 naxes = 2;
1003 if (naxes < 1 || naxes > 9) {
1004 naxes = wcs->naxes;
1005 wcs->naxis = naxes;
1006 }
1007 mem = naxes * naxes * sizeof(double);
1008 if (wcs->lin.piximg == NULL)
1009 wcs->lin.piximg = (double*)malloc(mem);
1010 if (wcs->lin.piximg != NULL) {
1011 if (wcs->lin.imgpix == NULL)
1012 wcs->lin.imgpix = (double*)malloc(mem);
1013 if (wcs->lin.imgpix != NULL) {
1014 wcs->lin.flag = LINSET;
1015 if (naxes == 2) {
1016 for (i = 0; i < 4; i++) {
1017 wcs->lin.piximg[i] = wcs->cd[i];
1018 }
1019 }
1020 else if (naxes == 3) {
1021 for (i = 0; i < 9; i++)
1022 wcs->lin.piximg[i] = 0.0;
1023 wcs->lin.piximg[0] = wcs->cd[0];
1024 wcs->lin.piximg[1] = wcs->cd[1];
1025 wcs->lin.piximg[3] = wcs->cd[2];
1026 wcs->lin.piximg[4] = wcs->cd[3];
1027 wcs->lin.piximg[8] = 1.0;
1028 }
1029 else if (naxes == 4) {
1030 for (i = 0; i < 16; i++)
1031 wcs->lin.piximg[i] = 0.0;
1032 wcs->lin.piximg[0] = wcs->cd[0];
1033 wcs->lin.piximg[1] = wcs->cd[1];
1034 wcs->lin.piximg[4] = wcs->cd[2];
1035 wcs->lin.piximg[5] = wcs->cd[3];
1036 wcs->lin.piximg[10] = 1.0;
1037 wcs->lin.piximg[15] = 1.0;
1038 }
1039 (void) matinv (naxes, wcs->lin.piximg, wcs->lin.imgpix);
1040 wcs->lin.crpix = wcs->crpix;
1041 wcs->lin.cdelt = wcs->cdelt;
1042 wcs->lin.pc = wcs->pc;
1043 wcs->lin.flag = LINSET;
1044 }
1045 }
1046 return;
1047 }
1048
1049
1050 /* Compute image rotation */
1051
1052 void
wcsrotset(wcs)1053 wcsrotset (wcs)
1054
1055 struct WorldCoor *wcs; /* World coordinate system structure */
1056 {
1057 int off;
1058 double cra, cdec, xc, xn, xe, yc, yn, ye;
1059
1060 /* If image is one-dimensional, leave rotation angle alone */
1061 if (wcs->nxpix < 1.5 || wcs->nypix < 1.5) {
1062 wcs->imrot = wcs->rot;
1063 wcs->pa_north = wcs->rot + 90.0;
1064 wcs->pa_east = wcs->rot + 180.0;
1065 return;
1066 }
1067
1068
1069 /* Do not try anything if image is LINEAR (not Cartesian projection) */
1070 if (wcs->syswcs == WCS_LINEAR)
1071 return;
1072
1073 wcs->xinc = fabs (wcs->xinc);
1074 wcs->yinc = fabs (wcs->yinc);
1075
1076 /* Compute position angles of North and East in image */
1077 xc = wcs->xrefpix;
1078 yc = wcs->yrefpix;
1079 pix2wcs (wcs, xc, yc, &cra, &cdec);
1080 if (wcs->coorflip) {
1081 wcs2pix (wcs, cra+wcs->yinc, cdec, &xe, &ye, &off);
1082 wcs2pix (wcs, cra, cdec+wcs->xinc, &xn, &yn, &off);
1083 }
1084 else {
1085 wcs2pix (wcs, cra+wcs->xinc, cdec, &xe, &ye, &off);
1086 wcs2pix (wcs, cra, cdec+wcs->yinc, &xn, &yn, &off);
1087 }
1088 wcs->pa_north = raddeg (atan2 (yn-yc, xn-xc));
1089 if (wcs->pa_north < -90.0)
1090 wcs->pa_north = wcs->pa_north + 360.0;
1091 wcs->pa_east = raddeg (atan2 (ye-yc, xe-xc));
1092 if (wcs->pa_east < -90.0)
1093 wcs->pa_east = wcs->pa_east + 360.0;
1094
1095 /* Compute image rotation angle from North */
1096 if (wcs->pa_north < -90.0)
1097 wcs->imrot = 270.0 + wcs->pa_north;
1098 else
1099 wcs->imrot = wcs->pa_north - 90.0;
1100
1101 /* Compute CROTA */
1102 if (wcs->coorflip) {
1103 wcs->rot = wcs->imrot + 90.0;
1104 if (wcs->rot < 0.0)
1105 wcs->rot = wcs->rot + 360.0;
1106 }
1107 else
1108 wcs->rot = wcs->imrot;
1109 if (wcs->rot < 0.0)
1110 wcs->rot = wcs->rot + 360.0;
1111 if (wcs->rot >= 360.0)
1112 wcs->rot = wcs->rot - 360.0;
1113
1114 /* Set image mirror flag based on axis orientation */
1115 wcs->imflip = 0;
1116 if (wcs->pa_east - wcs->pa_north < -80.0 &&
1117 wcs->pa_east - wcs->pa_north > -100.0)
1118 wcs->imflip = 1;
1119 if (wcs->pa_east - wcs->pa_north < 280.0 &&
1120 wcs->pa_east - wcs->pa_north > 260.0)
1121 wcs->imflip = 1;
1122 if (wcs->pa_north - wcs->pa_east > 80.0 &&
1123 wcs->pa_north - wcs->pa_east < 100.0)
1124 wcs->imflip = 1;
1125 if (wcs->coorflip) {
1126 if (wcs->imflip)
1127 wcs->yinc = -wcs->yinc;
1128 }
1129 else {
1130 if (!wcs->imflip)
1131 wcs->xinc = -wcs->xinc;
1132 }
1133
1134 return;
1135 }
1136
1137
1138 /* Return 1 if WCS structure is filled, else 0 */
1139
1140 int
iswcs(wcs)1141 iswcs (wcs)
1142
1143 struct WorldCoor *wcs; /* World coordinate system structure */
1144
1145 {
1146 if (wcs == NULL)
1147 return (0);
1148 else
1149 return (wcs->wcson);
1150 }
1151
1152
1153 /* Return 0 if WCS structure is filled, else 1 */
1154
1155 int
nowcs(wcs)1156 nowcs (wcs)
1157
1158 struct WorldCoor *wcs; /* World coordinate system structure */
1159
1160 {
1161 if (wcs == NULL)
1162 return (1);
1163 else
1164 return (!wcs->wcson);
1165 }
1166
1167
1168 /* Reset the center of a WCS structure */
1169
1170 void
wcsshift(wcs,rra,rdec,coorsys)1171 wcsshift (wcs,rra,rdec,coorsys)
1172
1173 struct WorldCoor *wcs; /* World coordinate system structure */
1174 double rra; /* Reference pixel right ascension in degrees */
1175 double rdec; /* Reference pixel declination in degrees */
1176 char *coorsys; /* FK4 or FK5 coordinates (1950 or 2000) */
1177
1178 {
1179 if (nowcs (wcs))
1180 return;
1181
1182 /* Approximate world coordinate system from a known plate scale */
1183 wcs->crval[0] = rra;
1184 wcs->crval[1] = rdec;
1185 wcs->xref = wcs->crval[0];
1186 wcs->yref = wcs->crval[1];
1187
1188
1189 /* Coordinate reference frame */
1190 strcpy (wcs->radecsys,coorsys);
1191 wcs->syswcs = wcscsys (coorsys);
1192 if (wcs->syswcs == WCS_B1950)
1193 wcs->equinox = 1950.0;
1194 else
1195 wcs->equinox = 2000.0;
1196
1197 return;
1198 }
1199
1200 /* Print position of WCS center, if WCS is set */
1201
1202 void
wcscent(wcs)1203 wcscent (wcs)
1204
1205 struct WorldCoor *wcs; /* World coordinate system structure */
1206
1207 {
1208 double xpix,ypix, xpos1, xpos2, ypos1, ypos2;
1209 char wcstring[32];
1210 double width, height, secpix, secpixh, secpixw;
1211 int lstr = 32;
1212
1213 if (nowcs (wcs))
1214 (void)fprintf (stderr,"No WCS information available\n");
1215 else {
1216 if (wcs->prjcode == WCS_DSS)
1217 (void)fprintf (stderr,"WCS plate center %s\n", wcs->center);
1218 xpix = 0.5 * wcs->nxpix;
1219 ypix = 0.5 * wcs->nypix;
1220 (void) pix2wcst (wcs,xpix,ypix,wcstring, lstr);
1221 (void)fprintf (stderr,"WCS center %s %s %s %s at pixel (%.2f,%.2f)\n",
1222 wcs->ctype[0],wcs->ctype[1],wcstring,wcs->ptype,xpix,ypix);
1223
1224 /* Image width */
1225 (void) pix2wcs (wcs,1.0,ypix,&xpos1,&ypos1);
1226 (void) pix2wcs (wcs,wcs->nxpix,ypix,&xpos2,&ypos2);
1227 if (wcs->syswcs == WCS_LINEAR) {
1228 width = xpos2 - xpos1;
1229 if (width < 100.0)
1230 (void)fprintf (stderr, "WCS width = %.5f %s ",width, wcs->units[0]);
1231 else
1232 (void)fprintf (stderr, "WCS width = %.3f %s ",width, wcs->units[0]);
1233 }
1234 else {
1235 width = wcsdist (xpos1,ypos1,xpos2,ypos2);
1236 if (width < 1/60.0)
1237 (void)fprintf (stderr, "WCS width = %.2f arcsec ",width*3600.0);
1238 else if (width < 1.0)
1239 (void)fprintf (stderr, "WCS width = %.2f arcmin ",width*60.0);
1240 else
1241 (void)fprintf (stderr, "WCS width = %.3f degrees ",width);
1242 }
1243 secpixw = width / (wcs->nxpix - 1.0);
1244
1245 /* Image height */
1246 (void) pix2wcs (wcs,xpix,1.0,&xpos1,&ypos1);
1247 (void) pix2wcs (wcs,xpix,wcs->nypix,&xpos2,&ypos2);
1248 if (wcs->syswcs == WCS_LINEAR) {
1249 height = ypos2 - ypos1;
1250 if (height < 100.0)
1251 (void)fprintf (stderr, " height = %.5f %s ",height, wcs->units[1]);
1252 else
1253 (void)fprintf (stderr, " height = %.3f %s ",height, wcs->units[1]);
1254 }
1255 else {
1256 height = wcsdist (xpos1,ypos1,xpos2,ypos2);
1257 if (height < 1/60.0)
1258 (void) fprintf (stderr, " height = %.2f arcsec",height*3600.0);
1259 else if (height < 1.0)
1260 (void) fprintf (stderr, " height = %.2f arcmin",height*60.0);
1261 else
1262 (void) fprintf (stderr, " height = %.3f degrees",height);
1263 }
1264 secpixh = height / (wcs->nypix - 1.0);
1265
1266 /* Image scale */
1267 if (wcs->syswcs == WCS_LINEAR) {
1268 (void) fprintf (stderr,"\n");
1269 (void) fprintf (stderr,"WCS %.5f %s/pixel, %.5f %s/pixel\n",
1270 wcs->xinc,wcs->units[0],wcs->yinc,wcs->units[1]);
1271 }
1272 else {
1273 if (wcs->xinc != 0.0 && wcs->yinc != 0.0)
1274 secpix = (fabs(wcs->xinc) + fabs(wcs->yinc)) * 0.5 * 3600.0;
1275 else if (secpixh > 0.0 && secpixw > 0.0)
1276 secpix = (secpixw + secpixh) * 0.5 * 3600.0;
1277 else if (wcs->xinc != 0.0 || wcs->yinc != 0.0)
1278 secpix = (fabs(wcs->xinc) + fabs(wcs->yinc)) * 3600.0;
1279 else
1280 secpix = (secpixw + secpixh) * 3600.0;
1281 if (secpix < 100.0)
1282 (void) fprintf (stderr, " %.3f arcsec/pixel\n",secpix);
1283 else if (secpix < 3600.0)
1284 (void) fprintf (stderr, " %.3f arcmin/pixel\n",secpix/60.0);
1285 else
1286 (void) fprintf (stderr, " %.3f degrees/pixel\n",secpix/3600.0);
1287 }
1288 }
1289 return;
1290 }
1291
1292 /* Return RA and Dec of image center, plus size in RA and Dec */
1293
1294 void
wcssize(wcs,cra,cdec,dra,ddec)1295 wcssize (wcs, cra, cdec, dra, ddec)
1296
1297 struct WorldCoor *wcs; /* World coordinate system structure */
1298 double *cra; /* Right ascension of image center (deg) (returned) */
1299 double *cdec; /* Declination of image center (deg) (returned) */
1300 double *dra; /* Half-width in right ascension (deg) (returned) */
1301 double *ddec; /* Half-width in declination (deg) (returned) */
1302
1303 {
1304 double width, height;
1305
1306 /* Find right ascension and declination of coordinates */
1307 if (iswcs(wcs)) {
1308 wcsfull (wcs, cra, cdec, &width, &height);
1309 *dra = 0.5 * width / cos (degrad (*cdec));
1310 *ddec = 0.5 * height;
1311 }
1312 else {
1313 *cra = 0.0;
1314 *cdec = 0.0;
1315 *dra = 0.0;
1316 *ddec = 0.0;
1317 }
1318 return;
1319 }
1320
1321
1322 /* Return RA and Dec of image center, plus size in degrees */
1323
1324 void
wcsfull(wcs,cra,cdec,width,height)1325 wcsfull (wcs, cra, cdec, width, height)
1326
1327 struct WorldCoor *wcs; /* World coordinate system structure */
1328 double *cra; /* Right ascension of image center (deg) (returned) */
1329 double *cdec; /* Declination of image center (deg) (returned) */
1330 double *width; /* Width in degrees (returned) */
1331 double *height; /* Height in degrees (returned) */
1332
1333 {
1334 double xpix, ypix, xpos1, xpos2, ypos1, ypos2, xcpix, ycpix;
1335 double xcent, ycent;
1336
1337 /* Find right ascension and declination of coordinates */
1338 if (iswcs(wcs)) {
1339 xcpix = (0.5 * wcs->nxpix) + 0.5;
1340 ycpix = (0.5 * wcs->nypix) + 0.5;
1341 (void) pix2wcs (wcs,xcpix,ycpix,&xcent, &ycent);
1342 *cra = xcent;
1343 *cdec = ycent;
1344
1345 /* Compute image width in degrees */
1346 xpix = 0.500001;
1347 (void) pix2wcs (wcs,xpix,ycpix,&xpos1,&ypos1);
1348 xpix = wcs->nxpix + 0.499999;
1349 (void) pix2wcs (wcs,xpix,ycpix,&xpos2,&ypos2);
1350 if (strncmp (wcs->ptype,"LIN",3) &&
1351 strncmp (wcs->ptype,"PIX",3)) {
1352 *width = wcsdist (xpos1,ypos1,xpos2,ypos2);
1353 }
1354 else
1355 *width = sqrt (((ypos2-ypos1) * (ypos2-ypos1)) +
1356 ((xpos2-xpos1) * (xpos2-xpos1)));
1357
1358 /* Compute image height in degrees */
1359 ypix = 0.5;
1360 (void) pix2wcs (wcs,xcpix,ypix,&xpos1,&ypos1);
1361 ypix = wcs->nypix + 0.5;
1362 (void) pix2wcs (wcs,xcpix,ypix,&xpos2,&ypos2);
1363 if (strncmp (wcs->ptype,"LIN",3) &&
1364 strncmp (wcs->ptype,"PIX",3))
1365 *height = wcsdist (xpos1,ypos1,xpos2,ypos2);
1366 else
1367 *height = sqrt (((ypos2-ypos1) * (ypos2-ypos1)) +
1368 ((xpos2-xpos1) * (xpos2-xpos1)));
1369 }
1370
1371 else {
1372 *cra = 0.0;
1373 *cdec = 0.0;
1374 *width = 0.0;
1375 *height = 0.0;
1376 }
1377
1378 return;
1379 }
1380
1381
1382 /* Return minimum and maximum RA and Dec of image in degrees */
1383
1384 void
wcsrange(wcs,ra1,ra2,dec1,dec2)1385 wcsrange (wcs, ra1, ra2, dec1, dec2)
1386
1387 struct WorldCoor *wcs; /* World coordinate system structure */
1388 double *ra1; /* Minimum right ascension of image (deg) (returned) */
1389 double *ra2; /* Maximum right ascension of image (deg) (returned) */
1390 double *dec1; /* Minimum declination of image (deg) (returned) */
1391 double *dec2; /* Maximum declination of image (deg) (returned) */
1392
1393 {
1394 double xpos1, xpos2, xpos3, xpos4, ypos1, ypos2, ypos3, ypos4, temp;
1395
1396 if (iswcs(wcs)) {
1397
1398 /* Compute image corner coordinates in degrees */
1399 (void) pix2wcs (wcs,1.0,1.0,&xpos1,&ypos1);
1400 (void) pix2wcs (wcs,1.0,wcs->nypix,&xpos2,&ypos2);
1401 (void) pix2wcs (wcs,wcs->nxpix,1.0,&xpos3,&ypos3);
1402 (void) pix2wcs (wcs,wcs->nxpix,wcs->nypix,&xpos4,&ypos4);
1403
1404 /* Find minimum right ascension or longitude */
1405 *ra1 = xpos1;
1406 if (xpos2 < *ra1) *ra1 = xpos2;
1407 if (xpos3 < *ra1) *ra1 = xpos3;
1408 if (xpos4 < *ra1) *ra1 = xpos4;
1409
1410 /* Find maximum right ascension or longitude */
1411 *ra2 = xpos1;
1412 if (xpos2 > *ra2) *ra2 = xpos2;
1413 if (xpos3 > *ra2) *ra2 = xpos3;
1414 if (xpos4 > *ra2) *ra2 = xpos4;
1415
1416 if (wcs->syswcs != WCS_LINEAR && wcs->syswcs != WCS_XY) {
1417 if (*ra2 - *ra1 > 180.0) {
1418 temp = *ra1;
1419 *ra1 = *ra2;
1420 *ra2 = temp;
1421 }
1422 }
1423
1424 /* Find minimum declination or latitude */
1425 *dec1 = ypos1;
1426 if (ypos2 < *dec1) *dec1 = ypos2;
1427 if (ypos3 < *dec1) *dec1 = ypos3;
1428 if (ypos4 < *dec1) *dec1 = ypos4;
1429
1430 /* Find maximum declination or latitude */
1431 *dec2 = ypos1;
1432 if (ypos2 > *dec2) *dec2 = ypos2;
1433 if (ypos3 > *dec2) *dec2 = ypos3;
1434 if (ypos4 > *dec2) *dec2 = ypos4;
1435 }
1436
1437 else {
1438 *ra1 = 0.0;
1439 *ra2 = 0.0;
1440 *dec1 = 0.0;
1441 *dec2 = 0.0;
1442 }
1443
1444 return;
1445 }
1446
1447
1448 /* Compute distance in degrees between two sky coordinates */
1449
1450 double
wcsdist(x1,y1,x2,y2)1451 wcsdist (x1,y1,x2,y2)
1452
1453 double x1,y1; /* (RA,Dec) or (Long,Lat) in degrees */
1454 double x2,y2; /* (RA,Dec) or (Long,Lat) in degrees */
1455
1456 {
1457 double r, diffi;
1458 double pos1[3], pos2[3], w, diff;
1459 int i;
1460
1461 /* Convert two vectors to direction cosines */
1462 r = 1.0;
1463 d2v3 (x1, y1, r, pos1);
1464 d2v3 (x2, y2, r, pos2);
1465
1466 /* Modulus squared of half the difference vector */
1467 w = 0.0;
1468 for (i = 0; i < 3; i++) {
1469 diffi = pos1[i] - pos2[i];
1470 w = w + (diffi * diffi);
1471 }
1472 w = w / 4.0;
1473 if (w > 1.0) w = 1.0;
1474
1475 /* Angle beween the vectors */
1476 diff = 2.0 * atan2 (sqrt (w), sqrt (1.0 - w));
1477 diff = raddeg (diff);
1478 return (diff);
1479 }
1480
1481
1482
1483 /* Compute distance in degrees between two sky coordinates */
1484
1485 double
wcsdist1(x1,y1,x2,y2)1486 wcsdist1 (x1,y1,x2,y2)
1487
1488 double x1,y1; /* (RA,Dec) or (Long,Lat) in degrees */
1489 double x2,y2; /* (RA,Dec) or (Long,Lat) in degrees */
1490
1491 {
1492 double d1, d2, r;
1493 double pos1[3], pos2[3], w, diff;
1494 int i;
1495
1496 /* Convert two vectors to direction cosines */
1497 r = 1.0;
1498 d2v3 (x1, y1, r, pos1);
1499 d2v3 (x2, y2, r, pos2);
1500
1501 w = 0.0;
1502 d1 = 0.0;
1503 d2 = 0.0;
1504 for (i = 0; i < 3; i++) {
1505 w = w + (pos1[i] * pos2[i]);
1506 d1 = d1 + (pos1[i] * pos1[i]);
1507 d2 = d2 + (pos2[i] * pos2[i]);
1508 }
1509 diff = acosdeg (w / (sqrt (d1) * sqrt (d2)));
1510 return (diff);
1511 }
1512
1513
1514 /* Compute distance in degrees between two sky coordinates away from pole */
1515
1516 double
wcsdiff(x1,y1,x2,y2)1517 wcsdiff (x1,y1,x2,y2)
1518
1519 double x1,y1; /* (RA,Dec) or (Long,Lat) in degrees */
1520 double x2,y2; /* (RA,Dec) or (Long,Lat) in degrees */
1521
1522 {
1523 double xdiff, ydiff, ycos, diff;
1524
1525 ycos = cos (degrad ((y2 + y1) / 2.0));
1526 xdiff = x2 - x1;
1527 if (xdiff > 180.0)
1528 xdiff = xdiff - 360.0;
1529 if (xdiff < -180.0)
1530 xdiff = xdiff + 360.0;
1531 xdiff = xdiff / ycos;
1532 ydiff = (y2 - y1);
1533 diff = sqrt ((xdiff * xdiff) + (ydiff * ydiff));
1534 return (diff);
1535 }
1536
1537
1538 /* Initialize catalog search command set by -wcscom */
1539
1540 void
wcscominit(wcs,i,command)1541 wcscominit (wcs, i, command)
1542
1543 struct WorldCoor *wcs; /* World coordinate system structure */
1544 int i; /* Number of command (0-9) to initialize */
1545 char *command; /* command with %s where coordinates will go */
1546
1547 {
1548 int lcom,icom;
1549
1550 if (iswcs(wcs)) {
1551 lcom = strlen (command);
1552 if (lcom > 0) {
1553 if (wcs->command_format[i] != NULL)
1554 free (wcs->command_format[i]);
1555 wcs->command_format[i] = (char *) calloc (lcom+2, 1);
1556 if (wcs->command_format[i] == NULL)
1557 return;
1558 for (icom = 0; icom < lcom; icom++) {
1559 if (command[icom] == '_')
1560 wcs->command_format[i][icom] = ' ';
1561 else
1562 wcs->command_format[i][icom] = command[icom];
1563 }
1564 wcs->command_format[i][lcom] = 0;
1565 }
1566 }
1567 return;
1568 }
1569
1570
1571 /* Execute Unix command with world coordinates (from x,y) and/or filename */
1572
1573 void
wcscom(wcs,i,filename,xfile,yfile,wcstring)1574 wcscom ( wcs, i, filename, xfile, yfile, wcstring )
1575
1576 struct WorldCoor *wcs; /* World coordinate system structure */
1577 int i; /* Number of command (0-9) to execute */
1578 char *filename; /* Image file name */
1579 double xfile,yfile; /* Image pixel coordinates for WCS command */
1580 char *wcstring; /* WCS String from pix2wcst() */
1581 {
1582 char command[120];
1583 char comform[120];
1584 char xystring[32];
1585 char *fileform, *posform, *imform;
1586 int ier;
1587
1588 if (nowcs (wcs)) {
1589 (void)fprintf(stderr,"WCSCOM: no WCS\n");
1590 return;
1591 }
1592
1593 if (wcs->command_format[i] != NULL)
1594 strcpy (comform, wcs->command_format[i]);
1595 else
1596 strcpy (comform, "sgsc -ah %s");
1597
1598 if (comform[0] > 0) {
1599
1600 /* Create and execute search command */
1601 fileform = strsrch (comform,"%f");
1602 imform = strsrch (comform,"%x");
1603 posform = strsrch (comform,"%s");
1604 if (imform != NULL) {
1605 *(imform+1) = 's';
1606 (void)sprintf (xystring, "%.2f %.2f", xfile, yfile);
1607 if (fileform != NULL) {
1608 *(fileform+1) = 's';
1609 if (posform == NULL) {
1610 if (imform < fileform)
1611 (void)sprintf(command, comform, xystring, filename);
1612 else
1613 (void)sprintf(command, comform, filename, xystring);
1614 }
1615 else if (fileform < posform) {
1616 if (imform < fileform)
1617 (void)sprintf(command, comform, xystring, filename,
1618 wcstring);
1619 else if (imform < posform)
1620 (void)sprintf(command, comform, filename, xystring,
1621 wcstring);
1622 else
1623 (void)sprintf(command, comform, filename, wcstring,
1624 xystring);
1625 }
1626 else
1627 if (imform < posform)
1628 (void)sprintf(command, comform, xystring, wcstring,
1629 filename);
1630 else if (imform < fileform)
1631 (void)sprintf(command, comform, wcstring, xystring,
1632 filename);
1633 else
1634 (void)sprintf(command, comform, wcstring, filename,
1635 xystring);
1636 }
1637 else if (posform == NULL)
1638 (void)sprintf(command, comform, xystring);
1639 else if (imform < posform)
1640 (void)sprintf(command, comform, xystring, wcstring);
1641 else
1642 (void)sprintf(command, comform, wcstring, xystring);
1643 }
1644 else if (fileform != NULL) {
1645 *(fileform+1) = 's';
1646 if (posform == NULL)
1647 (void)sprintf(command, comform, filename);
1648 else if (fileform < posform)
1649 (void)sprintf(command, comform, filename, wcstring);
1650 else
1651 (void)sprintf(command, comform, wcstring, filename);
1652 }
1653 else
1654 (void)sprintf(command, comform, wcstring);
1655 ier = system (command);
1656 if (ier)
1657 (void)fprintf(stderr,"WCSCOM: %s failed %d\n",command,ier);
1658 }
1659 return;
1660 }
1661
1662 /* Initialize WCS output coordinate system for use by PIX2WCS() */
1663
1664 void
wcsoutinit(wcs,coorsys)1665 wcsoutinit (wcs, coorsys)
1666
1667 struct WorldCoor *wcs; /* World coordinate system structure */
1668 char *coorsys; /* Input world coordinate system:
1669 FK4, FK5, B1950, J2000, GALACTIC, ECLIPTIC
1670 fk4, fk5, b1950, j2000, galactic, ecliptic */
1671 {
1672 int sysout, i;
1673
1674 if (nowcs (wcs))
1675 return;
1676
1677 /* If argument is null, set to image system and equinox */
1678 if (coorsys == NULL || strlen (coorsys) < 1 ||
1679 !strcmp(coorsys,"IMSYS") || !strcmp(coorsys,"imsys")) {
1680 sysout = wcs->syswcs;
1681 strcpy (wcs->radecout, wcs->radecsys);
1682 wcs->eqout = wcs->equinox;
1683 if (sysout == WCS_B1950) {
1684 if (wcs->eqout != 1950.0) {
1685 wcs->radecout[0] = 'B';
1686 sprintf (wcs->radecout+1,"%.4f", wcs->equinox);
1687 i = strlen(wcs->radecout) - 1;
1688 if (wcs->radecout[i] == '0')
1689 wcs->radecout[i] = (char)0;
1690 i = strlen(wcs->radecout) - 1;
1691 if (wcs->radecout[i] == '0')
1692 wcs->radecout[i] = (char)0;
1693 i = strlen(wcs->radecout) - 1;
1694 if (wcs->radecout[i] == '0')
1695 wcs->radecout[i] = (char)0;
1696 }
1697 else
1698 strcpy (wcs->radecout, "B1950");
1699 }
1700 else if (sysout == WCS_J2000) {
1701 if (wcs->eqout != 2000.0) {
1702 wcs->radecout[0] = 'J';
1703 sprintf (wcs->radecout+1,"%.4f", wcs->equinox);
1704 i = strlen(wcs->radecout) - 1;
1705 if (wcs->radecout[i] == '0')
1706 wcs->radecout[i] = (char)0;
1707 i = strlen(wcs->radecout) - 1;
1708 if (wcs->radecout[i] == '0')
1709 wcs->radecout[i] = (char)0;
1710 i = strlen(wcs->radecout) - 1;
1711 if (wcs->radecout[i] == '0')
1712 wcs->radecout[i] = (char)0;
1713 }
1714 else
1715 strcpy (wcs->radecout, "J2000");
1716 }
1717 }
1718
1719 /* Ignore new coordinate system if it is not supported */
1720 else {
1721 if ((sysout = wcscsys (coorsys)) < 0)
1722 return;
1723
1724 /* Do not try to convert linear or alt-az coordinates */
1725 if (sysout != wcs->syswcs &&
1726 (wcs->syswcs == WCS_LINEAR || wcs->syswcs == WCS_ALTAZ))
1727 return;
1728
1729 strcpy (wcs->radecout, coorsys);
1730 wcs->eqout = wcsceq (coorsys);
1731 }
1732
1733 wcs->sysout = sysout;
1734 if (wcs->wcson) {
1735
1736 /* Set output in degrees flag and number of decimal places */
1737 if (wcs->sysout == WCS_GALACTIC || wcs->sysout == WCS_ECLIPTIC ||
1738 wcs->sysout == WCS_PLANET) {
1739 wcs->degout = 1;
1740 wcs->ndec = 5;
1741 }
1742 else if (wcs->sysout == WCS_ALTAZ) {
1743 wcs->degout = 1;
1744 wcs->ndec = 5;
1745 }
1746 else if (wcs->sysout == WCS_NPOLE || wcs->sysout == WCS_SPA) {
1747 wcs->degout = 1;
1748 wcs->ndec = 5;
1749 }
1750 else {
1751 wcs->degout = 0;
1752 wcs->ndec = 3;
1753 }
1754 }
1755 return;
1756 }
1757
1758
1759 /* Return current value of WCS output coordinate system set by -wcsout */
1760 char *
getwcsout(wcs)1761 getwcsout(wcs)
1762 struct WorldCoor *wcs; /* World coordinate system structure */
1763 {
1764 if (nowcs (wcs))
1765 return (NULL);
1766 else
1767 return(wcs->radecout);
1768 }
1769
1770
1771 /* Initialize WCS input coordinate system for use by WCS2PIX() */
1772
1773 void
wcsininit(wcs,coorsys)1774 wcsininit (wcs, coorsys)
1775
1776 struct WorldCoor *wcs; /* World coordinate system structure */
1777 char *coorsys; /* Input world coordinate system:
1778 FK4, FK5, B1950, J2000, GALACTIC, ECLIPTIC
1779 fk4, fk5, b1950, j2000, galactic, ecliptic */
1780 {
1781 int sysin, i;
1782
1783 if (nowcs (wcs))
1784 return;
1785
1786 /* If argument is null, set to image system and equinox */
1787 if (coorsys == NULL || strlen (coorsys) < 1) {
1788 wcs->sysin = wcs->syswcs;
1789 strcpy (wcs->radecin, wcs->radecsys);
1790 wcs->eqin = wcs->equinox;
1791 if (wcs->sysin == WCS_B1950) {
1792 if (wcs->eqin != 1950.0) {
1793 wcs->radecin[0] = 'B';
1794 sprintf (wcs->radecin+1,"%.4f", wcs->equinox);
1795 i = strlen(wcs->radecin) - 1;
1796 if (wcs->radecin[i] == '0')
1797 wcs->radecin[i] = (char)0;
1798 i = strlen(wcs->radecin) - 1;
1799 if (wcs->radecin[i] == '0')
1800 wcs->radecin[i] = (char)0;
1801 i = strlen(wcs->radecin) - 1;
1802 if (wcs->radecin[i] == '0')
1803 wcs->radecin[i] = (char)0;
1804 }
1805 else
1806 strcpy (wcs->radecin, "B1950");
1807 }
1808 else if (wcs->sysin == WCS_J2000) {
1809 if (wcs->eqin != 2000.0) {
1810 wcs->radecin[0] = 'J';
1811 sprintf (wcs->radecin+1,"%.4f", wcs->equinox);
1812 i = strlen(wcs->radecin) - 1;
1813 if (wcs->radecin[i] == '0')
1814 wcs->radecin[i] = (char)0;
1815 i = strlen(wcs->radecin) - 1;
1816 if (wcs->radecin[i] == '0')
1817 wcs->radecin[i] = (char)0;
1818 i = strlen(wcs->radecin) - 1;
1819 if (wcs->radecin[i] == '0')
1820 wcs->radecin[i] = (char)0;
1821 }
1822 else
1823 strcpy (wcs->radecin, "J2000");
1824 }
1825 }
1826
1827 /* Ignore new coordinate system if it is not supported */
1828 if ((sysin = wcscsys (coorsys)) < 0)
1829 return;
1830
1831 wcs->sysin = sysin;
1832 wcs->eqin = wcsceq (coorsys);
1833 strcpy (wcs->radecin, coorsys);
1834 return;
1835 }
1836
1837
1838 /* Return current value of WCS input coordinate system set by wcsininit */
1839 char *
getwcsin(wcs)1840 getwcsin (wcs)
1841 struct WorldCoor *wcs; /* World coordinate system structure */
1842 {
1843 if (nowcs (wcs))
1844 return (NULL);
1845 else
1846 return (wcs->radecin);
1847 }
1848
1849
1850 /* Set WCS output in degrees or hh:mm:ss dd:mm:ss, returning old flag value */
1851 int
setwcsdeg(wcs,new)1852 setwcsdeg(wcs, new)
1853 struct WorldCoor *wcs; /* World coordinate system structure */
1854 int new; /* 1: degrees, 0: h:m:s d:m:s */
1855 {
1856 int old;
1857
1858 if (nowcs (wcs))
1859 return (0);
1860 old = wcs->degout;
1861 wcs->degout = new;
1862 if (new == 1 && old == 0 && wcs->ndec == 3)
1863 wcs->ndec = 6;
1864 if (new == 0 && old == 1 && wcs->ndec == 5)
1865 wcs->ndec = 3;
1866 return(old);
1867 }
1868
1869
1870 /* Set number of decimal places in pix2wcst output string */
1871 int
wcsndec(wcs,ndec)1872 wcsndec (wcs, ndec)
1873 struct WorldCoor *wcs; /* World coordinate system structure */
1874 int ndec; /* Number of decimal places in output string */
1875 /* If < 0, return current unchanged value */
1876 {
1877 if (nowcs (wcs))
1878 return (0);
1879 else if (ndec >= 0)
1880 wcs->ndec = ndec;
1881 return (wcs->ndec);
1882 }
1883
1884
1885
1886 /* Return current value of coordinate system */
1887 char *
getradecsys(wcs)1888 getradecsys(wcs)
1889 struct WorldCoor *wcs; /* World coordinate system structure */
1890 {
1891 if (nowcs (wcs))
1892 return (NULL);
1893 else
1894 return (wcs->radecsys);
1895 }
1896
1897
1898 /* Set output string mode for LINEAR coordinates */
1899
1900 void
setwcslin(wcs,mode)1901 setwcslin (wcs, mode)
1902 struct WorldCoor *wcs; /* World coordinate system structure */
1903 int mode; /* mode = 0: x y linear
1904 mode = 1: x units x units
1905 mode = 2: x y linear units */
1906 {
1907 if (iswcs (wcs))
1908 wcs->linmode = mode;
1909 return;
1910 }
1911
1912
1913 /* Convert pixel coordinates to World Coordinate string */
1914
1915 int
pix2wcst(wcs,xpix,ypix,wcstring,lstr)1916 pix2wcst (wcs, xpix, ypix, wcstring, lstr)
1917
1918 struct WorldCoor *wcs; /* World coordinate system structure */
1919 double xpix,ypix; /* Image coordinates in pixels */
1920 char *wcstring; /* World coordinate string (returned) */
1921 int lstr; /* Length of world coordinate string (returned) */
1922 {
1923 double xpos,ypos;
1924 char rastr[32], decstr[32];
1925 int minlength, lunits, lstring;
1926
1927 if (nowcs (wcs)) {
1928 if (lstr > 0)
1929 wcstring[0] = 0;
1930 return(0);
1931 }
1932
1933 pix2wcs (wcs,xpix,ypix,&xpos,&ypos);
1934
1935 /* If point is off scale, set string accordingly */
1936 if (wcs->offscl) {
1937 (void)sprintf (wcstring,"Off map");
1938 return (1);
1939 }
1940
1941 /* Print coordinates in degrees */
1942 else if (wcs->degout == 1) {
1943 minlength = 9 + (2 * wcs->ndec);
1944 if (lstr > minlength) {
1945 deg2str (rastr, 32, xpos, wcs->ndec);
1946 deg2str (decstr, 32, ypos, wcs->ndec);
1947 if (wcs->tabsys)
1948 (void)sprintf (wcstring,"%s %s", rastr, decstr);
1949 else
1950 (void)sprintf (wcstring,"%s %s", rastr, decstr);
1951 lstr = lstr - minlength;
1952 }
1953 else {
1954 if (wcs->tabsys)
1955 strncpy (wcstring,"********* **********",lstr);
1956 else
1957 strncpy (wcstring,"*******************",lstr);
1958 lstr = 0;
1959 }
1960 }
1961
1962 /* print coordinates in sexagesimal notation */
1963 else if (wcs->degout == 0) {
1964 minlength = 18 + (2 * wcs->ndec);
1965 if (lstr > minlength) {
1966 if (wcs->sysout == WCS_J2000 || wcs->sysout == WCS_B1950) {
1967 ra2str (rastr, 32, xpos, wcs->ndec);
1968 dec2str (decstr, 32, ypos, wcs->ndec-1);
1969 }
1970 else {
1971 dec2str (rastr, 32, xpos, wcs->ndec);
1972 dec2str (decstr, 32, ypos, wcs->ndec);
1973 }
1974 if (wcs->tabsys) {
1975 (void)sprintf (wcstring,"%s %s", rastr, decstr);
1976 }
1977 else {
1978 (void)sprintf (wcstring,"%s %s", rastr, decstr);
1979 }
1980 lstr = lstr - minlength;
1981 }
1982 else {
1983 if (wcs->tabsys) {
1984 strncpy (wcstring,"************* *************",lstr);
1985 }
1986 else {
1987 strncpy (wcstring,"**************************",lstr);
1988 }
1989 lstr = 0;
1990 }
1991 }
1992
1993 /* Label galactic coordinates */
1994 if (wcs->sysout == WCS_GALACTIC) {
1995 if (lstr > 9 && wcs->printsys) {
1996 if (wcs->tabsys)
1997 strcat (wcstring," galactic");
1998 else
1999 strcat (wcstring," galactic");
2000 }
2001 }
2002
2003 /* Label ecliptic coordinates */
2004 else if (wcs->sysout == WCS_ECLIPTIC) {
2005 if (lstr > 9 && wcs->printsys) {
2006 if (wcs->tabsys)
2007 strcat (wcstring," ecliptic");
2008 else
2009 strcat (wcstring," ecliptic");
2010 }
2011 }
2012
2013 /* Label planet coordinates */
2014 else if (wcs->sysout == WCS_PLANET) {
2015 if (lstr > 9 && wcs->printsys) {
2016 if (wcs->tabsys)
2017 strcat (wcstring," planet");
2018 else
2019 strcat (wcstring," planet");
2020 }
2021 }
2022
2023 /* Label alt-az coordinates */
2024 else if (wcs->sysout == WCS_ALTAZ) {
2025 if (lstr > 7 && wcs->printsys) {
2026 if (wcs->tabsys)
2027 strcat (wcstring," alt-az");
2028 else
2029 strcat (wcstring," alt-az");
2030 }
2031 }
2032
2033 /* Label north pole angle coordinates */
2034 else if (wcs->sysout == WCS_NPOLE) {
2035 if (lstr > 7 && wcs->printsys) {
2036 if (wcs->tabsys)
2037 strcat (wcstring," long-npa");
2038 else
2039 strcat (wcstring," long-npa");
2040 }
2041 }
2042
2043 /* Label south pole angle coordinates */
2044 else if (wcs->sysout == WCS_SPA) {
2045 if (lstr > 7 && wcs->printsys) {
2046 if (wcs->tabsys)
2047 strcat (wcstring," long-spa");
2048 else
2049 strcat (wcstring," long-spa");
2050 }
2051 }
2052
2053 /* Label equatorial coordinates */
2054 else if (wcs->sysout==WCS_B1950 || wcs->sysout==WCS_J2000) {
2055 if (lstr > (int) strlen(wcs->radecout)+1 && wcs->printsys) {
2056 if (wcs->tabsys)
2057 strcat (wcstring," ");
2058 else
2059 strcat (wcstring," ");
2060 strcat (wcstring, wcs->radecout);
2061 }
2062 }
2063
2064 /* Output linear coordinates */
2065 else {
2066 num2str (rastr, xpos, 0, wcs->ndec);
2067 num2str (decstr, ypos, 0, wcs->ndec);
2068 lstring = strlen (rastr) + strlen (decstr) + 1;
2069 lunits = strlen (wcs->units[0]) + strlen (wcs->units[1]) + 2;
2070 if (wcs->syswcs == WCS_LINEAR && wcs->linmode == 1) {
2071 if (lstr > lstring + lunits) {
2072 if (strlen (wcs->units[0]) > 0) {
2073 strcat (rastr, " ");
2074 strcat (rastr, wcs->units[0]);
2075 }
2076 if (strlen (wcs->units[1]) > 0) {
2077 strcat (decstr, " ");
2078 strcat (decstr, wcs->units[1]);
2079 }
2080 lstring = lstring + lunits;
2081 }
2082 }
2083 if (lstr > lstring) {
2084 if (wcs->tabsys)
2085 (void)sprintf (wcstring,"%s %s", rastr, decstr);
2086 else
2087 (void)sprintf (wcstring,"%s %s", rastr, decstr);
2088 }
2089 else {
2090 if (wcs->tabsys)
2091 strncpy (wcstring,"********** *********",lstr);
2092 else
2093 strncpy (wcstring,"*******************",lstr);
2094 }
2095 if (wcs->syswcs == WCS_LINEAR && wcs->linmode != 1 &&
2096 lstr > lstring + 7)
2097 strcat (wcstring, " linear");
2098 if (wcs->syswcs == WCS_LINEAR && wcs->linmode == 2 &&
2099 lstr > lstring + lunits + 7) {
2100 if (strlen (wcs->units[0]) > 0) {
2101 strcat (wcstring, " ");
2102 strcat (wcstring, wcs->units[0]);
2103 }
2104 if (strlen (wcs->units[1]) > 0) {
2105 strcat (wcstring, " ");
2106 strcat (wcstring, wcs->units[1]);
2107 }
2108
2109 }
2110 }
2111 return (1);
2112 }
2113
2114
2115 /* Convert pixel coordinates to World Coordinates */
2116
2117 void
pix2wcs(wcs,xpix,ypix,xpos,ypos)2118 pix2wcs (wcs,xpix,ypix,xpos,ypos)
2119
2120 struct WorldCoor *wcs; /* World coordinate system structure */
2121 double xpix,ypix; /* x and y image coordinates in pixels */
2122 double *xpos,*ypos; /* RA and Dec in degrees (returned) */
2123 {
2124 double xpi, ypi, xp, yp;
2125 double eqin, eqout;
2126 int wcspos();
2127
2128 if (nowcs (wcs))
2129 return;
2130 wcs->xpix = xpix;
2131 wcs->ypix = ypix;
2132 wcs->zpix = zpix;
2133 wcs->offscl = 0;
2134
2135 /* If this WCS is converted from another WCS rather than pixels, convert now */
2136 if (wcs->wcs != NULL) {
2137 pix2wcs (wcs->wcs, xpix, ypix, &xpi, &ypi);
2138 }
2139 else {
2140 pix2foc (wcs, xpix, ypix, &xpi, &ypi);
2141 }
2142
2143 /* Convert image coordinates to sky coordinates */
2144
2145 /* Use Digitized Sky Survey plate fit */
2146 if (wcs->prjcode == WCS_DSS) {
2147 if (dsspos (xpi, ypi, wcs, &xp, &yp))
2148 wcs->offscl = 1;
2149 }
2150
2151 /* Use SAO plate fit */
2152 else if (wcs->prjcode == WCS_PLT) {
2153 if (platepos (xpi, ypi, wcs, &xp, &yp))
2154 wcs->offscl = 1;
2155 }
2156
2157 /* Use NOAO IRAF corrected plane tangent projection */
2158 else if (wcs->prjcode == WCS_TNX) {
2159 if (tnxpos (xpi, ypi, wcs, &xp, &yp))
2160 wcs->offscl = 1;
2161 }
2162
2163 /* Use NOAO IRAF corrected zenithal projection */
2164 else if (wcs->prjcode == WCS_ZPX) {
2165 if (zpxpos (xpi, ypi, wcs, &xp, &yp))
2166 wcs->offscl = 1;
2167 }
2168
2169 /* Use Classic AIPS projections */
2170 else if (wcs->wcsproj == WCS_OLD || wcs->prjcode <= 0) {
2171 if (worldpos (xpi, ypi, wcs, &xp, &yp))
2172 wcs->offscl = 1;
2173 }
2174
2175 /* Use Mark Calabretta's WCSLIB projections */
2176 else if (wcspos (xpi, ypi, wcs, &xp, &yp))
2177 wcs->offscl = 1;
2178
2179
2180 /* Do not change coordinates if offscale */
2181 if (wcs->offscl) {
2182 *xpos = 0.0;
2183 *ypos = 0.0;
2184 }
2185 else {
2186
2187 /* Convert coordinates to output system, if not LINEAR */
2188 if (wcs->prjcode > 0) {
2189
2190 /* Convert coordinates to desired output system */
2191 eqin = wcs->equinox;
2192 eqout = wcs->eqout;
2193 wcscon (wcs->syswcs,wcs->sysout,eqin,eqout,&xp,&yp,wcs->epoch);
2194 }
2195 if (wcs->latbase == 90)
2196 yp = 90.0 - yp;
2197 else if (wcs->latbase == -90)
2198 yp = yp - 90.0;
2199 wcs->xpos = xp;
2200 wcs->ypos = yp;
2201 *xpos = xp;
2202 *ypos = yp;
2203 }
2204
2205 /* Keep RA/longitude within range if spherical coordinate output
2206 (Not LINEAR or XY) */
2207 if (wcs->sysout > 0 && wcs->sysout != 6 && wcs->sysout != 10) {
2208 if (*xpos < 0.0)
2209 *xpos = *xpos + 360.0;
2210 else if (*xpos > 360.0)
2211 *xpos = *xpos - 360.0;
2212 }
2213
2214 return;
2215 }
2216
2217
2218 /* Convert World Coordinates to pixel coordinates */
2219
2220 void
wcs2pix(wcs,xpos,ypos,xpix,ypix,offscl)2221 wcs2pix (wcs, xpos, ypos, xpix, ypix, offscl)
2222
2223 struct WorldCoor *wcs; /* World coordinate system structure */
2224 double xpos,ypos; /* World coordinates in degrees */
2225 double *xpix,*ypix; /* Image coordinates in pixels */
2226 int *offscl; /* 0 if within bounds, else off scale */
2227 {
2228 wcsc2pix (wcs, xpos, ypos, wcs->radecin, xpix, ypix, offscl);
2229 return;
2230 }
2231
2232 /* Convert World Coordinates to pixel coordinates */
2233
2234 void
wcsc2pix(wcs,xpos,ypos,coorsys,xpix,ypix,offscl)2235 wcsc2pix (wcs, xpos, ypos, coorsys, xpix, ypix, offscl)
2236
2237 struct WorldCoor *wcs; /* World coordinate system structure */
2238 double xpos,ypos; /* World coordinates in degrees */
2239 char *coorsys; /* Input world coordinate system:
2240 FK4, FK5, B1950, J2000, GALACTIC, ECLIPTIC
2241 fk4, fk5, b1950, j2000, galactic, ecliptic
2242 * If NULL, use image WCS */
2243 double *xpix,*ypix; /* Image coordinates in pixels */
2244 int *offscl; /* 0 if within bounds, else off scale */
2245 {
2246 double xp, yp, xpi, ypi;
2247 double eqin, eqout;
2248 int sysin;
2249 int wcspix();
2250
2251 if (nowcs (wcs))
2252 return;
2253
2254 *offscl = 0;
2255 xp = xpos;
2256 yp = ypos;
2257 if (wcs->latbase == 90)
2258 yp = 90.0 - yp;
2259 else if (wcs->latbase == -90)
2260 yp = yp - 90.0;
2261 if (coorsys == NULL) {
2262 sysin = wcs->syswcs;
2263 eqin = wcs->equinox;
2264 }
2265 else {
2266 sysin = wcscsys (coorsys);
2267 eqin = wcsceq (coorsys);
2268 }
2269 wcs->zpix = 1.0;
2270
2271 /* Convert coordinates to same system as image */
2272 if (sysin > 0 && sysin != 6 && sysin != 10) {
2273 eqout = wcs->equinox;
2274 wcscon (sysin, wcs->syswcs, eqin, eqout, &xp, &yp, wcs->epoch);
2275 }
2276
2277 /* Convert sky coordinates to image coordinates */
2278
2279 /* Use Digitized Sky Survey plate fit */
2280 if (wcs->prjcode == WCS_DSS) {
2281 if (dsspix (xp, yp, wcs, &xpi, &ypi))
2282 *offscl = 1;
2283 }
2284
2285 /* Use SAO polynomial plate fit */
2286 else if (wcs->prjcode == WCS_PLT) {
2287 if (platepix (xp, yp, wcs, &xpi, &ypi))
2288 *offscl = 1;
2289 }
2290
2291 /* Use NOAO IRAF corrected plane tangent projection */
2292 else if (wcs->prjcode == WCS_TNX) {
2293 if (tnxpix (xp, yp, wcs, &xpi, &ypi))
2294 *offscl = 1;
2295 }
2296
2297 /* Use NOAO IRAF corrected zenithal projection */
2298 else if (wcs->prjcode == WCS_ZPX) {
2299 if (zpxpix (xp, yp, wcs, &xpi, &ypi))
2300 *offscl = 1;
2301 }
2302
2303 /* Use Classic AIPS projections */
2304 else if (wcs->wcsproj == WCS_OLD || wcs->prjcode <= 0) {
2305 if (worldpix (xp, yp, wcs, &xpi, &ypi))
2306 *offscl = 1;
2307 }
2308
2309 /* Use Mark Calabretta's WCSLIB projections */
2310 else if (wcspix (xp, yp, wcs, &xpi, &ypi)) {
2311 *offscl = 1;
2312 }
2313
2314 /* If this WCS is converted from another WCS rather than pixels, convert now */
2315 if (wcs->wcs != NULL) {
2316 wcsc2pix (wcs->wcs, xpi, ypi, NULL, xpix, ypix, offscl);
2317 }
2318 else {
2319 foc2pix (wcs, xpi, ypi, xpix, ypix);
2320
2321 /* Set off-scale flag to 2 if off image but within bounds of projection */
2322 if (!*offscl) {
2323 if (*xpix < 0.5 || *ypix < 0.5)
2324 *offscl = 2;
2325 else if (*xpix > wcs->nxpix + 0.5 || *ypix > wcs->nypix + 0.5)
2326 *offscl = 2;
2327 }
2328 }
2329
2330 wcs->offscl = *offscl;
2331 wcs->xpos = xpos;
2332 wcs->ypos = ypos;
2333 wcs->xpix = *xpix;
2334 wcs->ypix = *ypix;
2335
2336 return;
2337 }
2338
2339
2340 int
wcspos(xpix,ypix,wcs,xpos,ypos)2341 wcspos (xpix, ypix, wcs, xpos, ypos)
2342
2343 /* Input: */
2344 double xpix; /* x pixel number (RA or long without rotation) */
2345 double ypix; /* y pixel number (dec or lat without rotation) */
2346 struct WorldCoor *wcs; /* WCS parameter structure */
2347
2348 /* Output: */
2349 double *xpos; /* x (RA) coordinate (deg) */
2350 double *ypos; /* y (dec) coordinate (deg) */
2351 {
2352 int offscl;
2353 int i;
2354 int wcsrev();
2355 double wcscrd[4], imgcrd[4], pixcrd[4];
2356 double phi, theta;
2357
2358 *xpos = 0.0;
2359 *ypos = 0.0;
2360
2361 pixcrd[0] = xpix;
2362 pixcrd[1] = ypix;
2363 if (wcs->prjcode == WCS_CSC || wcs->prjcode == WCS_QSC ||
2364 wcs->prjcode == WCS_TSC)
2365 pixcrd[2] = (double) (izpix + 1);
2366 else
2367 pixcrd[2] = zpix;
2368 pixcrd[3] = 1.0;
2369 for (i = 0; i < 4; i++)
2370 imgcrd[i] = 0.0;
2371 offscl = wcsrev ((void *)&wcs->ctype, &wcs->wcsl, pixcrd, &wcs->lin, imgcrd,
2372 &wcs->prj, &phi, &theta, wcs->crval, &wcs->cel, wcscrd);
2373 if (offscl == 0) {
2374 *xpos = wcscrd[wcs->wcsl.lng];
2375 *ypos = wcscrd[wcs->wcsl.lat];
2376 }
2377
2378 return (offscl);
2379 }
2380
2381 int
wcspix(xpos,ypos,wcs,xpix,ypix)2382 wcspix (xpos, ypos, wcs, xpix, ypix)
2383
2384 /* Input: */
2385 double xpos; /* x (RA) coordinate (deg) */
2386 double ypos; /* y (dec) coordinate (deg) */
2387 struct WorldCoor *wcs; /* WCS parameter structure */
2388
2389 /* Output: */
2390 double *xpix; /* x pixel number (RA or long without rotation) */
2391 double *ypix; /* y pixel number (dec or lat without rotation) */
2392
2393 {
2394 int offscl;
2395 int wcsfwd();
2396 double wcscrd[4], imgcrd[4], pixcrd[4];
2397 double phi, theta;
2398
2399 *xpix = 0.0;
2400 *ypix = 0.0;
2401 if (wcs->wcsl.flag != WCSSET) {
2402 if (wcsset (wcs->lin.naxis, (void *)&wcs->ctype, &wcs->wcsl) )
2403 return (1);
2404 }
2405
2406 /* Set input for WCSLIB subroutines */
2407 wcscrd[0] = 0.0;
2408 wcscrd[1] = 0.0;
2409 wcscrd[2] = 0.0;
2410 wcscrd[3] = 0.0;
2411 wcscrd[wcs->wcsl.lng] = xpos;
2412 wcscrd[wcs->wcsl.lat] = ypos;
2413
2414 /* Initialize output for WCSLIB subroutines */
2415 pixcrd[0] = 0.0;
2416 pixcrd[1] = 0.0;
2417 pixcrd[2] = 1.0;
2418 pixcrd[3] = 1.0;
2419 imgcrd[0] = 0.0;
2420 imgcrd[1] = 0.0;
2421 imgcrd[2] = 1.0;
2422 imgcrd[3] = 1.0;
2423
2424 /* Invoke WCSLIB subroutines for coordinate conversion */
2425 offscl = wcsfwd ((void *)&wcs->ctype, &wcs->wcsl, wcscrd, wcs->crval, &wcs->cel,
2426 &phi, &theta, &wcs->prj, imgcrd, &wcs->lin, pixcrd);
2427
2428 if (!offscl) {
2429 *xpix = pixcrd[0];
2430 *ypix = pixcrd[1];
2431 if (wcs->prjcode == WCS_CSC || wcs->prjcode == WCS_QSC ||
2432 wcs->prjcode == WCS_TSC)
2433 wcs->zpix = pixcrd[2] - 1.0;
2434 else
2435 wcs->zpix = pixcrd[2];
2436 }
2437 return (offscl);
2438 }
2439
2440
2441 /* Set third dimension for cube projections */
2442
2443 int
wcszin(izpix0)2444 wcszin (izpix0)
2445
2446 int izpix0; /* coordinate in third dimension
2447 (if < 0, return current value without changing it */
2448 {
2449 if (izpix0 > -1) {
2450 izpix = izpix0;
2451 zpix = (double) izpix0;
2452 }
2453 return (izpix);
2454 }
2455
2456
2457 /* Return third dimension for cube projections */
2458
2459 int
wcszout(wcs)2460 wcszout (wcs)
2461
2462 struct WorldCoor *wcs; /* WCS parameter structure */
2463 {
2464 return ((int) wcs->zpix);
2465 }
2466
2467 /* Set file name for error messages */
2468 void
setwcsfile(filename)2469 setwcsfile (filename)
2470 char *filename; /* FITS or IRAF file with WCS */
2471 { if (strlen (filename) < 256)
2472 strcpy (wcsfile, filename);
2473 else
2474 strncpy (wcsfile, filename, 255);
2475 return; }
2476
2477 /* Set error message */
2478 void
setwcserr(errmsg)2479 setwcserr (errmsg)
2480 char *errmsg; /* Error mesage < 80 char */
2481 { strcpy (wcserrmsg, errmsg); return; }
2482
2483 /* Print error message */
2484 void
wcserr()2485 wcserr ()
2486 { if (strlen (wcsfile) > 0)
2487 fprintf (stderr, "%s in file %s\n",wcserrmsg, wcsfile);
2488 else
2489 fprintf (stderr, "%s\n",wcserrmsg);
2490 return; }
2491
2492
2493 /* Flag to use AIPS WCS subroutines instead of WCSLIB */
2494 void
setdefwcs(wp)2495 setdefwcs (wp)
2496 int wp;
2497 { wcsproj0 = wp; return; }
2498
2499 int
getdefwcs()2500 getdefwcs ()
2501 { return (wcsproj0); }
2502
2503 /* Save output default coordinate system */
2504 static char wcscoor0[16];
2505
2506 void
savewcscoor(wcscoor)2507 savewcscoor (wcscoor)
2508 char *wcscoor;
2509 { strcpy (wcscoor0, wcscoor); return; }
2510
2511 /* Return preset output default coordinate system */
2512 char *
getwcscoor()2513 getwcscoor ()
2514 { return (wcscoor0); }
2515
2516
2517 /* Save default commands */
2518 static char *wcscom0[10];
2519
2520 void
savewcscom(i,wcscom)2521 savewcscom (i, wcscom)
2522 int i;
2523 char *wcscom;
2524 {
2525 int lcom;
2526 if (i < 0) i = 0;
2527 else if (i > 9) i = 9;
2528 lcom = strlen (wcscom) + 2;
2529 wcscom0[i] = (char *) calloc (lcom, 1);
2530 if (wcscom0[i] != NULL)
2531 strcpy (wcscom0[i], wcscom);
2532 return;
2533 }
2534
2535 void
setwcscom(wcs)2536 setwcscom (wcs)
2537 struct WorldCoor *wcs; /* WCS parameter structure */
2538 {
2539 char envar[16];
2540 int i;
2541 char *str;
2542 if (nowcs(wcs))
2543 return;
2544 for (i = 0; i < 10; i++) {
2545 if (i == 0)
2546 strcpy (envar, "WCS_COMMAND");
2547 else
2548 sprintf (envar, "WCS_COMMAND%d", i);
2549 if (wcscom0[i] != NULL)
2550 wcscominit (wcs, i, wcscom0[i]);
2551 else if ((str = getenv (envar)) != NULL)
2552 wcscominit (wcs, i, str);
2553 else if (i == 1)
2554 wcscominit (wcs, i, "sua2 -ah %s"); /* F1= Search USNO-A2.0 Catalog */
2555 else if (i == 2)
2556 wcscominit (wcs, i, "sgsc -ah %s"); /* F2= Search HST GSC */
2557 else if (i == 3)
2558 wcscominit (wcs, i, "sty2 -ah %s"); /* F3= Search Tycho-2 Catalog */
2559 else if (i == 4)
2560 wcscominit (wcs, i, "sppm -ah %s"); /* F4= Search PPM Catalog */
2561 else if (i == 5)
2562 wcscominit (wcs, i, "ssao -ah %s"); /* F5= Search SAO Catalog */
2563 else
2564 wcs->command_format[i] = NULL;
2565 }
2566 return;
2567 }
2568
2569 char *
getwcscom(i)2570 getwcscom (i)
2571 int i;
2572 { return (wcscom0[i]); }
2573
2574
2575 void
freewcscom(wcs)2576 freewcscom (wcs)
2577 struct WorldCoor *wcs; /* WCS parameter structure */
2578 {
2579 int i;
2580 for (i = 0; i < 10; i++) {
2581 if (wcscom0[i] != NULL) {
2582 free (wcscom0[i]);
2583 wcscom0[i] = NULL;
2584 }
2585 }
2586 if (iswcs (wcs)) {
2587 for (i = 0; i < 10; i++) {
2588 if (wcs->command_format[i] != NULL) {
2589 free (wcs->command_format[i]);
2590 }
2591 }
2592 }
2593 return;
2594 }
2595
2596 int
cpwcs(header,cwcs)2597 cpwcs (header, cwcs)
2598
2599 char **header; /* Pointer to start of FITS header */
2600 char *cwcs; /* Keyword suffix character for output WCS */
2601 {
2602 double tnum;
2603 int dkwd[MAXNKWD];
2604 int i, maxnkwd, ikwd, nleft, lbuff, lhead, nkwd, nbytes;
2605 int nkwdw;
2606 char **kwd;
2607 char *newhead, *oldhead;
2608 char kwdc[16], keyword[16];
2609 char tstr[80];
2610
2611 /* Allocate array of keywords to be transferred */
2612 maxnkwd = MAXNKWD;
2613 kwd = (char **)calloc (maxnkwd, sizeof(char *));
2614 for (ikwd = 0; ikwd < maxnkwd; ikwd++)
2615 kwd[ikwd] = (char *) calloc (16, 1);
2616
2617 /* Make list of all possible keywords to be transferred */
2618 nkwd = 0;
2619 strcpy (kwd[++nkwd], "EPOCH");
2620 dkwd[nkwd] = 1;
2621 strcpy (kwd[++nkwd], "EQUINOX");
2622 dkwd[nkwd] = 1;
2623 strcpy (kwd[++nkwd], "RADECSYS");
2624 dkwd[nkwd] = 0;
2625 strcpy (kwd[++nkwd], "CTYPE1");
2626 dkwd[nkwd] = 0;
2627 strcpy (kwd[++nkwd], "CTYPE2");
2628 dkwd[nkwd] = 0;
2629 strcpy (kwd[++nkwd], "CRVAL1");
2630 dkwd[nkwd] = 1;
2631 strcpy (kwd[++nkwd], "CRVAL2");
2632 dkwd[nkwd] = 1;
2633 strcpy (kwd[++nkwd], "CDELT1");
2634 dkwd[nkwd] = 1;
2635 strcpy (kwd[++nkwd], "CDELT2");
2636 dkwd[nkwd] = 1;
2637 strcpy (kwd[++nkwd], "CRPIX1");
2638 dkwd[nkwd] = 1;
2639 strcpy (kwd[++nkwd], "CRPIX2");
2640 dkwd[nkwd] = 1;
2641 strcpy (kwd[++nkwd], "CROTA1");
2642 dkwd[nkwd] = 1;
2643 strcpy (kwd[++nkwd], "CROTA2");
2644 dkwd[nkwd] = 1;
2645 strcpy (kwd[++nkwd], "CD1_1");
2646 dkwd[nkwd] = 1;
2647 strcpy (kwd[++nkwd], "CD1_2");
2648 dkwd[nkwd] = 1;
2649 strcpy (kwd[++nkwd], "CD2_1");
2650 dkwd[nkwd] = 1;
2651 strcpy (kwd[++nkwd], "CD2_2");
2652 dkwd[nkwd] = 1;
2653 strcpy (kwd[++nkwd], "PC1_1");
2654 dkwd[nkwd] = 1;
2655 strcpy (kwd[++nkwd], "PC1_2");
2656 dkwd[nkwd] = 1;
2657 strcpy (kwd[++nkwd], "PC2_1");
2658 dkwd[nkwd] = 1;
2659 strcpy (kwd[++nkwd], "PC2_2");
2660 dkwd[nkwd] = 1;
2661 strcpy (kwd[++nkwd], "PC001001");
2662 dkwd[nkwd] = 1;
2663 strcpy (kwd[++nkwd], "PC001002");
2664 dkwd[nkwd] = 1;
2665 strcpy (kwd[++nkwd], "PC002001");
2666 dkwd[nkwd] = 1;
2667 strcpy (kwd[++nkwd], "PC002002");
2668 dkwd[nkwd] = 1;
2669 strcpy (kwd[++nkwd], "LATPOLE");
2670 dkwd[nkwd] = 1;
2671 strcpy (kwd[++nkwd], "LONPOLE");
2672 dkwd[nkwd] = 1;
2673 for (i = 1; i < 13; i++) {
2674 sprintf (keyword,"CO1_%d", i);
2675 strcpy (kwd[++nkwd], keyword);
2676 dkwd[nkwd] = 1;
2677 }
2678 for (i = 1; i < 13; i++) {
2679 sprintf (keyword,"CO2_%d", i);
2680 strcpy (kwd[++nkwd], keyword);
2681 dkwd[nkwd] = 1;
2682 }
2683 for (i = 0; i < 10; i++) {
2684 sprintf (keyword,"PROJP%d", i);
2685 strcpy (kwd[++nkwd], keyword);
2686 dkwd[nkwd] = 1;
2687 }
2688 for (i = 0; i < MAXPV; i++) {
2689 sprintf (keyword,"PV1_%d", i);
2690 strcpy (kwd[++nkwd], keyword);
2691 dkwd[nkwd] = 1;
2692 }
2693 for (i = 0; i < MAXPV; i++) {
2694 sprintf (keyword,"PV2_%d", i);
2695 strcpy (kwd[++nkwd], keyword);
2696 dkwd[nkwd] = 1;
2697 }
2698
2699 /* Allocate new header buffer if needed */
2700 lhead = (ksearch (*header, "END") - *header) + 80;
2701 lbuff = gethlength (*header);
2702 nleft = (lbuff - lhead) / 80;
2703 if (nleft < nkwd) {
2704 nbytes = lhead + (nkwd * 80) + 400;
2705 newhead = (char *) calloc (1, nbytes);
2706 strncpy (newhead, *header, lhead);
2707 oldhead = *header;
2708 header = &newhead;
2709 free (oldhead);
2710 }
2711
2712 /* Copy keywords to new WCS ID in header */
2713 nkwdw = 0;
2714 for (i = 0; i < nkwd; i++) {
2715 if (dkwd[i]) {
2716 if (hgetr8 (*header, kwd[i], &tnum)) {
2717 nkwdw++;
2718 if (!strncmp (kwd[i], "PC0", 3)) {
2719 if (!strcmp (kwd[i], "PC001001"))
2720 strcpy (kwdc, "PC1_1");
2721 else if (!strcmp (kwd[i], "PC001002"))
2722 strcpy (kwdc, "PC1_2");
2723 else if (!strcmp (kwd[i], "PC002001"))
2724 strcpy (kwdc, "PC2_1");
2725 else
2726 strcpy (kwdc, "PC2_2");
2727 }
2728 else
2729 strcpy (kwdc, kwd[i]);
2730 strncat (kwdc, cwcs, 1);
2731 (void)hputr8 (*header, kwdc, tnum);
2732 }
2733 }
2734 else {
2735 if (hgets (*header, kwd[i], 80, tstr)) {
2736 nkwdw++;
2737 if (!strncmp (kwd[i], "RADECSYS", 8))
2738 strcpy (kwdc, "RADECSY");
2739 else
2740 strcpy (kwdc, kwd[i]);
2741 strncat (kwdc, cwcs, 1);
2742 hputs (*header, kwdc, tstr);
2743 }
2744 }
2745 }
2746
2747 /* Free keyword list array */
2748 for (ikwd = 0; ikwd < maxnkwd; ikwd++)
2749 free (kwd[ikwd]);
2750 free (kwd);
2751 kwd = NULL;
2752 return (nkwdw);
2753 }
2754
2755
2756 /* Oct 28 1994 new program
2757 * Dec 21 1994 Implement CD rotation matrix
2758 * Dec 22 1994 Allow RA and DEC to be either x,y or y,x
2759 *
2760 * Mar 6 1995 Add Digital Sky Survey plate fit
2761 * May 2 1995 Add prototype of PIX2WCST to WCSCOM
2762 * May 25 1995 Print leading zero for hours and degrees
2763 * Jun 21 1995 Add WCS2PIX to get pixels from WCS
2764 * Jun 21 1995 Read plate scale from FITS header for plate solution
2765 * Jul 6 1995 Pass WCS structure as argument; malloc it in WCSINIT
2766 * Jul 6 1995 Check string lengths in PIX2WCST
2767 * Aug 16 1995 Add galactic coordinate conversion to PIX2WCST
2768 * Aug 17 1995 Return 0 from iswcs if wcs structure is not yet set
2769 * Sep 8 1995 Do not include malloc.h if VMS
2770 * Sep 8 1995 Check for legal WCS before trying anything
2771 * Sep 8 1995 Do not try to set WCS if missing key keywords
2772 * Oct 18 1995 Add WCSCENT and WCSDIST to print center and size of image
2773 * Nov 6 1995 Include stdlib.h instead of malloc.h
2774 * Dec 6 1995 Fix format statement in PIX2WCST
2775 * Dec 19 1995 Change MALLOC to CALLOC to initialize array to zeroes
2776 * Dec 19 1995 Explicitly initialize rotation matrix and yinc
2777 * Dec 22 1995 If SECPIX is set, use approximate WCS
2778 * Dec 22 1995 Always print coordinate system
2779 *
2780 * Jan 12 1996 Use plane-tangent, not linear, projection if SECPIX is set
2781 * Jan 12 1996 Add WCSSET to set WCS without an image
2782 * Feb 15 1996 Replace all calls to HGETC with HGETS
2783 * Feb 20 1996 Add tab table output from PIX2WCST
2784 * Apr 2 1996 Convert all equinoxes to B1950 or J2000
2785 * Apr 26 1996 Get and use image epoch for accurate FK4/FK5 conversions
2786 * May 16 1996 Clean up internal documentation
2787 * May 17 1996 Return width in right ascension degrees, not sky degrees
2788 * May 24 1996 Remove extraneous print command from WCSSIZE
2789 * May 28 1996 Add NOWCS and WCSSHIFT subroutines
2790 * Jun 11 1996 Drop unused variables after running lint
2791 * Jun 12 1996 Set equinox as well as system in WCSSHIFT
2792 * Jun 14 1996 Make DSS keyword searches more robust
2793 * Jul 1 1996 Allow for SECPIX1 and SECPIX2 keywords
2794 * Jul 2 1996 Test for CTYPE1 instead of CRVAL1
2795 * Jul 5 1996 Declare all subroutines in wcs.h
2796 * Jul 19 1996 Add subroutine WCSFULL to return real image size
2797 * Aug 12 1996 Allow systemless coordinates which cannot be converted
2798 * Aug 15 1996 Allow LINEAR WCS to pass numbers through transparently
2799 * Aug 15 1996 Add WCSERR to print error message under calling program control
2800 * Aug 16 1996 Add latitude and longitude as image coordinate types
2801 * Aug 26 1996 Fix arguments to HLENGTH in WCSNINIT
2802 * Aug 28 1996 Explicitly set OFFSCL in WCS2PIX if coordinates outside image
2803 * Sep 3 1996 Return computed pixel values even if they are offscale
2804 * Sep 6 1996 Allow filename to be passed by WCSCOM
2805 * Oct 8 1996 Default to 2000 for EQUINOX and EPOCH and FK5 for RADECSYS
2806 * Oct 8 1996 If EPOCH is 0 and EQUINOX is not set, default to 1950 and FK4
2807 * Oct 15 1996 Add comparison when testing an assignment
2808 * Oct 16 1996 Allow PIXEL CTYPE which means WCS is same as image coordinates
2809 * Oct 21 1996 Add WCS_COMMAND environment variable
2810 * Oct 25 1996 Add image scale to WCSCENT
2811 * Oct 30 1996 Fix bugs in WCS2PIX
2812 * Oct 31 1996 Fix CD matrix rotation angle computation
2813 * Oct 31 1996 Use inline degree <-> radian conversion functions
2814 * Nov 1 1996 Add option to change number of decimal places in PIX2WCST
2815 * Nov 5 1996 Set wcs->crot to 1 if rotation matrix is used
2816 * Dec 2 1996 Add altitide/azimuth coordinates
2817 * Dec 13 1996 Fix search format setting from environment
2818 *
2819 * Jan 22 1997 Add ifdef for Eric Mandel (SAOtng)
2820 * Feb 5 1997 Add wcsout for Eric Mandel
2821 * Mar 20 1997 Drop unused variable STR in WCSCOM
2822 * May 21 1997 Do not make pixel coordinates mod 360 in PIX2WCST
2823 * May 22 1997 Add PIXEL prjcode = -1;
2824 * Jul 11 1997 Get center pixel x and y from header even if no WCS
2825 * Aug 7 1997 Add NOAO PIXSCALi keywords for default WCS
2826 * Oct 15 1997 Do not reset reference pixel in WCSSHIFT
2827 * Oct 20 1997 Set chip rotation
2828 * Oct 24 1997 Keep longitudes between 0 and 360, not -180 and +180
2829 * Nov 5 1997 Do no compute crot and srot; they are now computed in worldpos
2830 * Dec 16 1997 Set rotation and axis increments from CD matrix
2831 *
2832 * Jan 6 1998 Deal with J2000 and B1950 as EQUINOX values (from ST)
2833 * Jan 7 1998 Read INSTRUME and DETECTOR header keywords
2834 * Jan 7 1998 Fix tab-separated output
2835 * Jan 9 1998 Precess coordinates if either FITS projection or *DSS plate*
2836 * Jan 16 1998 Change PTYPE to not include initial hyphen
2837 * Jan 16 1998 Change WCSSET to WCSXINIT to avoid conflict with Calabretta
2838 * Jan 23 1998 Change PCODE to PRJCODE to avoid conflict with Calabretta
2839 * Jan 27 1998 Add LATPOLE and LONGPOLE for Calabretta projections
2840 * Feb 5 1998 Make cd and dc into vectors; use matinv() to invert cd
2841 * Feb 5 1998 In wcsoutinit(), check that corsys is a valid pointer
2842 * Feb 18 1998 Fix bugs for Calabretta projections
2843 * Feb 19 1998 Add wcs structure access subroutines from Eric Mandel
2844 * Feb 19 1998 Add wcsreset() to make sure derived values are reset
2845 * Feb 19 1998 Always set oldwcs to 1 if NCP projection
2846 * Feb 19 1998 Add subroutine to set oldwcs default
2847 * Feb 20 1998 Initialize projection types one at a time for SunOS C
2848 * Feb 23 1998 Add TNX projection from NOAO; treat it as TAN
2849 * Feb 23 1998 Compute size based on max and min coordinates, not sides
2850 * Feb 26 1998 Add code to set center pixel if part of detector array
2851 * Mar 6 1998 Write 8-character values to RADECSYS
2852 * Mar 9 1998 Add naxis to WCS structure
2853 * Mar 16 1998 Use separate subroutine for IRAF TNX projection
2854 * Mar 20 1998 Set PC matrix if more than two axes and it's not in header
2855 * Mar 20 1998 Reset lin flag in WCSRESET if CDELTn
2856 * Mar 20 1998 Set CD matrix with CDELTs and CROTA in wcsinit and wcsreset
2857 * Mar 20 1998 Allow initialization of rotation angle alone
2858 * Mar 23 1998 Use dsspos() and dsspix() instead of platepos() and platepix()
2859 * Mar 24 1998 Set up PLT/PLATE plate polynomial fit using platepos() and platepix()
2860 * Mar 25 1998 Read plate fit coefficients from header in getwcscoeffs()
2861 * Mar 27 1998 Check for FITS WCS before DSS WCS
2862 * Mar 27 1998 Compute scale from edges if xinc and yinc not set in wcscent()
2863 * Apr 6 1998 Change plate coefficient keywords from PLTij to COi_j
2864 * Apr 6 1998 Read all coefficients in line instead of with subroutine
2865 * Apr 7 1998 Change amd_i_coeff to i_coeff
2866 * Apr 8 1998 Add wcseqset to change equinox after wcs has been set
2867 * Apr 10 1998 Use separate counters for x and y polynomial coefficients
2868 * Apr 13 1998 Use CD/CDELT+CDROTA if oldwcs is set
2869 * Apr 14 1998 Use codes instead of strings for various coordinate systems
2870 * Apr 14 1998 Separate input coordinate conversion from output conversion
2871 * Apr 14 1998 Use wcscon() for most coordinate conversion
2872 * Apr 17 1998 Always compute cdelt[]
2873 * Apr 17 1998 Deal with reversed axis more correctly
2874 * Apr 17 1998 Compute rotation angle and approximate CDELTn for polynomial
2875 * Apr 23 1998 Deprecate xref/yref in favor of crval[]
2876 * Apr 23 1998 Deprecate xrefpix/yrefpix in favor of crpix[]
2877 * Apr 23 1998 Add LINEAR to coordinate system types
2878 * Apr 23 1998 Always use AIPS subroutines for LINEAR or PIXEL
2879 * Apr 24 1998 Format linear coordinates better
2880 * Apr 28 1998 Change coordinate system flags to WCS_*
2881 * Apr 28 1998 Change projection flags to WCS_*
2882 * Apr 28 1998 Add subroutine wcsc2pix for coordinates to pixels with system
2883 * Apr 28 1998 Add setlinmode() to set output string mode for LINEAR coordinates
2884 * Apr 30 1998 Fix bug by setting degree flag for lat and long in wcsinit()
2885 * Apr 30 1998 Allow leading "-"s in projecting in wcsxinit()
2886 * May 1 1998 Assign new output coordinate system only if legitimate system
2887 * May 1 1998 Do not allow oldwcs=1 unless allowed projection
2888 * May 4 1998 Fix bug in units reading for LINEAR coordinates
2889 * May 6 1998 Initialize to no CD matrix
2890 * May 6 1998 Use TAN instead of TNX if oldwcs
2891 * May 12 1998 Set 3rd and 4th coordinates in wcspos()
2892 * May 12 1998 Return *xpos and *ypos = 0 in pix2wcs() if offscale
2893 * May 12 1998 Declare undeclared external subroutines after lint
2894 * May 13 1998 Add equinox conversion to specified output equinox
2895 * May 13 1998 Set output or input system to image with null argument
2896 * May 15 1998 Return reference pixel, cdelts, and rotation for DSS
2897 * May 20 1998 Fix bad bug so setlinmode() is no-op if wcs not set
2898 * May 20 1998 Fix bug so getwcsout() returns null pointer if wcs not set
2899 * May 27 1998 Change WCS_LPR back to WCS_LIN; allow CAR in oldwcs
2900 * May 28 1998 Go back to old WCSFULL, computing height and width from center
2901 * May 29 1998 Add wcskinit() to initialize WCS from arguments
2902 * May 29 1998 Add wcstype() to set projection from arguments
2903 * May 29 1998 Add wcscdset(), and wcsdeltset() to set scale from arguments
2904 * Jun 1 1998 In wcstype(), reconstruct ctype for WCS structure
2905 * Jun 11 1998 Split off header-dependent subroutines to wcsinit.c
2906 * Jun 18 1998 Add wcspcset() for PC matrix initialization
2907 * Jun 24 1998 Add string lengths to ra2str(), dec2str, and deg2str() calls
2908 * Jun 25 1998 Use AIPS software for CAR projection
2909 * Jun 25 1998 Add wcsndec to set number of decimal places in output string
2910 * Jul 6 1998 Add wcszin() and wcszout() to use third dimension of images
2911 * Jul 7 1998 Change setlinmode() to setwcslin(); setdegout() to setwcsdeg()
2912 * Jul 10 1998 Initialize matrices correctly for naxis > 2 in wcs<>set()
2913 * Jul 16 1998 Initialize coordinates to be returned in wcspos()
2914 * Jul 17 1998 Link lin structure arrays to wcs structure arrays
2915 * Jul 20 1998 In wcscdset() compute sign of xinc and yinc from CD1_1, CD 2_2
2916 * Jul 20 1998 In wcscdset() compute sign of rotation based on CD1_1, CD 1_2
2917 * Jul 22 1998 Add wcslibrot() to compute lin() rotation matrix
2918 * Jul 30 1998 Set wcs->naxes and lin.naxis in wcsxinit() and wcskinit()
2919 * Aug 5 1998 Use old WCS subroutines to deal with COE projection (for ESO)
2920 * Aug 14 1998 Add option to print image coordinates with wcscom()
2921 * Aug 14 1998 Add multiple command options to wcscom*()
2922 * Aug 31 1998 Declare undeclared arguments to wcspcset()
2923 * Sep 3 1998 Set CD rotation and cdelts from sky axis position angles
2924 * Sep 16 1998 Add option to use North Polar Angle instead of Latitude
2925 * Sep 29 1998 Initialize additional WCS commands from the environment
2926 * Oct 14 1998 Fix bug in wcssize() which didn't divide dra by cos(dec)
2927 * Nov 12 1998 Fix sign of CROTA when either axis is reflected
2928 * Dec 2 1998 Fix non-arcsecond scale factors in wcscent()
2929 * Dec 2 1998 Add PLANET coordinate system to pix2wcst()
2930
2931 * Jan 20 1999 Free lin.imgpix and lin.piximg in wcsfree()
2932 * Feb 22 1999 Fix bug setting latitude reference value of latbase != 0
2933 * Feb 22 1999 Fix bug so that quad cube faces are 0-5, not 1-6
2934 * Mar 16 1999 Always initialize all 4 imgcrds and pixcrds in wcspix()
2935 * Mar 16 1999 Always return (0,0) from wcs2pix if offscale
2936 * Apr 7 1999 Add code to put file name in error messages
2937 * Apr 7 1999 Document utility subroutines at end of file
2938 * May 6 1999 Fix bug printing height of LINEAR image
2939 * Jun 16 1999 Add wcsrange() to return image RA and Dec limits
2940 * Jul 8 1999 Always use FK5 and FK4 instead of J2000 and B1950 in RADECSYS
2941 * Aug 16 1999 Print dd:mm:ss dd:mm:ss if not J2000 or B1950 output
2942 * Aug 20 1999 Add WCS string argument to wcscom(); don't compute it there
2943 * Aug 20 1999 Change F3 WCS command default from Tycho to ACT
2944 * Oct 15 1999 Free wcs using wcsfree()
2945 * Oct 21 1999 Drop declarations of unused functions after lint
2946 * Oct 25 1999 Drop out of wcsfree() if wcs is null pointer
2947 * Nov 17 1999 Fix bug which caused software to miss NCP projection
2948 *
2949 * Jan 24 2000 Default to AIPS for NCP, CAR, and COE proj.; if -z use WCSLIB
2950 * Feb 24 2000 If coorsys is null in wcsc2pix, wcs->radecin is assumed
2951 * May 10 2000 In wcstype(), default to WCS_LIN, not error (after Bill Joye)
2952 * Jun 22 2000 In wcsrotset(), leave rotation angle alone in 1-d image
2953 * Jul 3 2000 Initialize wcscrd[] to zero in wcspix()
2954 *
2955 * Feb 20 2001 Add recursion to wcs2pix() and pix2wcs() for dependent WCS's
2956 * Mar 20 2001 Add braces to avoid ambiguity in if/else groupings
2957 * Mar 22 2001 Free WCS structure in wcsfree even if it is not filled
2958 * Sep 12 2001 Fix bug which omitted tab in pix2wcst() galactic coord output
2959 *
2960 * Mar 7 2002 Fix bug which gave wrong pa's and rotation for reflected RA
2961 * (but correct WCS conversions!)
2962 * Mar 28 2002 Add SZP projection
2963 * Apr 3 2002 Synchronize projection types with other subroutines
2964 * Apr 3 2002 Drop special cases of projections
2965 * Apr 9 2002 Implement inversion of multiple WCSs in wcsc2pix()
2966 * Apr 25 2002 Use Tycho-2 catalog instead of ACT in setwcscom()
2967 * May 13 2002 Free WCSNAME in wcsfree()
2968 *
2969 * Mar 31 2003 Add distcode to wcstype()
2970 * Apr 1 2003 Add calls to foc2pix() in wcs2pix() and pix2foc() in pix2wcs()
2971 * May 20 2003 Declare argument i in savewcscom()
2972 * Sep 29 2003 Fix bug to compute width and height correctly in wcsfull()
2973 * Sep 29 2003 Fix bug to deal with all-sky images orrectly in wcsfull()
2974 * Oct 1 2003 Rename wcs->naxes to wcs->naxis to match WCSLIB 3.2
2975 * Nov 3 2003 Set distortion code by calling setdistcode() in wcstype()
2976 * Dec 3 2003 Add back wcs->naxes for compatibility
2977 * Dec 3 2003 Add braces in if...else in pix2wcst()
2978 *
2979 * Sep 17 2004 If spherical coordinate output, keep 0 < long/RA < 360
2980 * Sep 17 2004 Fix bug in wcsfull() when wrapping around RA=0:00
2981 * Nov 1 2004 Keep wcs->rot between 0 and 360
2982 *
2983 * Mar 9 2005 Fix bug in wcsrotset() which set rot > 360 to 360
2984 * Jun 27 2005 Fix ctype in calls to wcs subroutines
2985 * Jul 21 2005 Fix bug in wcsrange() at RA ~= 0.0
2986 *
2987 * Apr 24 2006 Always set inverse CD matrix to 2 dimensions in wcspcset()
2988 * May 3 2006 (void *) pointers so arguments match type, from Robert Lupton
2989 * Jun 30 2006 Set only 2-dimensional PC matrix; that is all lin* can deal with
2990 * Oct 30 2006 In pix2wcs(), do not limit x to between 0 and 360 if XY or LINEAR
2991 * Oct 30 2006 In wcsc2pix(), do not precess LINEAR or XY coordinates
2992 * Dec 21 2006 Add cpwcs() to copy WCS keywords to new suffix
2993 *
2994 * Jan 4 2007 Fix pointer to header in cpwcs()
2995 * Jan 5 2007 Drop declarations of wcscon(); it is in wcs.h
2996 * Jan 9 2006 Drop declarations of fk425e() and fk524e(); moved to wcs.h
2997 * Jan 9 2006 Drop *pix() and *pos() external declarations; moved to wcs.h
2998 * Jan 9 2006 Drop matinv() external declaration; it is already in wcslib.h
2999 * Feb 15 2007 If CTYPEi contains DET, set to WCS_PIX projection
3000 * Feb 23 2007 Fix bug when checking for "DET" in CTYPEi
3001 * Apr 2 2007 Fix PC to CD matrix conversion
3002 * Jul 25 2007 Compute distance between two coordinates using d2v3()
3003 *
3004 * Apr 7 2010 In wcstype() set number of WCS projections from NWCSTYPE
3005 *
3006 * Mar 11 2011 Add NOAO ZPX projection (Frank Valdes)
3007 * Mar 14 2011 Delete j<=MAXPV PVi_j parameters (for SCAMP polynomials via Ed Los)
3008 * Mar 17 2011 Fix WCSDEP bug found by Ed Los
3009 * May 9 2011 Free WCS structure recursively if WCSDEP is used
3010 * Sep 1 2011 Add TPV projection type for SCAMP TAN with PVs
3011 *
3012 * Oct 19 2012 Drop d1 and d2 from wcsdist(); diffi from wcsdist1()
3013 * Oct 19 2012 Drop depwcs; it's in main wcs structure
3014 *
3015 * Jun 8 2016 Increase ctype, ctype1, and ctype2 to 16 characters for distortion
3016 * Jun 23 2016 Set initial allocation of keyword arrays to MAXNKWD instead of 100 in cpwcs()
3017 * Jun 24 2016 wcs->ptype contains only 3-letter projection code
3018 */
3019