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