/*** File libwcs/wcs.c
*** October 19, 2012
*** By Jessica Mink, jmink@cfa.harvard.edu
*** Harvard-Smithsonian Center for Astrophysics
*** Copyright (C) 1994-2012
*** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Correspondence concerning WCSTools should be addressed as follows:
Internet email: jmink@cfa.harvard.edu
Postal address: Jessica Mink
Smithsonian Astrophysical Observatory
60 Garden St.
Cambridge, MA 02138 USA
* Module: wcs.c (World Coordinate Systems)
* Purpose: Convert FITS WCS to pixels and vice versa:
* Subroutine: wcsxinit (cra,cdec,secpix,xrpix,yrpix,nxpix,nypix,rotate,equinox,epoch,proj)
* sets a WCS structure from arguments
* Subroutine: wcskinit (nxpix,nypix,ctype1,ctype2,crpix1,crpix2,crval1,crval2,
cd,cdelt1,cdelt2,crota,equinox,epoch)
* sets a WCS structure from keyword-based arguments
* Subroutine: wcsreset (wcs,crpix1,crpix2,crval1,crval2,cdelt1,cdelt2,crota,cd, equinox)
* resets an existing WCS structure from arguments
* Subroutine: wcsdeltset (wcs,cdelt1,cdelt2,crota) sets rotation and scaling
* Subroutine: wcscdset (wcs, cd) sets rotation and scaling from a CD matrix
* Subroutine: wcspcset (wcs,cdelt1,cdelt2,pc) sets rotation and scaling
* Subroutine: wcseqset (wcs, equinox) resets an existing WCS structure to new equinox
* Subroutine: iswcs(wcs) returns 1 if WCS structure is filled, else 0
* Subroutine: nowcs(wcs) returns 0 if WCS structure is filled, else 1
* Subroutine: wcscent (wcs) prints the image center and size in WCS units
* Subroutine: wcssize (wcs, cra, cdec, dra, ddec) returns image center and size
* Subroutine: wcsfull (wcs, cra, cdec, width, height) returns image center and size
* Subroutine: wcsrange (wcs, ra1, ra2, dec1, dec2) returns image coordinate limits
* Subroutine: wcsshift (wcs,cra,cdec) resets the center of a WCS structure
* Subroutine: wcsdist (x1,y1,x2,y2) compute angular distance between ra/dec or lat/long
* Subroutine: wcsdiff (x1,y1,x2,y2) compute angular distance between ra/dec or lat/long
* Subroutine: wcscominit (wcs,command) sets up a command format for execution by wcscom
* Subroutine: wcsoutinit (wcs,coor) sets up the coordinate system used by pix2wcs
* Subroutine: getwcsout (wcs) returns current output coordinate system used by pix2wcs
* Subroutine: wcsininit (wcs,coor) sets up the coordinate system used by wcs2pix
* Subroutine: getwcsin (wcs) returns current input coordinate system used by wcs2pix
* Subroutine: setwcsdeg(wcs, new) sets WCS output in degrees or hh:mm:ss
* Subroutine: getradecsys(wcs) returns current coordinate system type
* Subroutine: wcscom (wcs,file,x,y,wcstr) executes a command using the current world coordinates
* Subroutine: setwcslin (wcs, mode) sets output string mode for LINEAR
* Subroutine: pix2wcst (wcs,xpix,ypix,wcstring,lstr) pixels -> sky coordinate string
* Subroutine: pix2wcs (wcs,xpix,ypix,xpos,ypos) pixel coordinates -> sky coordinates
* Subroutine: wcsc2pix (wcs,xpos,ypos,coorsys,xpix,ypix,offscl) sky coordinates -> pixel coordinates
* Subroutine: wcs2pix (wcs,xpos,ypos,xpix,ypix,offscl) sky coordinates -> pixel coordinates
* Subroutine: wcszin (izpix) sets third dimension for pix2wcs() and pix2wcst()
* Subroutine: wcszout (wcs) returns third dimension from wcs2pix()
* Subroutine: setwcsfile (filename) Set file name for error messages
* Subroutine: setwcserr (errmsg) Set error message
* Subroutine: wcserr() Print error message
* Subroutine: setdefwcs (wcsproj) Set flag to choose AIPS or WCSLIB WCS subroutines
* Subroutine: getdefwcs() Get flag to switch between AIPS and WCSLIB subroutines
* Subroutine: savewcscoor (wcscoor)
* Subroutine: getwcscoor() Return preset output default coordinate system
* Subroutine: savewcscom (i, wcscom) Save specified WCS command
* Subroutine: setwcscom (wcs) Initialize WCS commands
* Subroutine: getwcscom (i) Return specified WCS command
* Subroutine: wcsfree (wcs) Free storage used by WCS structure
* Subroutine: freewcscom (wcs) Free storage used by WCS commands
* Subroutine: cpwcs (&header, cwcs)
*/
#include <string.h> /* strstr, NULL */
#include <stdio.h> /* stderr */
#include <math.h>
#include "wcs.h"
#ifndef VMS
#include <stdlib.h>
#endif
static char wcserrmsg[80];
static char wcsfile[256]={""};
static void wcslibrot();
void wcsrotset();
static int wcsproj0 = 0;
static int izpix = 0;
static double zpix = 0.0;
void
wcsfree (wcs)
struct WorldCoor *wcs; /* WCS structure */
{
if (nowcs (wcs)) {
/* Free WCS structure if allocated but not filled */
if (wcs)
free (wcs);
return;
}
/* Free WCS on which this WCS depends */
if (wcs->wcs) {
wcsfree (wcs->wcs);
wcs->wcs = NULL;
}
freewcscom (wcs);
if (wcs->wcsname != NULL)
free (wcs->wcsname);
if (wcs->lin.imgpix != NULL)
free (wcs->lin.imgpix);
if (wcs->lin.piximg != NULL)
free (wcs->lin.piximg);
if (wcs->inv_x != NULL)
poly_end (wcs->inv_x);
if (wcs->inv_y != NULL)
poly_end (wcs->inv_y);
free (wcs);
return;
}
/* Set up a WCS structure from subroutine arguments */
struct WorldCoor *
wcsxinit (cra,cdec,secpix,xrpix,yrpix,nxpix,nypix,rotate,equinox,epoch,proj)
double cra; /* Center right ascension in degrees */
double cdec; /* Center declination in degrees */
double secpix; /* Number of arcseconds per pixel */
double xrpix; /* Reference pixel X coordinate */
double yrpix; /* Reference pixel X coordinate */
int nxpix; /* Number of pixels along x-axis */
int nypix; /* Number of pixels along y-axis */
double rotate; /* Rotation angle (clockwise positive) in degrees */
int equinox; /* Equinox of coordinates, 1950 and 2000 supported */
double epoch; /* Epoch of coordinates, used for FK4/FK5 conversion
* no effect if 0 */
char *proj; /* Projection */
{
struct WorldCoor *wcs;
double cdelt1, cdelt2;
wcs = (struct WorldCoor *) calloc (1, sizeof(struct WorldCoor));
/* Set WCSLIB flags so that structures will be reinitialized */
wcs->cel.flag = 0;
wcs->lin.flag = 0;
wcs->wcsl.flag = 0;
/* Image dimensions */
wcs->naxis = 2;
wcs->naxes = 2;
wcs->lin.naxis = 2;
wcs->nxpix = nxpix;
wcs->nypix = nypix;
wcs->wcsproj = wcsproj0;
wcs->crpix[0] = xrpix;
wcs->crpix[1] = yrpix;
wcs->xrefpix = wcs->crpix[0];
wcs->yrefpix = wcs->crpix[1];
wcs->lin.crpix = wcs->crpix;
wcs->crval[0] = cra;
wcs->crval[1] = cdec;
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
wcs->cel.ref[0] = wcs->crval[0];
wcs->cel.ref[1] = wcs->crval[1];
wcs->cel.ref[2] = 999.0;
strcpy (wcs->c1type,"RA");
strcpy (wcs->c2type,"DEC");
/* Allan Brighton: 28.4.98: for backward compat., remove leading "--" */
while (proj && *proj == '-')
proj++;
strcpy (wcs->ptype,proj);
strcpy (wcs->ctype[0],"RA---");
strcpy (wcs->ctype[1],"DEC--");
strcat (wcs->ctype[0],proj);
strcat (wcs->ctype[1],proj);
if (wcstype (wcs, wcs->ctype[0], wcs->ctype[1])) {
wcsfree (wcs);
return (NULL);
}
/* Approximate world coordinate system from a known plate scale */
cdelt1 = -secpix / 3600.0;
cdelt2 = secpix / 3600.0;
wcsdeltset (wcs, cdelt1, cdelt2, rotate);
wcs->lin.cdelt = wcs->cdelt;
wcs->lin.pc = wcs->pc;
/* Coordinate reference frame and equinox */
wcs->equinox = (double) equinox;
if (equinox > 1980)
strcpy (wcs->radecsys,"FK5");
else
strcpy (wcs->radecsys,"FK4");
if (epoch > 0)
wcs->epoch = epoch;
else
wcs->epoch = 0.0;
wcs->wcson = 1;
wcs->syswcs = wcscsys (wcs->radecsys);
wcsoutinit (wcs, wcs->radecsys);
wcsininit (wcs, wcs->radecsys);
wcs->eqout = 0.0;
wcs->printsys = 1;
wcs->tabsys = 0;
/* Initialize special WCS commands */
setwcscom (wcs);
return (wcs);
}
/* Set up a WCS structure from subroutine arguments based on FITS keywords */
struct WorldCoor *
wcskinit (naxis1, naxis2, ctype1, ctype2, crpix1, crpix2, crval1, crval2,
cd, cdelt1, cdelt2, crota, equinox, epoch)
int naxis1; /* Number of pixels along x-axis */
int naxis2; /* Number of pixels along y-axis */
char *ctype1; /* FITS WCS projection for axis 1 */
char *ctype2; /* FITS WCS projection for axis 2 */
double crpix1, crpix2; /* Reference pixel coordinates */
double crval1, crval2; /* Coordinates at reference pixel in degrees */
double *cd; /* Rotation matrix, used if not NULL */
double cdelt1, cdelt2; /* scale in degrees/pixel, ignored if cd is not NULL */
double crota; /* Rotation angle in degrees, ignored if cd is not NULL */
int equinox; /* Equinox of coordinates, 1950 and 2000 supported */
double epoch; /* Epoch of coordinates, used for FK4/FK5 conversion
* no effect if 0 */
{
struct WorldCoor *wcs;
wcs = (struct WorldCoor *) calloc (1, sizeof(struct WorldCoor));
/* Set WCSLIB flags so that structures will be reinitialized */
wcs->cel.flag = 0;
wcs->lin.flag = 0;
wcs->wcsl.flag = 0;
/* Image dimensions */
wcs->naxis = 2;
wcs->naxes = 2;
wcs->lin.naxis = 2;
wcs->nxpix = naxis1;
wcs->nypix = naxis2;
wcs->wcsproj = wcsproj0;
wcs->crpix[0] = crpix1;
wcs->crpix[1] = crpix2;
wcs->xrefpix = wcs->crpix[0];
wcs->yrefpix = wcs->crpix[1];
wcs->lin.crpix = wcs->crpix;
if (wcstype (wcs, ctype1, ctype2)) {
wcsfree (wcs);
return (NULL);
}
if (wcs->latbase == 90)
crval2 = 90.0 - crval2;
else if (wcs->latbase == -90)
crval2 = crval2 - 90.0;
wcs->crval[0] = crval1;
wcs->crval[1] = crval2;
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
wcs->cel.ref[0] = wcs->crval[0];
wcs->cel.ref[1] = wcs->crval[1];
wcs->cel.ref[2] = 999.0;
if (cd != NULL)
wcscdset (wcs, cd);
else if (cdelt1 != 0.0)
wcsdeltset (wcs, cdelt1, cdelt2, crota);
else {
wcsdeltset (wcs, 1.0, 1.0, crota);
setwcserr ("WCSRESET: setting CDELT to 1");
}
wcs->lin.cdelt = wcs->cdelt;
wcs->lin.pc = wcs->pc;
/* Coordinate reference frame and equinox */
wcs->equinox = (double) equinox;
if (equinox > 1980)
strcpy (wcs->radecsys,"FK5");
else
strcpy (wcs->radecsys,"FK4");
if (epoch > 0)
wcs->epoch = epoch;
else
wcs->epoch = 0.0;
wcs->wcson = 1;
strcpy (wcs->radecout, wcs->radecsys);
wcs->syswcs = wcscsys (wcs->radecsys);
wcsoutinit (wcs, wcs->radecsys);
wcsininit (wcs, wcs->radecsys);
wcs->eqout = 0.0;
wcs->printsys = 1;
wcs->tabsys = 0;
/* Initialize special WCS commands */
setwcscom (wcs);
return (wcs);
}
/* Set projection in WCS structure from FITS keyword values */
int
wcstype (wcs, ctype1, ctype2)
struct WorldCoor *wcs; /* World coordinate system structure */
char *ctype1; /* FITS WCS projection for axis 1 */
char *ctype2; /* FITS WCS projection for axis 2 */
{
int i, iproj;
int nctype = NWCSTYPE;
char ctypes[NWCSTYPE][4];
char dtypes[10][4];
/* Initialize projection types */
strcpy (ctypes[0], "LIN");
strcpy (ctypes[1], "AZP");
strcpy (ctypes[2], "SZP");
strcpy (ctypes[3], "TAN");
strcpy (ctypes[4], "SIN");
strcpy (ctypes[5], "STG");
strcpy (ctypes[6], "ARC");
strcpy (ctypes[7], "ZPN");
strcpy (ctypes[8], "ZEA");
strcpy (ctypes[9], "AIR");
strcpy (ctypes[10], "CYP");
strcpy (ctypes[11], "CAR");
strcpy (ctypes[12], "MER");
strcpy (ctypes[13], "CEA");
strcpy (ctypes[14], "COP");
strcpy (ctypes[15], "COD");
strcpy (ctypes[16], "COE");
strcpy (ctypes[17], "COO");
strcpy (ctypes[18], "BON");
strcpy (ctypes[19], "PCO");
strcpy (ctypes[20], "SFL");
strcpy (ctypes[21], "PAR");
strcpy (ctypes[22], "AIT");
strcpy (ctypes[23], "MOL");
strcpy (ctypes[24], "CSC");
strcpy (ctypes[25], "QSC");
strcpy (ctypes[26], "TSC");
strcpy (ctypes[27], "NCP");
strcpy (ctypes[28], "GLS");
strcpy (ctypes[29], "DSS");
strcpy (ctypes[30], "PLT");
strcpy (ctypes[31], "TNX");
strcpy (ctypes[32], "ZPX");
strcpy (ctypes[33], "TPV");
/* Initialize distortion types */
strcpy (dtypes[1], "SIP");
if (!strncmp (ctype1, "LONG",4))
strncpy (ctype1, "XLON",4);
strcpy (wcs->ctype[0], ctype1);
strcpy (wcs->c1type, ctype1);
strcpy (wcs->ptype, ctype1);
/* Linear coordinates */
if (!strncmp (ctype1,"LINEAR",6))
wcs->prjcode = WCS_LIN;
/* Pixel coordinates */
else if (!strncmp (ctype1,"PIXEL",6))
wcs->prjcode = WCS_PIX;
/*Detector pixel coordinates */
else if (strsrch (ctype1,"DET"))
wcs->prjcode = WCS_PIX;
/* Set up right ascension, declination, latitude, or longitude */
else if (ctype1[0] == 'R' || ctype1[0] == 'D' ||
ctype1[0] == 'A' || ctype1[1] == 'L') {
wcs->c1type[0] = ctype1[0];
wcs->c1type[1] = ctype1[1];
if (ctype1[2] == '-') {
wcs->c1type[2] = 0;
iproj = 3;
}
else {
wcs->c1type[2] = ctype1[2];
iproj = 4;
if (ctype1[3] == '-') {
wcs->c1type[3] = 0;
}
else {
wcs->c1type[3] = ctype1[3];
wcs->c1type[4] = 0;
}
}
if (ctype1[iproj] == '-') iproj = iproj + 1;
if (ctype1[iproj] == '-') iproj = iproj + 1;
if (ctype1[iproj] == '-') iproj = iproj + 1;
if (ctype1[iproj] == '-') iproj = iproj + 1;
wcs->ptype[0] = ctype1[iproj];
wcs->ptype[1] = ctype1[iproj+1];
wcs->ptype[2] = ctype1[iproj+2];
wcs->ptype[3] = 0;
sprintf (wcs->ctype[0],"%-4s%4s",wcs->c1type,wcs->ptype);
for (i = 0; i < 8; i++)
if (wcs->ctype[0][i] == ' ') wcs->ctype[0][i] = '-';
/* Find projection type */
wcs->prjcode = 0; /* default type is linear */
for (i = 1; i < nctype; i++) {
if (!strncmp(wcs->ptype, ctypes[i], 3))
wcs->prjcode = i;
}
/* Handle "obsolete" NCP projection (now WCSLIB should be OK)
if (wcs->prjcode == WCS_NCP) {
if (wcs->wcsproj == WCS_BEST)
wcs->wcsproj = WCS_OLD;
else if (wcs->wcsproj == WCS_ALT)
wcs->wcsproj = WCS_NEW;
} */
/* Work around bug in WCSLIB handling of CAR projection
else if (wcs->prjcode == WCS_CAR) {
if (wcs->wcsproj == WCS_BEST)
wcs->wcsproj = WCS_OLD;
else if (wcs->wcsproj == WCS_ALT)
wcs->wcsproj = WCS_NEW;
} */
/* Work around bug in WCSLIB handling of COE projection
else if (wcs->prjcode == WCS_COE) {
if (wcs->wcsproj == WCS_BEST)
wcs->wcsproj = WCS_OLD;
else if (wcs->wcsproj == WCS_ALT)
wcs->wcsproj = WCS_NEW;
}
else if (wcs->wcsproj == WCS_BEST) */
if (wcs->wcsproj == WCS_BEST)
wcs->wcsproj = WCS_NEW;
else if (wcs->wcsproj == WCS_ALT)
wcs->wcsproj = WCS_OLD;
/* if (wcs->wcsproj == WCS_OLD && (
wcs->prjcode != WCS_STG && wcs->prjcode != WCS_AIT &&
wcs->prjcode != WCS_MER && wcs->prjcode != WCS_GLS &&
wcs->prjcode != WCS_ARC && wcs->prjcode != WCS_TAN &&
wcs->prjcode != WCS_TNX && wcs->prjcode != WCS_SIN &&
wcs->prjcode != WCS_PIX && wcs->prjcode != WCS_LIN &&
wcs->prjcode != WCS_CAR && wcs->prjcode != WCS_COE &&
wcs->prjcode != WCS_NCP && wcs->prjcode != WCS_ZPX))
wcs->wcsproj = WCS_NEW; */
/* Handle NOAO corrected TNX as uncorrected TAN if oldwcs is set */
if (wcs->wcsproj == WCS_OLD && wcs->prjcode == WCS_TNX) {
wcs->ctype[0][6] = 'A';
wcs->ctype[0][7] = 'N';
wcs->prjcode = WCS_TAN;
}
/* Handle NOAO corrected ZPX as uncorrected ZPN if oldwcs is set */
if (wcs->wcsproj == WCS_OLD && wcs->prjcode == WCS_ZPX) {
wcs->ctype[0][6] = 'P';
wcs->ctype[0][7] = 'N';
wcs->prjcode = WCS_ZPN;
}
}
/* If not sky coordinates, assume linear */
else {
wcs->prjcode = WCS_LIN;
return (0);
}
/* Second coordinate type */
if (!strncmp (ctype2, "NPOL",4)) {
ctype2[0] = ctype1[0];
strncpy (ctype2+1, "LAT",3);
wcs->latbase = 90;
strcpy (wcs->radecsys,"NPOLE");
wcs->syswcs = WCS_NPOLE;
}
else if (!strncmp (ctype2, "SPA-",4)) {
ctype2[0] = ctype1[0];
strncpy (ctype2+1, "LAT",3);
wcs->latbase = -90;
strcpy (wcs->radecsys,"SPA");
wcs->syswcs = WCS_SPA;
}
else
wcs->latbase = 0;
strcpy (wcs->ctype[1], ctype2);
strcpy (wcs->c2type, ctype2);
/* Linear coordinates */
if (!strncmp (ctype2,"LINEAR",6))
wcs->prjcode = WCS_LIN;
/* Pixel coordinates */
else if (!strncmp (ctype2,"PIXEL",6))
wcs->prjcode = WCS_PIX;
/* Set up right ascension, declination, latitude, or longitude */
else if (ctype2[0] == 'R' || ctype2[0] == 'D' ||
ctype2[0] == 'A' || ctype2[1] == 'L') {
wcs->c2type[0] = ctype2[0];
wcs->c2type[1] = ctype2[1];
if (ctype2[2] == '-') {
wcs->c2type[2] = 0;
iproj = 3;
}
else {
wcs->c2type[2] = ctype2[2];
iproj = 4;
if (ctype2[3] == '-') {
wcs->c2type[3] = 0;
}
else {
wcs->c2type[3] = ctype2[3];
wcs->c2type[4] = 0;
}
}
if (ctype2[iproj] == '-') iproj = iproj + 1;
if (ctype2[iproj] == '-') iproj = iproj + 1;
if (ctype2[iproj] == '-') iproj = iproj + 1;
if (ctype2[iproj] == '-') iproj = iproj + 1;
wcs->ptype[0] = ctype2[iproj];
wcs->ptype[1] = ctype2[iproj+1];
wcs->ptype[2] = ctype2[iproj+2];
wcs->ptype[3] = 0;
if (!strncmp (ctype1, "DEC", 3) ||
!strncmp (ctype1+1, "LAT", 3))
wcs->coorflip = 1;
else
wcs->coorflip = 0;
if (ctype2[1] == 'L' || ctype2[0] == 'A') {
wcs->degout = 1;
wcs->ndec = 5;
}
else {
wcs->degout = 0;
wcs->ndec = 3;
}
sprintf (wcs->ctype[1],"%-4s%4s",wcs->c2type,wcs->ptype);
for (i = 0; i < 8; i++)
if (wcs->ctype[1][i] == ' ') wcs->ctype[1][i] = '-';
}
/* If not sky coordinates, assume linear */
else {
wcs->prjcode = WCS_LIN;
}
/* Set distortion code from CTYPE1 extension */
setdistcode (wcs, ctype1);
return (0);
}
int
wcsreset (wcs, crpix1, crpix2, crval1, crval2, cdelt1, cdelt2, crota, cd)
struct WorldCoor *wcs; /* World coordinate system data structure */
double crpix1, crpix2; /* Reference pixel coordinates */
double crval1, crval2; /* Coordinates at reference pixel in degrees */
double cdelt1, cdelt2; /* scale in degrees/pixel, ignored if cd is not NULL */
double crota; /* Rotation angle in degrees, ignored if cd is not NULL */
double *cd; /* Rotation matrix, used if not NULL */
{
if (nowcs (wcs))
return (-1);
/* Set WCSLIB flags so that structures will be reinitialized */
wcs->cel.flag = 0;
wcs->lin.flag = 0;
wcs->wcsl.flag = 0;
/* Reference pixel coordinates and WCS value */
wcs->crpix[0] = crpix1;
wcs->crpix[1] = crpix2;
wcs->xrefpix = wcs->crpix[0];
wcs->yrefpix = wcs->crpix[1];
wcs->lin.crpix = wcs->crpix;
wcs->crval[0] = crval1;
wcs->crval[1] = crval2;
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
if (wcs->coorflip) {
wcs->cel.ref[1] = wcs->crval[0];
wcs->cel.ref[0] = wcs->crval[1];
}
else {
wcs->cel.ref[0] = wcs->crval[0];
wcs->cel.ref[1] = wcs->crval[1];
}
/* Keep ref[2] and ref[3] from input */
/* Initialize to no plate fit */
wcs->ncoeff1 = 0;
wcs->ncoeff2 = 0;
if (cd != NULL)
wcscdset (wcs, cd);
else if (cdelt1 != 0.0)
wcsdeltset (wcs, cdelt1, cdelt2, crota);
else {
wcs->xinc = 1.0;
wcs->yinc = 1.0;
setwcserr ("WCSRESET: setting CDELT to 1");
}
/* Coordinate reference frame, equinox, and epoch */
if (!strncmp (wcs->ptype,"LINEAR",6) ||
!strncmp (wcs->ptype,"PIXEL",5))
wcs->degout = -1;
wcs->wcson = 1;
return (0);
}
void
wcseqset (wcs, equinox)
struct WorldCoor *wcs; /* World coordinate system data structure */
double equinox; /* Desired equinox as fractional year */
{
if (nowcs (wcs))
return;
/* Leave WCS alone if already at desired equinox */
if (wcs->equinox == equinox)
return;
/* Convert center from B1950 (FK4) to J2000 (FK5) */
if (equinox == 2000.0 && wcs->equinox == 1950.0) {
if (wcs->coorflip) {
fk425e (&wcs->crval[1], &wcs->crval[0], wcs->epoch);
wcs->cel.ref[1] = wcs->crval[0];
wcs->cel.ref[0] = wcs->crval[1];
}
else {
fk425e (&wcs->crval[0], &wcs->crval[1], wcs->epoch);
wcs->cel.ref[0] = wcs->crval[0];
wcs->cel.ref[1] = wcs->crval[1];
}
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
wcs->equinox = 2000.0;
strcpy (wcs->radecsys, "FK5");
wcs->syswcs = WCS_J2000;
wcs->cel.flag = 0;
wcs->wcsl.flag = 0;
}
/* Convert center from J2000 (FK5) to B1950 (FK4) */
else if (equinox == 1950.0 && wcs->equinox == 2000.0) {
if (wcs->coorflip) {
fk524e (&wcs->crval[1], &wcs->crval[0], wcs->epoch);
wcs->cel.ref[1] = wcs->crval[0];
wcs->cel.ref[0] = wcs->crval[1];
}
else {
fk524e (&wcs->crval[0], &wcs->crval[1], wcs->epoch);
wcs->cel.ref[0] = wcs->crval[0];
wcs->cel.ref[1] = wcs->crval[1];
}
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
wcs->equinox = 1950.0;
strcpy (wcs->radecsys, "FK4");
wcs->syswcs = WCS_B1950;
wcs->cel.flag = 0;
wcs->wcsl.flag = 0;
}
wcsoutinit (wcs, wcs->radecsys);
wcsininit (wcs, wcs->radecsys);
return;
}
/* Set scale and rotation in WCS structure */
void
wcscdset (wcs, cd)
struct WorldCoor *wcs; /* World coordinate system structure */
double *cd; /* CD matrix, ignored if NULL */
{
double tcd;
if (cd == NULL)
return;
wcs->rotmat = 1;
wcs->cd[0] = cd[0];
wcs->cd[1] = cd[1];
wcs->cd[2] = cd[2];
wcs->cd[3] = cd[3];
(void) matinv (2, wcs->cd, wcs->dc);
/* Compute scale */
wcs->xinc = sqrt (cd[0]*cd[0] + cd[2]*cd[2]);
wcs->yinc = sqrt (cd[1]*cd[1] + cd[3]*cd[3]);
/* Deal with x=Dec/y=RA case */
if (wcs->coorflip) {
tcd = cd[1];
cd[1] = -cd[2];
cd[2] = -tcd;
}
wcslibrot (wcs);
wcs->wcson = 1;
/* Compute image rotation */
wcsrotset (wcs);
wcs->cdelt[0] = wcs->xinc;
wcs->cdelt[1] = wcs->yinc;
return;
}
/* Set scale and rotation in WCS structure from axis scale and rotation */
void
wcsdeltset (wcs, cdelt1, cdelt2, crota)
struct WorldCoor *wcs; /* World coordinate system structure */
double cdelt1; /* degrees/pixel in first axis (or both axes) */
double cdelt2; /* degrees/pixel in second axis if nonzero */
double crota; /* Rotation counterclockwise in degrees */
{
double *pci;
double crot, srot;
int i, j, naxes;
naxes = wcs->naxis;
if (naxes > 2)
naxes = 2;
wcs->cdelt[0] = cdelt1;
if (cdelt2 != 0.0)
wcs->cdelt[1] = cdelt2;
else
wcs->cdelt[1] = cdelt1;
wcs->xinc = wcs->cdelt[0];
wcs->yinc = wcs->cdelt[1];
pci = wcs->pc;
for (i = 0; i < naxes; i++) {
for (j = 0; j < naxes; j++) {
if (i ==j)
*pci = 1.0;
else
*pci = 0.0;
pci++;
}
}
wcs->rotmat = 0;
/* If image is reversed, value of CROTA is flipped, too */
wcs->rot = crota;
if (wcs->rot < 0.0)
wcs->rot = wcs->rot + 360.0;
if (wcs->rot >= 360.0)
wcs->rot = wcs->rot - 360.0;
crot = cos (degrad(wcs->rot));
if (cdelt1 * cdelt2 > 0)
srot = sin (-degrad(wcs->rot));
else
srot = sin (degrad(wcs->rot));
/* Set CD matrix */
wcs->cd[0] = wcs->cdelt[0] * crot;
if (wcs->cdelt[0] < 0)
wcs->cd[1] = -fabs (wcs->cdelt[1]) * srot;
else
wcs->cd[1] = fabs (wcs->cdelt[1]) * srot;
if (wcs->cdelt[1] < 0)
wcs->cd[2] = fabs (wcs->cdelt[0]) * srot;
else
wcs->cd[2] = -fabs (wcs->cdelt[0]) * srot;
wcs->cd[3] = wcs->cdelt[1] * crot;
(void) matinv (2, wcs->cd, wcs->dc);
/* Set rotation matrix */
wcslibrot (wcs);
/* Set image rotation and mirroring */
if (wcs->coorflip) {
if (wcs->cdelt[0] < 0 && wcs->cdelt[1] > 0) {
wcs->imflip = 1;
wcs->imrot = wcs->rot - 90.0;
if (wcs->imrot < -180.0) wcs->imrot = wcs->imrot + 360.0;
wcs->pa_north = wcs->rot;
wcs->pa_east = wcs->rot - 90.0;
if (wcs->pa_east < -180.0) wcs->pa_east = wcs->pa_east + 360.0;
}
else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] < 0) {
wcs->imflip = 1;
wcs->imrot = wcs->rot + 90.0;
if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
wcs->pa_north = wcs->rot;
wcs->pa_east = wcs->rot - 90.0;
if (wcs->pa_east < -180.0) wcs->pa_east = wcs->pa_east + 360.0;
}
else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] > 0) {
wcs->imflip = 0;
wcs->imrot = wcs->rot + 90.0;
if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
wcs->pa_north = wcs->imrot;
wcs->pa_east = wcs->rot + 90.0;
if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
}
else if (wcs->cdelt[0] < 0 && wcs->cdelt[1] < 0) {
wcs->imflip = 0;
wcs->imrot = wcs->rot - 90.0;
if (wcs->imrot < -180.0) wcs->imrot = wcs->imrot + 360.0;
wcs->pa_north = wcs->imrot;
wcs->pa_east = wcs->rot + 90.0;
if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
}
}
else {
if (wcs->cdelt[0] < 0 && wcs->cdelt[1] > 0) {
wcs->imflip = 0;
wcs->imrot = wcs->rot;
wcs->pa_north = wcs->rot + 90.0;
if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
wcs->pa_east = wcs->rot + 180.0;
if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
}
else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] < 0) {
wcs->imflip = 0;
wcs->imrot = wcs->rot + 180.0;
if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
wcs->pa_north = wcs->imrot + 90.0;
if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
wcs->pa_east = wcs->imrot + 180.0;
if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
}
else if (wcs->cdelt[0] > 0 && wcs->cdelt[1] > 0) {
wcs->imflip = 1;
wcs->imrot = -wcs->rot;
wcs->pa_north = wcs->imrot + 90.0;
if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
wcs->pa_east = wcs->rot;
}
else if (wcs->cdelt[0] < 0 && wcs->cdelt[1] < 0) {
wcs->imflip = 1;
wcs->imrot = wcs->rot + 180.0;
if (wcs->imrot > 180.0) wcs->imrot = wcs->imrot - 360.0;
wcs->pa_north = wcs->imrot + 90.0;
if (wcs->pa_north > 180.0) wcs->pa_north = wcs->pa_north - 360.0;
wcs->pa_east = wcs->rot + 90.0;
if (wcs->pa_east > 180.0) wcs->pa_east = wcs->pa_east - 360.0;
}
}
return;
}
/* Set scale and rotation in WCS structure */
void
wcspcset (wcs, cdelt1, cdelt2, pc)
struct WorldCoor *wcs; /* World coordinate system structure */
double cdelt1; /* degrees/pixel in first axis (or both axes) */
double cdelt2; /* degrees/pixel in second axis if nonzero */
double *pc; /* Rotation matrix, ignored if NULL */
{
double *pci, *pc0i;
int i, j, naxes;
if (pc == NULL)
return;
naxes = wcs->naxis;
/* if (naxes > 2)
naxes = 2; */
if (naxes < 1 || naxes > 9) {
naxes = wcs->naxes;
wcs->naxis = naxes;
}
wcs->cdelt[0] = cdelt1;
if (cdelt2 != 0.0)
wcs->cdelt[1] = cdelt2;
else
wcs->cdelt[1] = cdelt1;
wcs->xinc = wcs->cdelt[0];
wcs->yinc = wcs->cdelt[1];
/* Set rotation matrix */
pci = wcs->pc;
pc0i = pc;
for (i = 0; i < naxes; i++) {
for (j = 0; j < naxes; j++) {
*pci = *pc0i;
pci++;
pc0i++;
}
}
/* Set CD matrix */
if (naxes > 1) {
wcs->cd[0] = pc[0] * wcs->cdelt[0];
wcs->cd[1] = pc[1] * wcs->cdelt[0];
wcs->cd[2] = pc[naxes] * wcs->cdelt[1];
wcs->cd[3] = pc[naxes+1] * wcs->cdelt[1];
}
else {
wcs->cd[0] = pc[0] * wcs->cdelt[0];
wcs->cd[1] = 0.0;
wcs->cd[2] = 0.0;
wcs->cd[3] = 1.0;
}
(void) matinv (2, wcs->cd, wcs->dc);
wcs->rotmat = 1;
(void)linset (&wcs->lin);
wcs->wcson = 1;
wcsrotset (wcs);
return;
}
/* Set up rotation matrix for WCSLIB projection subroutines */
static void
wcslibrot (wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
int i, mem, naxes;
naxes = wcs->naxis;
if (naxes > 2)
naxes = 2;
if (naxes < 1 || naxes > 9) {
naxes = wcs->naxes;
wcs->naxis = naxes;
}
mem = naxes * naxes * sizeof(double);
if (wcs->lin.piximg == NULL)
wcs->lin.piximg = (double*)malloc(mem);
if (wcs->lin.piximg != NULL) {
if (wcs->lin.imgpix == NULL)
wcs->lin.imgpix = (double*)malloc(mem);
if (wcs->lin.imgpix != NULL) {
wcs->lin.flag = LINSET;
if (naxes == 2) {
for (i = 0; i < 4; i++) {
wcs->lin.piximg[i] = wcs->cd[i];
}
}
else if (naxes == 3) {
for (i = 0; i < 9; i++)
wcs->lin.piximg[i] = 0.0;
wcs->lin.piximg[0] = wcs->cd[0];
wcs->lin.piximg[1] = wcs->cd[1];
wcs->lin.piximg[3] = wcs->cd[2];
wcs->lin.piximg[4] = wcs->cd[3];
wcs->lin.piximg[8] = 1.0;
}
else if (naxes == 4) {
for (i = 0; i < 16; i++)
wcs->lin.piximg[i] = 0.0;
wcs->lin.piximg[0] = wcs->cd[0];
wcs->lin.piximg[1] = wcs->cd[1];
wcs->lin.piximg[4] = wcs->cd[2];
wcs->lin.piximg[5] = wcs->cd[3];
wcs->lin.piximg[10] = 1.0;
wcs->lin.piximg[15] = 1.0;
}
(void) matinv (naxes, wcs->lin.piximg, wcs->lin.imgpix);
wcs->lin.crpix = wcs->crpix;
wcs->lin.cdelt = wcs->cdelt;
wcs->lin.pc = wcs->pc;
wcs->lin.flag = LINSET;
}
}
return;
}
/* Compute image rotation */
void
wcsrotset (wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
int off;
double cra, cdec, xc, xn, xe, yc, yn, ye;
/* If image is one-dimensional, leave rotation angle alone */
if (wcs->nxpix < 1.5 || wcs->nypix < 1.5) {
wcs->imrot = wcs->rot;
wcs->pa_north = wcs->rot + 90.0;
wcs->pa_east = wcs->rot + 180.0;
return;
}
/* Do not try anything if image is LINEAR (not Cartesian projection) */
if (wcs->syswcs == WCS_LINEAR)
return;
wcs->xinc = fabs (wcs->xinc);
wcs->yinc = fabs (wcs->yinc);
/* Compute position angles of North and East in image */
xc = wcs->xrefpix;
yc = wcs->yrefpix;
pix2wcs (wcs, xc, yc, &cra, &cdec);
if (wcs->coorflip) {
wcs2pix (wcs, cra+wcs->yinc, cdec, &xe, &ye, &off);
wcs2pix (wcs, cra, cdec+wcs->xinc, &xn, &yn, &off);
}
else {
wcs2pix (wcs, cra+wcs->xinc, cdec, &xe, &ye, &off);
wcs2pix (wcs, cra, cdec+wcs->yinc, &xn, &yn, &off);
}
wcs->pa_north = raddeg (atan2 (yn-yc, xn-xc));
if (wcs->pa_north < -90.0)
wcs->pa_north = wcs->pa_north + 360.0;
wcs->pa_east = raddeg (atan2 (ye-yc, xe-xc));
if (wcs->pa_east < -90.0)
wcs->pa_east = wcs->pa_east + 360.0;
/* Compute image rotation angle from North */
if (wcs->pa_north < -90.0)
wcs->imrot = 270.0 + wcs->pa_north;
else
wcs->imrot = wcs->pa_north - 90.0;
/* Compute CROTA */
if (wcs->coorflip) {
wcs->rot = wcs->imrot + 90.0;
if (wcs->rot < 0.0)
wcs->rot = wcs->rot + 360.0;
}
else
wcs->rot = wcs->imrot;
if (wcs->rot < 0.0)
wcs->rot = wcs->rot + 360.0;
if (wcs->rot >= 360.0)
wcs->rot = wcs->rot - 360.0;
/* Set image mirror flag based on axis orientation */
wcs->imflip = 0;
if (wcs->pa_east - wcs->pa_north < -80.0 &&
wcs->pa_east - wcs->pa_north > -100.0)
wcs->imflip = 1;
if (wcs->pa_east - wcs->pa_north < 280.0 &&
wcs->pa_east - wcs->pa_north > 260.0)
wcs->imflip = 1;
if (wcs->pa_north - wcs->pa_east > 80.0 &&
wcs->pa_north - wcs->pa_east < 100.0)
wcs->imflip = 1;
if (wcs->coorflip) {
if (wcs->imflip)
wcs->yinc = -wcs->yinc;
}
else {
if (!wcs->imflip)
wcs->xinc = -wcs->xinc;
}
return;
}
/* Return 1 if WCS structure is filled, else 0 */
int
iswcs (wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
if (wcs == NULL)
return (0);
else
return (wcs->wcson);
}
/* Return 0 if WCS structure is filled, else 1 */
int
nowcs (wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
if (wcs == NULL)
return (1);
else
return (!wcs->wcson);
}
/* Reset the center of a WCS structure */
void
wcsshift (wcs,rra,rdec,coorsys)
struct WorldCoor *wcs; /* World coordinate system structure */
double rra; /* Reference pixel right ascension in degrees */
double rdec; /* Reference pixel declination in degrees */
char *coorsys; /* FK4 or FK5 coordinates (1950 or 2000) */
{
if (nowcs (wcs))
return;
/* Approximate world coordinate system from a known plate scale */
wcs->crval[0] = rra;
wcs->crval[1] = rdec;
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
/* Coordinate reference frame */
strcpy (wcs->radecsys,coorsys);
wcs->syswcs = wcscsys (coorsys);
if (wcs->syswcs == WCS_B1950)
wcs->equinox = 1950.0;
else
wcs->equinox = 2000.0;
return;
}
/* Print position of WCS center, if WCS is set */
void
wcscent (wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
double xpix,ypix, xpos1, xpos2, ypos1, ypos2;
char wcstring[32];
double width, height, secpix, secpixh, secpixw;
int lstr = 32;
if (nowcs (wcs))
(void)fprintf (stderr,"No WCS information available\n");
else {
if (wcs->prjcode == WCS_DSS)
(void)fprintf (stderr,"WCS plate center %s\n", wcs->center);
xpix = 0.5 * wcs->nxpix;
ypix = 0.5 * wcs->nypix;
(void) pix2wcst (wcs,xpix,ypix,wcstring, lstr);
(void)fprintf (stderr,"WCS center %s %s %s %s at pixel (%.2f,%.2f)\n",
wcs->ctype[0],wcs->ctype[1],wcstring,wcs->ptype,xpix,ypix);
/* Image width */
(void) pix2wcs (wcs,1.0,ypix,&xpos1,&ypos1);
(void) pix2wcs (wcs,wcs->nxpix,ypix,&xpos2,&ypos2);
if (wcs->syswcs == WCS_LINEAR) {
width = xpos2 - xpos1;
if (width < 100.0)
(void)fprintf (stderr, "WCS width = %.5f %s ",width, wcs->units[0]);
else
(void)fprintf (stderr, "WCS width = %.3f %s ",width, wcs->units[0]);
}
else {
width = wcsdist (xpos1,ypos1,xpos2,ypos2);
if (width < 1/60.0)
(void)fprintf (stderr, "WCS width = %.2f arcsec ",width*3600.0);
else if (width < 1.0)
(void)fprintf (stderr, "WCS width = %.2f arcmin ",width*60.0);
else
(void)fprintf (stderr, "WCS width = %.3f degrees ",width);
}
secpixw = width / (wcs->nxpix - 1.0);
/* Image height */
(void) pix2wcs (wcs,xpix,1.0,&xpos1,&ypos1);
(void) pix2wcs (wcs,xpix,wcs->nypix,&xpos2,&ypos2);
if (wcs->syswcs == WCS_LINEAR) {
height = ypos2 - ypos1;
if (height < 100.0)
(void)fprintf (stderr, " height = %.5f %s ",height, wcs->units[1]);
else
(void)fprintf (stderr, " height = %.3f %s ",height, wcs->units[1]);
}
else {
height = wcsdist (xpos1,ypos1,xpos2,ypos2);
if (height < 1/60.0)
(void) fprintf (stderr, " height = %.2f arcsec",height*3600.0);
else if (height < 1.0)
(void) fprintf (stderr, " height = %.2f arcmin",height*60.0);
else
(void) fprintf (stderr, " height = %.3f degrees",height);
}
secpixh = height / (wcs->nypix - 1.0);
/* Image scale */
if (wcs->syswcs == WCS_LINEAR) {
(void) fprintf (stderr,"\n");
(void) fprintf (stderr,"WCS %.5f %s/pixel, %.5f %s/pixel\n",
wcs->xinc,wcs->units[0],wcs->yinc,wcs->units[1]);
}
else {
if (wcs->xinc != 0.0 && wcs->yinc != 0.0)
secpix = (fabs(wcs->xinc) + fabs(wcs->yinc)) * 0.5 * 3600.0;
else if (secpixh > 0.0 && secpixw > 0.0)
secpix = (secpixw + secpixh) * 0.5 * 3600.0;
else if (wcs->xinc != 0.0 || wcs->yinc != 0.0)
secpix = (fabs(wcs->xinc) + fabs(wcs->yinc)) * 3600.0;
else
secpix = (secpixw + secpixh) * 3600.0;
if (secpix < 100.0)
(void) fprintf (stderr, " %.3f arcsec/pixel\n",secpix);
else if (secpix < 3600.0)
(void) fprintf (stderr, " %.3f arcmin/pixel\n",secpix/60.0);
else
(void) fprintf (stderr, " %.3f degrees/pixel\n",secpix/3600.0);
}
}
return;
}
/* Return RA and Dec of image center, plus size in RA and Dec */
void
wcssize (wcs, cra, cdec, dra, ddec)
struct WorldCoor *wcs; /* World coordinate system structure */
double *cra; /* Right ascension of image center (deg) (returned) */
double *cdec; /* Declination of image center (deg) (returned) */
double *dra; /* Half-width in right ascension (deg) (returned) */
double *ddec; /* Half-width in declination (deg) (returned) */
{
double width, height;
/* Find right ascension and declination of coordinates */
if (iswcs(wcs)) {
wcsfull (wcs, cra, cdec, &width, &height);
*dra = 0.5 * width / cos (degrad (*cdec));
*ddec = 0.5 * height;
}
else {
*cra = 0.0;
*cdec = 0.0;
*dra = 0.0;
*ddec = 0.0;
}
return;
}
/* Return RA and Dec of image center, plus size in degrees */
void
wcsfull (wcs, cra, cdec, width, height)
struct WorldCoor *wcs; /* World coordinate system structure */
double *cra; /* Right ascension of image center (deg) (returned) */
double *cdec; /* Declination of image center (deg) (returned) */
double *width; /* Width in degrees (returned) */
double *height; /* Height in degrees (returned) */
{
double xpix, ypix, xpos1, xpos2, ypos1, ypos2, xcpix, ycpix;
double xcent, ycent;
/* Find right ascension and declination of coordinates */
if (iswcs(wcs)) {
xcpix = (0.5 * wcs->nxpix) + 0.5;
ycpix = (0.5 * wcs->nypix) + 0.5;
(void) pix2wcs (wcs,xcpix,ycpix,&xcent, &ycent);
*cra = xcent;
*cdec = ycent;
/* Compute image width in degrees */
xpix = 0.500001;
(void) pix2wcs (wcs,xpix,ycpix,&xpos1,&ypos1);
xpix = wcs->nxpix + 0.499999;
(void) pix2wcs (wcs,xpix,ycpix,&xpos2,&ypos2);
if (strncmp (wcs->ptype,"LINEAR",6) &&
strncmp (wcs->ptype,"PIXEL",5)) {
*width = wcsdist (xpos1,ypos1,xpos2,ypos2);
}
else
*width = sqrt (((ypos2-ypos1) * (ypos2-ypos1)) +
((xpos2-xpos1) * (xpos2-xpos1)));
/* Compute image height in degrees */
ypix = 0.5;
(void) pix2wcs (wcs,xcpix,ypix,&xpos1,&ypos1);
ypix = wcs->nypix + 0.5;
(void) pix2wcs (wcs,xcpix,ypix,&xpos2,&ypos2);
if (strncmp (wcs->ptype,"LINEAR",6) &&
strncmp (wcs->ptype,"PIXEL",5))
*height = wcsdist (xpos1,ypos1,xpos2,ypos2);
else
*height = sqrt (((ypos2-ypos1) * (ypos2-ypos1)) +
((xpos2-xpos1) * (xpos2-xpos1)));
}
else {
*cra = 0.0;
*cdec = 0.0;
*width = 0.0;
*height = 0.0;
}
return;
}
/* Return minimum and maximum RA and Dec of image in degrees */
void
wcsrange (wcs, ra1, ra2, dec1, dec2)
struct WorldCoor *wcs; /* World coordinate system structure */
double *ra1; /* Minimum right ascension of image (deg) (returned) */
double *ra2; /* Maximum right ascension of image (deg) (returned) */
double *dec1; /* Minimum declination of image (deg) (returned) */
double *dec2; /* Maximum declination of image (deg) (returned) */
{
double xpos1, xpos2, xpos3, xpos4, ypos1, ypos2, ypos3, ypos4, temp;
if (iswcs(wcs)) {
/* Compute image corner coordinates in degrees */
(void) pix2wcs (wcs,1.0,1.0,&xpos1,&ypos1);
(void) pix2wcs (wcs,1.0,wcs->nypix,&xpos2,&ypos2);
(void) pix2wcs (wcs,wcs->nxpix,1.0,&xpos3,&ypos3);
(void) pix2wcs (wcs,wcs->nxpix,wcs->nypix,&xpos4,&ypos4);
/* Find minimum right ascension or longitude */
*ra1 = xpos1;
if (xpos2 < *ra1) *ra1 = xpos2;
if (xpos3 < *ra1) *ra1 = xpos3;
if (xpos4 < *ra1) *ra1 = xpos4;
/* Find maximum right ascension or longitude */
*ra2 = xpos1;
if (xpos2 > *ra2) *ra2 = xpos2;
if (xpos3 > *ra2) *ra2 = xpos3;
if (xpos4 > *ra2) *ra2 = xpos4;
if (wcs->syswcs != WCS_LINEAR && wcs->syswcs != WCS_XY) {
if (*ra2 - *ra1 > 180.0) {
temp = *ra1;
*ra1 = *ra2;
*ra2 = temp;
}
}
/* Find minimum declination or latitude */
*dec1 = ypos1;
if (ypos2 < *dec1) *dec1 = ypos2;
if (ypos3 < *dec1) *dec1 = ypos3;
if (ypos4 < *dec1) *dec1 = ypos4;
/* Find maximum declination or latitude */
*dec2 = ypos1;
if (ypos2 > *dec2) *dec2 = ypos2;
if (ypos3 > *dec2) *dec2 = ypos3;
if (ypos4 > *dec2) *dec2 = ypos4;
}
else {
*ra1 = 0.0;
*ra2 = 0.0;
*dec1 = 0.0;
*dec2 = 0.0;
}
return;
}
/* Compute distance in degrees between two sky coordinates */
double
wcsdist (x1,y1,x2,y2)
double x1,y1; /* (RA,Dec) or (Long,Lat) in degrees */
double x2,y2; /* (RA,Dec) or (Long,Lat) in degrees */
{
double r, diffi;
double pos1[3], pos2[3], w, diff;
int i;
/* Convert two vectors to direction cosines */
r = 1.0;
d2v3 (x1, y1, r, pos1);
d2v3 (x2, y2, r, pos2);
/* Modulus squared of half the difference vector */
w = 0.0;
for (i = 0; i < 3; i++) {
diffi = pos1[i] - pos2[i];
w = w + (diffi * diffi);
}
w = w / 4.0;
if (w > 1.0) w = 1.0;
/* Angle beween the vectors */
diff = 2.0 * atan2 (sqrt (w), sqrt (1.0 - w));
diff = raddeg (diff);
return (diff);
}
/* Compute distance in degrees between two sky coordinates */
double
wcsdist1 (x1,y1,x2,y2)
double x1,y1; /* (RA,Dec) or (Long,Lat) in degrees */
double x2,y2; /* (RA,Dec) or (Long,Lat) in degrees */
{
double d1, d2, r;
double pos1[3], pos2[3], w, diff;
int i;
/* Convert two vectors to direction cosines */
r = 1.0;
d2v3 (x1, y1, r, pos1);
d2v3 (x2, y2, r, pos2);
w = 0.0;
d1 = 0.0;
d2 = 0.0;
for (i = 0; i < 3; i++) {
w = w + (pos1[i] * pos2[i]);
d1 = d1 + (pos1[i] * pos1[i]);
d2 = d2 + (pos2[i] * pos2[i]);
}
diff = acosdeg (w / (sqrt (d1) * sqrt (d2)));
return (diff);
}
/* Compute distance in degrees between two sky coordinates away from pole */
double
wcsdiff (x1,y1,x2,y2)
double x1,y1; /* (RA,Dec) or (Long,Lat) in degrees */
double x2,y2; /* (RA,Dec) or (Long,Lat) in degrees */
{
double xdiff, ydiff, ycos, diff;
ycos = cos (degrad ((y2 + y1) / 2.0));
xdiff = x2 - x1;
if (xdiff > 180.0)
xdiff = xdiff - 360.0;
if (xdiff < -180.0)
xdiff = xdiff + 360.0;
xdiff = xdiff / ycos;
ydiff = (y2 - y1);
diff = sqrt ((xdiff * xdiff) + (ydiff * ydiff));
return (diff);
}
/* Initialize catalog search command set by -wcscom */
void
wcscominit (wcs, i, command)
struct WorldCoor *wcs; /* World coordinate system structure */
int i; /* Number of command (0-9) to initialize */
char *command; /* command with %s where coordinates will go */
{
int lcom,icom;
if (iswcs(wcs)) {
lcom = strlen (command);
if (lcom > 0) {
if (wcs->command_format[i] != NULL)
free (wcs->command_format[i]);
wcs->command_format[i] = (char *) calloc (lcom+2, 1);
if (wcs->command_format[i] == NULL)
return;
for (icom = 0; icom < lcom; icom++) {
if (command[icom] == '_')
wcs->command_format[i][icom] = ' ';
else
wcs->command_format[i][icom] = command[icom];
}
wcs->command_format[i][lcom] = 0;
}
}
return;
}
/* Execute Unix command with world coordinates (from x,y) and/or filename */
void
wcscom ( wcs, i, filename, xfile, yfile, wcstring )
struct WorldCoor *wcs; /* World coordinate system structure */
int i; /* Number of command (0-9) to execute */
char *filename; /* Image file name */
double xfile,yfile; /* Image pixel coordinates for WCS command */
char *wcstring; /* WCS String from pix2wcst() */
{
char command[120];
char comform[120];
char xystring[32];
char *fileform, *posform, *imform;
int ier;
if (nowcs (wcs)) {
(void)fprintf(stderr,"WCSCOM: no WCS\n");
return;
}
if (wcs->command_format[i] != NULL)
strcpy (comform, wcs->command_format[i]);
else
strcpy (comform, "sgsc -ah %s");
if (comform[0] > 0) {
/* Create and execute search command */
fileform = strsrch (comform,"%f");
imform = strsrch (comform,"%x");
posform = strsrch (comform,"%s");
if (imform != NULL) {
*(imform+1) = 's';
(void)sprintf (xystring, "%.2f %.2f", xfile, yfile);
if (fileform != NULL) {
*(fileform+1) = 's';
if (posform == NULL) {
if (imform < fileform)
(void)sprintf(command, comform, xystring, filename);
else
(void)sprintf(command, comform, filename, xystring);
}
else if (fileform < posform) {
if (imform < fileform)
(void)sprintf(command, comform, xystring, filename,
wcstring);
else if (imform < posform)
(void)sprintf(command, comform, filename, xystring,
wcstring);
else
(void)sprintf(command, comform, filename, wcstring,
xystring);
}
else
if (imform < posform)
(void)sprintf(command, comform, xystring, wcstring,
filename);
else if (imform < fileform)
(void)sprintf(command, comform, wcstring, xystring,
filename);
else
(void)sprintf(command, comform, wcstring, filename,
xystring);
}
else if (posform == NULL)
(void)sprintf(command, comform, xystring);
else if (imform < posform)
(void)sprintf(command, comform, xystring, wcstring);
else
(void)sprintf(command, comform, wcstring, xystring);
}
else if (fileform != NULL) {
*(fileform+1) = 's';
if (posform == NULL)
(void)sprintf(command, comform, filename);
else if (fileform < posform)
(void)sprintf(command, comform, filename, wcstring);
else
(void)sprintf(command, comform, wcstring, filename);
}
else
(void)sprintf(command, comform, wcstring);
ier = system (command);
if (ier)
(void)fprintf(stderr,"WCSCOM: %s failed %d\n",command,ier);
}
return;
}
/* Initialize WCS output coordinate system for use by PIX2WCS() */
void
wcsoutinit (wcs, coorsys)
struct WorldCoor *wcs; /* World coordinate system structure */
char *coorsys; /* Input world coordinate system:
FK4, FK5, B1950, J2000, GALACTIC, ECLIPTIC
fk4, fk5, b1950, j2000, galactic, ecliptic */
{
int sysout, i;
if (nowcs (wcs))
return;
/* If argument is null, set to image system and equinox */
if (coorsys == NULL || strlen (coorsys) < 1 ||
!strcmp(coorsys,"IMSYS") || !strcmp(coorsys,"imsys")) {
sysout = wcs->syswcs;
strcpy (wcs->radecout, wcs->radecsys);
wcs->eqout = wcs->equinox;
if (sysout == WCS_B1950) {
if (wcs->eqout != 1950.0) {
wcs->radecout[0] = 'B';
sprintf (wcs->radecout+1,"%.4f", wcs->equinox);
i = strlen(wcs->radecout) - 1;
if (wcs->radecout[i] == '0')
wcs->radecout[i] = (char)0;
i = strlen(wcs->radecout) - 1;
if (wcs->radecout[i] == '0')
wcs->radecout[i] = (char)0;
i = strlen(wcs->radecout) - 1;
if (wcs->radecout[i] == '0')
wcs->radecout[i] = (char)0;
}
else
strcpy (wcs->radecout, "B1950");
}
else if (sysout == WCS_J2000) {
if (wcs->eqout != 2000.0) {
wcs->radecout[0] = 'J';
sprintf (wcs->radecout+1,"%.4f", wcs->equinox);
i = strlen(wcs->radecout) - 1;
if (wcs->radecout[i] == '0')
wcs->radecout[i] = (char)0;
i = strlen(wcs->radecout) - 1;
if (wcs->radecout[i] == '0')
wcs->radecout[i] = (char)0;
i = strlen(wcs->radecout) - 1;
if (wcs->radecout[i] == '0')
wcs->radecout[i] = (char)0;
}
else
strcpy (wcs->radecout, "J2000");
}
}
/* Ignore new coordinate system if it is not supported */
else {
if ((sysout = wcscsys (coorsys)) < 0)
return;
/* Do not try to convert linear or alt-az coordinates */
if (sysout != wcs->syswcs &&
(wcs->syswcs == WCS_LINEAR || wcs->syswcs == WCS_ALTAZ))
return;
strcpy (wcs->radecout, coorsys);
wcs->eqout = wcsceq (coorsys);
}
wcs->sysout = sysout;
if (wcs->wcson) {
/* Set output in degrees flag and number of decimal places */
if (wcs->sysout == WCS_GALACTIC || wcs->sysout == WCS_ECLIPTIC ||
wcs->sysout == WCS_PLANET) {
wcs->degout = 1;
wcs->ndec = 5;
}
else if (wcs->sysout == WCS_ALTAZ) {
wcs->degout = 1;
wcs->ndec = 5;
}
else if (wcs->sysout == WCS_NPOLE || wcs->sysout == WCS_SPA) {
wcs->degout = 1;
wcs->ndec = 5;
}
else {
wcs->degout = 0;
wcs->ndec = 3;
}
}
return;
}
/* Return current value of WCS output coordinate system set by -wcsout */
char *
getwcsout(wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
if (nowcs (wcs))
return (NULL);
else
return(wcs->radecout);
}
/* Initialize WCS input coordinate system for use by WCS2PIX() */
void
wcsininit (wcs, coorsys)
struct WorldCoor *wcs; /* World coordinate system structure */
char *coorsys; /* Input world coordinate system:
FK4, FK5, B1950, J2000, GALACTIC, ECLIPTIC
fk4, fk5, b1950, j2000, galactic, ecliptic */
{
int sysin, i;
if (nowcs (wcs))
return;
/* If argument is null, set to image system and equinox */
if (coorsys == NULL || strlen (coorsys) < 1) {
wcs->sysin = wcs->syswcs;
strcpy (wcs->radecin, wcs->radecsys);
wcs->eqin = wcs->equinox;
if (wcs->sysin == WCS_B1950) {
if (wcs->eqin != 1950.0) {
wcs->radecin[0] = 'B';
sprintf (wcs->radecin+1,"%.4f", wcs->equinox);
i = strlen(wcs->radecin) - 1;
if (wcs->radecin[i] == '0')
wcs->radecin[i] = (char)0;
i = strlen(wcs->radecin) - 1;
if (wcs->radecin[i] == '0')
wcs->radecin[i] = (char)0;
i = strlen(wcs->radecin) - 1;
if (wcs->radecin[i] == '0')
wcs->radecin[i] = (char)0;
}
else
strcpy (wcs->radecin, "B1950");
}
else if (wcs->sysin == WCS_J2000) {
if (wcs->eqin != 2000.0) {
wcs->radecin[0] = 'J';
sprintf (wcs->radecin+1,"%.4f", wcs->equinox);
i = strlen(wcs->radecin) - 1;
if (wcs->radecin[i] == '0')
wcs->radecin[i] = (char)0;
i = strlen(wcs->radecin) - 1;
if (wcs->radecin[i] == '0')
wcs->radecin[i] = (char)0;
i = strlen(wcs->radecin) - 1;
if (wcs->radecin[i] == '0')
wcs->radecin[i] = (char)0;
}
else
strcpy (wcs->radecin, "J2000");
}
}
/* Ignore new coordinate system if it is not supported */
if ((sysin = wcscsys (coorsys)) < 0)
return;
wcs->sysin = sysin;
wcs->eqin = wcsceq (coorsys);
strcpy (wcs->radecin, coorsys);
return;
}
/* Return current value of WCS input coordinate system set by wcsininit */
char *
getwcsin (wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
if (nowcs (wcs))
return (NULL);
else
return (wcs->radecin);
}
/* Set WCS output in degrees or hh:mm:ss dd:mm:ss, returning old flag value */
int
setwcsdeg(wcs, new)
struct WorldCoor *wcs; /* World coordinate system structure */
int new; /* 1: degrees, 0: h:m:s d:m:s */
{
int old;
if (nowcs (wcs))
return (0);
old = wcs->degout;
wcs->degout = new;
if (new == 1 && old == 0 && wcs->ndec == 3)
wcs->ndec = 6;
if (new == 0 && old == 1 && wcs->ndec == 5)
wcs->ndec = 3;
return(old);
}
/* Set number of decimal places in pix2wcst output string */
int
wcsndec (wcs, ndec)
struct WorldCoor *wcs; /* World coordinate system structure */
int ndec; /* Number of decimal places in output string */
/* If < 0, return current unchanged value */
{
if (nowcs (wcs))
return (0);
else if (ndec >= 0)
wcs->ndec = ndec;
return (wcs->ndec);
}
/* Return current value of coordinate system */
char *
getradecsys(wcs)
struct WorldCoor *wcs; /* World coordinate system structure */
{
if (nowcs (wcs))
return (NULL);
else
return (wcs->radecsys);
}
/* Set output string mode for LINEAR coordinates */
void
setwcslin (wcs, mode)
struct WorldCoor *wcs; /* World coordinate system structure */
int mode; /* mode = 0: x y linear
mode = 1: x units x units
mode = 2: x y linear units */
{
if (iswcs (wcs))
wcs->linmode = mode;
return;
}
/* Convert pixel coordinates to World Coordinate string */
int
pix2wcst (wcs, xpix, ypix, wcstring, lstr)
struct WorldCoor *wcs; /* World coordinate system structure */
double xpix,ypix; /* Image coordinates in pixels */
char *wcstring; /* World coordinate string (returned) */
int lstr; /* Length of world coordinate string (returned) */
{
double xpos,ypos;
char rastr[32], decstr[32];
int minlength, lunits, lstring;
if (nowcs (wcs)) {
if (lstr > 0)
wcstring[0] = 0;
return(0);
}
pix2wcs (wcs,xpix,ypix,&xpos,&ypos);
/* If point is off scale, set string accordingly */
if (wcs->offscl) {
(void)sprintf (wcstring,"Off map");
return (1);
}
/* Print coordinates in degrees */
else if (wcs->degout == 1) {
minlength = 9 + (2 * wcs->ndec);
if (lstr > minlength) {
deg2str (rastr, 32, xpos, wcs->ndec);
deg2str (decstr, 32, ypos, wcs->ndec);
if (wcs->tabsys)
(void)sprintf (wcstring,"%s %s", rastr, decstr);
else
(void)sprintf (wcstring,"%s %s", rastr, decstr);
lstr = lstr - minlength;
}
else {
if (wcs->tabsys)
strncpy (wcstring,"********* **********",lstr);
else
strncpy (wcstring,"*******************",lstr);
lstr = 0;
}
}
/* print coordinates in sexagesimal notation */
else if (wcs->degout == 0) {
minlength = 18 + (2 * wcs->ndec);
if (lstr > minlength) {
if (wcs->sysout == WCS_J2000 || wcs->sysout == WCS_B1950) {
ra2str (rastr, 32, xpos, wcs->ndec);
dec2str (decstr, 32, ypos, wcs->ndec-1);
}
else {
dec2str (rastr, 32, xpos, wcs->ndec);
dec2str (decstr, 32, ypos, wcs->ndec);
}
if (wcs->tabsys) {
(void)sprintf (wcstring,"%s %s", rastr, decstr);
}
else {
(void)sprintf (wcstring,"%s %s", rastr, decstr);
}
lstr = lstr - minlength;
}
else {
if (wcs->tabsys) {
strncpy (wcstring,"************* *************",lstr);
}
else {
strncpy (wcstring,"**************************",lstr);
}
lstr = 0;
}
}
/* Label galactic coordinates */
if (wcs->sysout == WCS_GALACTIC) {
if (lstr > 9 && wcs->printsys) {
if (wcs->tabsys)
strcat (wcstring," galactic");
else
strcat (wcstring," galactic");
}
}
/* Label ecliptic coordinates */
else if (wcs->sysout == WCS_ECLIPTIC) {
if (lstr > 9 && wcs->printsys) {
if (wcs->tabsys)
strcat (wcstring," ecliptic");
else
strcat (wcstring," ecliptic");
}
}
/* Label planet coordinates */
else if (wcs->sysout == WCS_PLANET) {
if (lstr > 9 && wcs->printsys) {
if (wcs->tabsys)
strcat (wcstring," planet");
else
strcat (wcstring," planet");
}
}
/* Label alt-az coordinates */
else if (wcs->sysout == WCS_ALTAZ) {
if (lstr > 7 && wcs->printsys) {
if (wcs->tabsys)
strcat (wcstring," alt-az");
else
strcat (wcstring," alt-az");
}
}
/* Label north pole angle coordinates */
else if (wcs->sysout == WCS_NPOLE) {
if (lstr > 7 && wcs->printsys) {
if (wcs->tabsys)
strcat (wcstring," long-npa");
else
strcat (wcstring," long-npa");
}
}
/* Label south pole angle coordinates */
else if (wcs->sysout == WCS_SPA) {
if (lstr > 7 && wcs->printsys) {
if (wcs->tabsys)
strcat (wcstring," long-spa");
else
strcat (wcstring," long-spa");
}
}
/* Label equatorial coordinates */
else if (wcs->sysout==WCS_B1950 || wcs->sysout==WCS_J2000) {
if (lstr > (int) strlen(wcs->radecout)+1 && wcs->printsys) {
if (wcs->tabsys)
strcat (wcstring," ");
else
strcat (wcstring," ");
strcat (wcstring, wcs->radecout);
}
}
/* Output linear coordinates */
else {
num2str (rastr, xpos, 0, wcs->ndec);
num2str (decstr, ypos, 0, wcs->ndec);
lstring = strlen (rastr) + strlen (decstr) + 1;
lunits = strlen (wcs->units[0]) + strlen (wcs->units[1]) + 2;
if (wcs->syswcs == WCS_LINEAR && wcs->linmode == 1) {
if (lstr > lstring + lunits) {
if (strlen (wcs->units[0]) > 0) {
strcat (rastr, " ");
strcat (rastr, wcs->units[0]);
}
if (strlen (wcs->units[1]) > 0) {
strcat (decstr, " ");
strcat (decstr, wcs->units[1]);
}
lstring = lstring + lunits;
}
}
if (lstr > lstring) {
if (wcs->tabsys)
(void)sprintf (wcstring,"%s %s", rastr, decstr);
else
(void)sprintf (wcstring,"%s %s", rastr, decstr);
}
else {
if (wcs->tabsys)
strncpy (wcstring,"********** *********",lstr);
else
strncpy (wcstring,"*******************",lstr);
}
if (wcs->syswcs == WCS_LINEAR && wcs->linmode != 1 &&
lstr > lstring + 7)
strcat (wcstring, " linear");
if (wcs->syswcs == WCS_LINEAR && wcs->linmode == 2 &&
lstr > lstring + lunits + 7) {
if (strlen (wcs->units[0]) > 0) {
strcat (wcstring, " ");
strcat (wcstring, wcs->units[0]);
}
if (strlen (wcs->units[1]) > 0) {
strcat (wcstring, " ");
strcat (wcstring, wcs->units[1]);
}
}
}
return (1);
}
/* Convert pixel coordinates to World Coordinates */
void
pix2wcs (wcs,xpix,ypix,xpos,ypos)
struct WorldCoor *wcs; /* World coordinate system structure */
double xpix,ypix; /* x and y image coordinates in pixels */
double *xpos,*ypos; /* RA and Dec in degrees (returned) */
{
double xpi, ypi, xp, yp;
double eqin, eqout;
int wcspos();
if (nowcs (wcs))
return;
wcs->xpix = xpix;
wcs->ypix = ypix;
wcs->zpix = zpix;
wcs->offscl = 0;
/* If this WCS is converted from another WCS rather than pixels, convert now */
if (wcs->wcs != NULL) {
pix2wcs (wcs->wcs, xpix, ypix, &xpi, &ypi);
}
else {
pix2foc (wcs, xpix, ypix, &xpi, &ypi);
}
/* Convert image coordinates to sky coordinates */
/* Use Digitized Sky Survey plate fit */
if (wcs->prjcode == WCS_DSS) {
if (dsspos (xpi, ypi, wcs, &xp, &yp))
wcs->offscl = 1;
}
/* Use SAO plate fit */
else if (wcs->prjcode == WCS_PLT) {
if (platepos (xpi, ypi, wcs, &xp, &yp))
wcs->offscl = 1;
}
/* Use NOAO IRAF corrected plane tangent projection */
else if (wcs->prjcode == WCS_TNX) {
if (tnxpos (xpi, ypi, wcs, &xp, &yp))
wcs->offscl = 1;
}
/* Use NOAO IRAF corrected zenithal projection */
else if (wcs->prjcode == WCS_ZPX) {
if (zpxpos (xpi, ypi, wcs, &xp, &yp))
wcs->offscl = 1;
}
/* Use Classic AIPS projections */
else if (wcs->wcsproj == WCS_OLD || wcs->prjcode <= 0) {
if (worldpos (xpi, ypi, wcs, &xp, &yp))
wcs->offscl = 1;
}
/* Use Mark Calabretta's WCSLIB projections */
else if (wcspos (xpi, ypi, wcs, &xp, &yp))
wcs->offscl = 1;
/* Do not change coordinates if offscale */
if (wcs->offscl) {
*xpos = 0.0;
*ypos = 0.0;
}
else {
/* Convert coordinates to output system, if not LINEAR */
if (wcs->prjcode > 0) {
/* Convert coordinates to desired output system */
eqin = wcs->equinox;
eqout = wcs->eqout;
wcscon (wcs->syswcs,wcs->sysout,eqin,eqout,&xp,&yp,wcs->epoch);
}
if (wcs->latbase == 90)
yp = 90.0 - yp;
else if (wcs->latbase == -90)
yp = yp - 90.0;
wcs->xpos = xp;
wcs->ypos = yp;
*xpos = xp;
*ypos = yp;
}
/* Keep RA/longitude within range if spherical coordinate output
(Not LINEAR or XY) */
if (wcs->sysout > 0 && wcs->sysout != 6 && wcs->sysout != 10) {
if (*xpos < 0.0)
*xpos = *xpos + 360.0;
else if (*xpos > 360.0)
*xpos = *xpos - 360.0;
}
return;
}
/* Convert World Coordinates to pixel coordinates */
void
wcs2pix (wcs, xpos, ypos, xpix, ypix, offscl)
struct WorldCoor *wcs; /* World coordinate system structure */
double xpos,ypos; /* World coordinates in degrees */
double *xpix,*ypix; /* Image coordinates in pixels */
int *offscl; /* 0 if within bounds, else off scale */
{
wcsc2pix (wcs, xpos, ypos, wcs->radecin, xpix, ypix, offscl);
return;
}
/* Convert World Coordinates to pixel coordinates */
void
wcsc2pix (wcs, xpos, ypos, coorsys, xpix, ypix, offscl)
struct WorldCoor *wcs; /* World coordinate system structure */
double xpos,ypos; /* World coordinates in degrees */
char *coorsys; /* Input world coordinate system:
FK4, FK5, B1950, J2000, GALACTIC, ECLIPTIC
fk4, fk5, b1950, j2000, galactic, ecliptic
* If NULL, use image WCS */
double *xpix,*ypix; /* Image coordinates in pixels */
int *offscl; /* 0 if within bounds, else off scale */
{
double xp, yp, xpi, ypi;
double eqin, eqout;
int sysin;
int wcspix();
if (nowcs (wcs))
return;
*offscl = 0;
xp = xpos;
yp = ypos;
if (wcs->latbase == 90)
yp = 90.0 - yp;
else if (wcs->latbase == -90)
yp = yp - 90.0;
if (coorsys == NULL) {
sysin = wcs->syswcs;
eqin = wcs->equinox;
}
else {
sysin = wcscsys (coorsys);
eqin = wcsceq (coorsys);
}
wcs->zpix = 1.0;
/* Convert coordinates to same system as image */
if (sysin > 0 && sysin != 6 && sysin != 10) {
eqout = wcs->equinox;
wcscon (sysin, wcs->syswcs, eqin, eqout, &xp, &yp, wcs->epoch);
}
/* Convert sky coordinates to image coordinates */
/* Use Digitized Sky Survey plate fit */
if (wcs->prjcode == WCS_DSS) {
if (dsspix (xp, yp, wcs, &xpi, &ypi))
*offscl = 1;
}
/* Use SAO polynomial plate fit */
else if (wcs->prjcode == WCS_PLT) {
if (platepix (xp, yp, wcs, &xpi, &ypi))
*offscl = 1;
}
/* Use NOAO IRAF corrected plane tangent projection */
else if (wcs->prjcode == WCS_TNX) {
if (tnxpix (xp, yp, wcs, &xpi, &ypi))
*offscl = 1;
}
/* Use NOAO IRAF corrected zenithal projection */
else if (wcs->prjcode == WCS_ZPX) {
if (zpxpix (xp, yp, wcs, &xpi, &ypi))
*offscl = 1;
}
/* Use Classic AIPS projections */
else if (wcs->wcsproj == WCS_OLD || wcs->prjcode <= 0) {
if (worldpix (xp, yp, wcs, &xpi, &ypi))
*offscl = 1;
}
/* Use Mark Calabretta's WCSLIB projections */
else if (wcspix (xp, yp, wcs, &xpi, &ypi)) {
*offscl = 1;
}
/* If this WCS is converted from another WCS rather than pixels, convert now */
if (wcs->wcs != NULL) {
wcsc2pix (wcs->wcs, xpi, ypi, NULL, xpix, ypix, offscl);
}
else {
foc2pix (wcs, xpi, ypi, xpix, ypix);
/* Set off-scale flag to 2 if off image but within bounds of projection */
if (!*offscl) {
if (*xpix < 0.5 || *ypix < 0.5)
*offscl = 2;
else if (*xpix > wcs->nxpix + 0.5 || *ypix > wcs->nypix + 0.5)
*offscl = 2;
}
}
wcs->offscl = *offscl;
wcs->xpos = xpos;
wcs->ypos = ypos;
wcs->xpix = *xpix;
wcs->ypix = *ypix;
return;
}
int
wcspos (xpix, ypix, wcs, xpos, ypos)
/* Input: */
double xpix; /* x pixel number (RA or long without rotation) */
double ypix; /* y pixel number (dec or lat without rotation) */
struct WorldCoor *wcs; /* WCS parameter structure */
/* Output: */
double *xpos; /* x (RA) coordinate (deg) */
double *ypos; /* y (dec) coordinate (deg) */
{
int offscl;
int i;
int wcsrev();
double wcscrd[4], imgcrd[4], pixcrd[4];
double phi, theta;
*xpos = 0.0;
*ypos = 0.0;
pixcrd[0] = xpix;
pixcrd[1] = ypix;
if (wcs->prjcode == WCS_CSC || wcs->prjcode == WCS_QSC ||
wcs->prjcode == WCS_TSC)
pixcrd[2] = (double) (izpix + 1);
else
pixcrd[2] = zpix;
pixcrd[3] = 1.0;
for (i = 0; i < 4; i++)
imgcrd[i] = 0.0;
offscl = wcsrev ((void *)&wcs->ctype, &wcs->wcsl, pixcrd, &wcs->lin, imgcrd,
&wcs->prj, &phi, &theta, wcs->crval, &wcs->cel, wcscrd);
if (offscl == 0) {
*xpos = wcscrd[wcs->wcsl.lng];
*ypos = wcscrd[wcs->wcsl.lat];
}
return (offscl);
}
int
wcspix (xpos, ypos, wcs, xpix, ypix)
/* Input: */
double xpos; /* x (RA) coordinate (deg) */
double ypos; /* y (dec) coordinate (deg) */
struct WorldCoor *wcs; /* WCS parameter structure */
/* Output: */
double *xpix; /* x pixel number (RA or long without rotation) */
double *ypix; /* y pixel number (dec or lat without rotation) */
{
int offscl;
int wcsfwd();
double wcscrd[4], imgcrd[4], pixcrd[4];
double phi, theta;
*xpix = 0.0;
*ypix = 0.0;
if (wcs->wcsl.flag != WCSSET) {
if (wcsset (wcs->lin.naxis, (void *)&wcs->ctype, &wcs->wcsl) )
return (1);
}
/* Set input for WCSLIB subroutines */
wcscrd[0] = 0.0;
wcscrd[1] = 0.0;
wcscrd[2] = 0.0;
wcscrd[3] = 0.0;
wcscrd[wcs->wcsl.lng] = xpos;
wcscrd[wcs->wcsl.lat] = ypos;
/* Initialize output for WCSLIB subroutines */
pixcrd[0] = 0.0;
pixcrd[1] = 0.0;
pixcrd[2] = 1.0;
pixcrd[3] = 1.0;
imgcrd[0] = 0.0;
imgcrd[1] = 0.0;
imgcrd[2] = 1.0;
imgcrd[3] = 1.0;
/* Invoke WCSLIB subroutines for coordinate conversion */
offscl = wcsfwd ((void *)&wcs->ctype, &wcs->wcsl, wcscrd, wcs->crval, &wcs->cel,
&phi, &theta, &wcs->prj, imgcrd, &wcs->lin, pixcrd);
if (!offscl) {
*xpix = pixcrd[0];
*ypix = pixcrd[1];
if (wcs->prjcode == WCS_CSC || wcs->prjcode == WCS_QSC ||
wcs->prjcode == WCS_TSC)
wcs->zpix = pixcrd[2] - 1.0;
else
wcs->zpix = pixcrd[2];
}
return (offscl);
}
/* Set third dimension for cube projections */
int
wcszin (izpix0)
int izpix0; /* coordinate in third dimension
(if < 0, return current value without changing it */
{
if (izpix0 > -1) {
izpix = izpix0;
zpix = (double) izpix0;
}
return (izpix);
}
/* Return third dimension for cube projections */
int
wcszout (wcs)
struct WorldCoor *wcs; /* WCS parameter structure */
{
return ((int) wcs->zpix);
}
/* Set file name for error messages */
void
setwcsfile (filename)
char *filename; /* FITS or IRAF file with WCS */
{ if (strlen (filename) < 256)
strcpy (wcsfile, filename);
else
strncpy (wcsfile, filename, 255);
return; }
/* Set error message */
void
setwcserr (errmsg)
char *errmsg; /* Error mesage < 80 char */
{ strcpy (wcserrmsg, errmsg); return; }
/* Print error message */
void
wcserr ()
{ if (strlen (wcsfile) > 0)
fprintf (stderr, "%s in file %s\n",wcserrmsg, wcsfile);
else
fprintf (stderr, "%s\n",wcserrmsg);
return; }
/* Flag to use AIPS WCS subroutines instead of WCSLIB */
void
setdefwcs (wp)
int wp;
{ wcsproj0 = wp; return; }
int
getdefwcs ()
{ return (wcsproj0); }
/* Save output default coordinate system */
static char wcscoor0[16];
void
savewcscoor (wcscoor)
char *wcscoor;
{ strcpy (wcscoor0, wcscoor); return; }
/* Return preset output default coordinate system */
char *
getwcscoor ()
{ return (wcscoor0); }
/* Save default commands */
static char *wcscom0[10];
void
savewcscom (i, wcscom)
int i;
char *wcscom;
{
int lcom;
if (i < 0) i = 0;
else if (i > 9) i = 9;
lcom = strlen (wcscom) + 2;
wcscom0[i] = (char *) calloc (lcom, 1);
if (wcscom0[i] != NULL)
strcpy (wcscom0[i], wcscom);
return;
}
void
setwcscom (wcs)
struct WorldCoor *wcs; /* WCS parameter structure */
{
char envar[16];
int i;
char *str;
if (nowcs(wcs))
return;
for (i = 0; i < 10; i++) {
if (i == 0)
strcpy (envar, "WCS_COMMAND");
else
sprintf (envar, "WCS_COMMAND%d", i);
if (wcscom0[i] != NULL)
wcscominit (wcs, i, wcscom0[i]);
else if ((str = getenv (envar)) != NULL)
wcscominit (wcs, i, str);
else if (i == 1)
wcscominit (wcs, i, "sua2 -ah %s"); /* F1= Search USNO-A2.0 Catalog */
else if (i == 2)
wcscominit (wcs, i, "sgsc -ah %s"); /* F2= Search HST GSC */
else if (i == 3)
wcscominit (wcs, i, "sty2 -ah %s"); /* F3= Search Tycho-2 Catalog */
else if (i == 4)
wcscominit (wcs, i, "sppm -ah %s"); /* F4= Search PPM Catalog */
else if (i == 5)
wcscominit (wcs, i, "ssao -ah %s"); /* F5= Search SAO Catalog */
else
wcs->command_format[i] = NULL;
}
return;
}
char *
getwcscom (i)
int i;
{ return (wcscom0[i]); }
void
freewcscom (wcs)
struct WorldCoor *wcs; /* WCS parameter structure */
{
int i;
for (i = 0; i < 10; i++) {
if (wcscom0[i] != NULL) {
free (wcscom0[i]);
wcscom0[i] = NULL;
}
}
if (iswcs (wcs)) {
for (i = 0; i < 10; i++) {
if (wcs->command_format[i] != NULL) {
free (wcs->command_format[i]);
}
}
}
return;
}
int
cpwcs (header, cwcs)
char **header; /* Pointer to start of FITS header */
char *cwcs; /* Keyword suffix character for output WCS */
{
double tnum;
int dkwd[100];
int i, maxnkwd, ikwd, nleft, lbuff, lhead, nkwd, nbytes;
int nkwdw;
char **kwd;
char *newhead, *oldhead;
char kwdc[16], keyword[16];
char tstr[80];
/* Allocate array of keywords to be transferred */
maxnkwd = 100;
kwd = (char **)calloc (maxnkwd, sizeof(char *));
for (ikwd = 0; ikwd < maxnkwd; ikwd++)
kwd[ikwd] = (char *) calloc (16, 1);
/* Make list of all possible keywords to be transferred */
nkwd = 0;
strcpy (kwd[++nkwd], "EPOCH");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "EQUINOX");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "RADECSYS");
dkwd[nkwd] = 0;
strcpy (kwd[++nkwd], "CTYPE1");
dkwd[nkwd] = 0;
strcpy (kwd[++nkwd], "CTYPE2");
dkwd[nkwd] = 0;
strcpy (kwd[++nkwd], "CRVAL1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CRVAL2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CDELT1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CDELT2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CRPIX1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CRPIX2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CROTA1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CROTA2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CD1_1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CD1_2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CD2_1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "CD2_2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC1_1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC1_2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC2_1");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC2_2");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC001001");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC001002");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC002001");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "PC002002");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "LATPOLE");
dkwd[nkwd] = 1;
strcpy (kwd[++nkwd], "LONPOLE");
dkwd[nkwd] = 1;
for (i = 1; i < 13; i++) {
sprintf (keyword,"CO1_%d", i);
strcpy (kwd[++nkwd], keyword);
dkwd[nkwd] = 1;
}
for (i = 1; i < 13; i++) {
sprintf (keyword,"CO2_%d", i);
strcpy (kwd[++nkwd], keyword);
dkwd[nkwd] = 1;
}
for (i = 0; i < 10; i++) {
sprintf (keyword,"PROJP%d", i);
strcpy (kwd[++nkwd], keyword);
dkwd[nkwd] = 1;
}
for (i = 0; i < MAXPV; i++) {
sprintf (keyword,"PV1_%d", i);
strcpy (kwd[++nkwd], keyword);
dkwd[nkwd] = 1;
}
for (i = 0; i < MAXPV; i++) {
sprintf (keyword,"PV2_%d", i);
strcpy (kwd[++nkwd], keyword);
dkwd[nkwd] = 1;
}
/* Allocate new header buffer if needed */
lhead = (ksearch (*header, "END") - *header) + 80;
lbuff = gethlength (*header);
nleft = (lbuff - lhead) / 80;
if (nleft < nkwd) {
nbytes = lhead + (nkwd * 80) + 400;
newhead = (char *) calloc (1, nbytes);
strncpy (newhead, *header, lhead);
oldhead = *header;
header = &newhead;
free (oldhead);
}
/* Copy keywords to new WCS ID in header */
nkwdw = 0;
for (i = 0; i < nkwd; i++) {
if (dkwd[i]) {
if (hgetr8 (*header, kwd[i], &tnum)) {
nkwdw++;
if (!strncmp (kwd[i], "PC0", 3)) {
if (!strcmp (kwd[i], "PC001001"))
strcpy (kwdc, "PC1_1");
else if (!strcmp (kwd[i], "PC001002"))
strcpy (kwdc, "PC1_2");
else if (!strcmp (kwd[i], "PC002001"))
strcpy (kwdc, "PC2_1");
else
strcpy (kwdc, "PC2_2");
}
else
strcpy (kwdc, kwd[i]);
strncat (kwdc, cwcs, 1);
(void)hputr8 (*header, kwdc, tnum);
}
}
else {
if (hgets (*header, kwd[i], 80, tstr)) {
nkwdw++;
if (!strncmp (kwd[i], "RADECSYS", 8))
strcpy (kwdc, "RADECSY");
else
strcpy (kwdc, kwd[i]);
strncat (kwdc, cwcs, 1);
hputs (*header, kwdc, tstr);
}
}
}
/* Free keyword list array */
for (ikwd = 0; ikwd < maxnkwd; ikwd++)
free (kwd[ikwd]);
free (kwd);
kwd = NULL;
return (nkwdw);
}
/* Oct 28 1994 new program
* Dec 21 1994 Implement CD rotation matrix
* Dec 22 1994 Allow RA and DEC to be either x,y or y,x
*
* Mar 6 1995 Add Digital Sky Survey plate fit
* May 2 1995 Add prototype of PIX2WCST to WCSCOM
* May 25 1995 Print leading zero for hours and degrees
* Jun 21 1995 Add WCS2PIX to get pixels from WCS
* Jun 21 1995 Read plate scale from FITS header for plate solution
* Jul 6 1995 Pass WCS structure as argument; malloc it in WCSINIT
* Jul 6 1995 Check string lengths in PIX2WCST
* Aug 16 1995 Add galactic coordinate conversion to PIX2WCST
* Aug 17 1995 Return 0 from iswcs if wcs structure is not yet set
* Sep 8 1995 Do not include malloc.h if VMS
* Sep 8 1995 Check for legal WCS before trying anything
* Sep 8 1995 Do not try to set WCS if missing key keywords
* Oct 18 1995 Add WCSCENT and WCSDIST to print center and size of image
* Nov 6 1995 Include stdlib.h instead of malloc.h
* Dec 6 1995 Fix format statement in PIX2WCST
* Dec 19 1995 Change MALLOC to CALLOC to initialize array to zeroes
* Dec 19 1995 Explicitly initialize rotation matrix and yinc
* Dec 22 1995 If SECPIX is set, use approximate WCS
* Dec 22 1995 Always print coordinate system
*
* Jan 12 1996 Use plane-tangent, not linear, projection if SECPIX is set
* Jan 12 1996 Add WCSSET to set WCS without an image
* Feb 15 1996 Replace all calls to HGETC with HGETS
* Feb 20 1996 Add tab table output from PIX2WCST
* Apr 2 1996 Convert all equinoxes to B1950 or J2000
* Apr 26 1996 Get and use image epoch for accurate FK4/FK5 conversions
* May 16 1996 Clean up internal documentation
* May 17 1996 Return width in right ascension degrees, not sky degrees
* May 24 1996 Remove extraneous print command from WCSSIZE
* May 28 1996 Add NOWCS and WCSSHIFT subroutines
* Jun 11 1996 Drop unused variables after running lint
* Jun 12 1996 Set equinox as well as system in WCSSHIFT
* Jun 14 1996 Make DSS keyword searches more robust
* Jul 1 1996 Allow for SECPIX1 and SECPIX2 keywords
* Jul 2 1996 Test for CTYPE1 instead of CRVAL1
* Jul 5 1996 Declare all subroutines in wcs.h
* Jul 19 1996 Add subroutine WCSFULL to return real image size
* Aug 12 1996 Allow systemless coordinates which cannot be converted
* Aug 15 1996 Allow LINEAR WCS to pass numbers through transparently
* Aug 15 1996 Add WCSERR to print error message under calling program control
* Aug 16 1996 Add latitude and longitude as image coordinate types
* Aug 26 1996 Fix arguments to HLENGTH in WCSNINIT
* Aug 28 1996 Explicitly set OFFSCL in WCS2PIX if coordinates outside image
* Sep 3 1996 Return computed pixel values even if they are offscale
* Sep 6 1996 Allow filename to be passed by WCSCOM
* Oct 8 1996 Default to 2000 for EQUINOX and EPOCH and FK5 for RADECSYS
* Oct 8 1996 If EPOCH is 0 and EQUINOX is not set, default to 1950 and FK4
* Oct 15 1996 Add comparison when testing an assignment
* Oct 16 1996 Allow PIXEL CTYPE which means WCS is same as image coordinates
* Oct 21 1996 Add WCS_COMMAND environment variable
* Oct 25 1996 Add image scale to WCSCENT
* Oct 30 1996 Fix bugs in WCS2PIX
* Oct 31 1996 Fix CD matrix rotation angle computation
* Oct 31 1996 Use inline degree <-> radian conversion functions
* Nov 1 1996 Add option to change number of decimal places in PIX2WCST
* Nov 5 1996 Set wcs->crot to 1 if rotation matrix is used
* Dec 2 1996 Add altitide/azimuth coordinates
* Dec 13 1996 Fix search format setting from environment
*
* Jan 22 1997 Add ifdef for Eric Mandel (SAOtng)
* Feb 5 1997 Add wcsout for Eric Mandel
* Mar 20 1997 Drop unused variable STR in WCSCOM
* May 21 1997 Do not make pixel coordinates mod 360 in PIX2WCST
* May 22 1997 Add PIXEL prjcode = -1;
* Jul 11 1997 Get center pixel x and y from header even if no WCS
* Aug 7 1997 Add NOAO PIXSCALi keywords for default WCS
* Oct 15 1997 Do not reset reference pixel in WCSSHIFT
* Oct 20 1997 Set chip rotation
* Oct 24 1997 Keep longitudes between 0 and 360, not -180 and +180
* Nov 5 1997 Do no compute crot and srot; they are now computed in worldpos
* Dec 16 1997 Set rotation and axis increments from CD matrix
*
* Jan 6 1998 Deal with J2000 and B1950 as EQUINOX values (from ST)
* Jan 7 1998 Read INSTRUME and DETECTOR header keywords
* Jan 7 1998 Fix tab-separated output
* Jan 9 1998 Precess coordinates if either FITS projection or *DSS plate*
* Jan 16 1998 Change PTYPE to not include initial hyphen
* Jan 16 1998 Change WCSSET to WCSXINIT to avoid conflict with Calabretta
* Jan 23 1998 Change PCODE to PRJCODE to avoid conflict with Calabretta
* Jan 27 1998 Add LATPOLE and LONGPOLE for Calabretta projections
* Feb 5 1998 Make cd and dc into vectors; use matinv() to invert cd
* Feb 5 1998 In wcsoutinit(), check that corsys is a valid pointer
* Feb 18 1998 Fix bugs for Calabretta projections
* Feb 19 1998 Add wcs structure access subroutines from Eric Mandel
* Feb 19 1998 Add wcsreset() to make sure derived values are reset
* Feb 19 1998 Always set oldwcs to 1 if NCP projection
* Feb 19 1998 Add subroutine to set oldwcs default
* Feb 20 1998 Initialize projection types one at a time for SunOS C
* Feb 23 1998 Add TNX projection from NOAO; treat it as TAN
* Feb 23 1998 Compute size based on max and min coordinates, not sides
* Feb 26 1998 Add code to set center pixel if part of detector array
* Mar 6 1998 Write 8-character values to RADECSYS
* Mar 9 1998 Add naxis to WCS structure
* Mar 16 1998 Use separate subroutine for IRAF TNX projection
* Mar 20 1998 Set PC matrix if more than two axes and it's not in header
* Mar 20 1998 Reset lin flag in WCSRESET if CDELTn
* Mar 20 1998 Set CD matrix with CDELTs and CROTA in wcsinit and wcsreset
* Mar 20 1998 Allow initialization of rotation angle alone
* Mar 23 1998 Use dsspos() and dsspix() instead of platepos() and platepix()
* Mar 24 1998 Set up PLT/PLATE plate polynomial fit using platepos() and platepix()
* Mar 25 1998 Read plate fit coefficients from header in getwcscoeffs()
* Mar 27 1998 Check for FITS WCS before DSS WCS
* Mar 27 1998 Compute scale from edges if xinc and yinc not set in wcscent()
* Apr 6 1998 Change plate coefficient keywords from PLTij to COi_j
* Apr 6 1998 Read all coefficients in line instead of with subroutine
* Apr 7 1998 Change amd_i_coeff to i_coeff
* Apr 8 1998 Add wcseqset to change equinox after wcs has been set
* Apr 10 1998 Use separate counters for x and y polynomial coefficients
* Apr 13 1998 Use CD/CDELT+CDROTA if oldwcs is set
* Apr 14 1998 Use codes instead of strings for various coordinate systems
* Apr 14 1998 Separate input coordinate conversion from output conversion
* Apr 14 1998 Use wcscon() for most coordinate conversion
* Apr 17 1998 Always compute cdelt[]
* Apr 17 1998 Deal with reversed axis more correctly
* Apr 17 1998 Compute rotation angle and approximate CDELTn for polynomial
* Apr 23 1998 Deprecate xref/yref in favor of crval[]
* Apr 23 1998 Deprecate xrefpix/yrefpix in favor of crpix[]
* Apr 23 1998 Add LINEAR to coordinate system types
* Apr 23 1998 Always use AIPS subroutines for LINEAR or PIXEL
* Apr 24 1998 Format linear coordinates better
* Apr 28 1998 Change coordinate system flags to WCS_*
* Apr 28 1998 Change projection flags to WCS_*
* Apr 28 1998 Add subroutine wcsc2pix for coordinates to pixels with system
* Apr 28 1998 Add setlinmode() to set output string mode for LINEAR coordinates
* Apr 30 1998 Fix bug by setting degree flag for lat and long in wcsinit()
* Apr 30 1998 Allow leading "-"s in projecting in wcsxinit()
* May 1 1998 Assign new output coordinate system only if legitimate system
* May 1 1998 Do not allow oldwcs=1 unless allowed projection
* May 4 1998 Fix bug in units reading for LINEAR coordinates
* May 6 1998 Initialize to no CD matrix
* May 6 1998 Use TAN instead of TNX if oldwcs
* May 12 1998 Set 3rd and 4th coordinates in wcspos()
* May 12 1998 Return *xpos and *ypos = 0 in pix2wcs() if offscale
* May 12 1998 Declare undeclared external subroutines after lint
* May 13 1998 Add equinox conversion to specified output equinox
* May 13 1998 Set output or input system to image with null argument
* May 15 1998 Return reference pixel, cdelts, and rotation for DSS
* May 20 1998 Fix bad bug so setlinmode() is no-op if wcs not set
* May 20 1998 Fix bug so getwcsout() returns null pointer if wcs not set
* May 27 1998 Change WCS_LPR back to WCS_LIN; allow CAR in oldwcs
* May 28 1998 Go back to old WCSFULL, computing height and width from center
* May 29 1998 Add wcskinit() to initialize WCS from arguments
* May 29 1998 Add wcstype() to set projection from arguments
* May 29 1998 Add wcscdset(), and wcsdeltset() to set scale from arguments
* Jun 1 1998 In wcstype(), reconstruct ctype for WCS structure
* Jun 11 1998 Split off header-dependent subroutines to wcsinit.c
* Jun 18 1998 Add wcspcset() for PC matrix initialization
* Jun 24 1998 Add string lengths to ra2str(), dec2str, and deg2str() calls
* Jun 25 1998 Use AIPS software for CAR projection
* Jun 25 1998 Add wcsndec to set number of decimal places in output string
* Jul 6 1998 Add wcszin() and wcszout() to use third dimension of images
* Jul 7 1998 Change setlinmode() to setwcslin(); setdegout() to setwcsdeg()
* Jul 10 1998 Initialize matrices correctly for naxis > 2 in wcs<>set()
* Jul 16 1998 Initialize coordinates to be returned in wcspos()
* Jul 17 1998 Link lin structure arrays to wcs structure arrays
* Jul 20 1998 In wcscdset() compute sign of xinc and yinc from CD1_1, CD 2_2
* Jul 20 1998 In wcscdset() compute sign of rotation based on CD1_1, CD 1_2
* Jul 22 1998 Add wcslibrot() to compute lin() rotation matrix
* Jul 30 1998 Set wcs->naxes and lin.naxis in wcsxinit() and wcskinit()
* Aug 5 1998 Use old WCS subroutines to deal with COE projection (for ESO)
* Aug 14 1998 Add option to print image coordinates with wcscom()
* Aug 14 1998 Add multiple command options to wcscom*()
* Aug 31 1998 Declare undeclared arguments to wcspcset()
* Sep 3 1998 Set CD rotation and cdelts from sky axis position angles
* Sep 16 1998 Add option to use North Polar Angle instead of Latitude
* Sep 29 1998 Initialize additional WCS commands from the environment
* Oct 14 1998 Fix bug in wcssize() which didn't divide dra by cos(dec)
* Nov 12 1998 Fix sign of CROTA when either axis is reflected
* Dec 2 1998 Fix non-arcsecond scale factors in wcscent()
* Dec 2 1998 Add PLANET coordinate system to pix2wcst()
* Jan 20 1999 Free lin.imgpix and lin.piximg in wcsfree()
* Feb 22 1999 Fix bug setting latitude reference value of latbase != 0
* Feb 22 1999 Fix bug so that quad cube faces are 0-5, not 1-6
* Mar 16 1999 Always initialize all 4 imgcrds and pixcrds in wcspix()
* Mar 16 1999 Always return (0,0) from wcs2pix if offscale
* Apr 7 1999 Add code to put file name in error messages
* Apr 7 1999 Document utility subroutines at end of file
* May 6 1999 Fix bug printing height of LINEAR image
* Jun 16 1999 Add wcsrange() to return image RA and Dec limits
* Jul 8 1999 Always use FK5 and FK4 instead of J2000 and B1950 in RADECSYS
* Aug 16 1999 Print dd:mm:ss dd:mm:ss if not J2000 or B1950 output
* Aug 20 1999 Add WCS string argument to wcscom(); don't compute it there
* Aug 20 1999 Change F3 WCS command default from Tycho to ACT
* Oct 15 1999 Free wcs using wcsfree()
* Oct 21 1999 Drop declarations of unused functions after lint
* Oct 25 1999 Drop out of wcsfree() if wcs is null pointer
* Nov 17 1999 Fix bug which caused software to miss NCP projection
*
* Jan 24 2000 Default to AIPS for NCP, CAR, and COE proj.; if -z use WCSLIB
* Feb 24 2000 If coorsys is null in wcsc2pix, wcs->radecin is assumed
* May 10 2000 In wcstype(), default to WCS_LIN, not error (after Bill Joye)
* Jun 22 2000 In wcsrotset(), leave rotation angle alone in 1-d image
* Jul 3 2000 Initialize wcscrd[] to zero in wcspix()
*
* Feb 20 2001 Add recursion to wcs2pix() and pix2wcs() for dependent WCS's
* Mar 20 2001 Add braces to avoid ambiguity in if/else groupings
* Mar 22 2001 Free WCS structure in wcsfree even if it is not filled
* Sep 12 2001 Fix bug which omitted tab in pix2wcst() galactic coord output
*
* Mar 7 2002 Fix bug which gave wrong pa's and rotation for reflected RA
* (but correct WCS conversions!)
* Mar 28 2002 Add SZP projection
* Apr 3 2002 Synchronize projection types with other subroutines
* Apr 3 2002 Drop special cases of projections
* Apr 9 2002 Implement inversion of multiple WCSs in wcsc2pix()
* Apr 25 2002 Use Tycho-2 catalog instead of ACT in setwcscom()
* May 13 2002 Free WCSNAME in wcsfree()
*
* Mar 31 2003 Add distcode to wcstype()
* Apr 1 2003 Add calls to foc2pix() in wcs2pix() and pix2foc() in pix2wcs()
* May 20 2003 Declare argument i in savewcscom()
* Sep 29 2003 Fix bug to compute width and height correctly in wcsfull()
* Sep 29 2003 Fix bug to deal with all-sky images orrectly in wcsfull()
* Oct 1 2003 Rename wcs->naxes to wcs->naxis to match WCSLIB 3.2
* Nov 3 2003 Set distortion code by calling setdistcode() in wcstype()
* Dec 3 2003 Add back wcs->naxes for compatibility
* Dec 3 2003 Add braces in if...else in pix2wcst()
*
* Sep 17 2004 If spherical coordinate output, keep 0 < long/RA < 360
* Sep 17 2004 Fix bug in wcsfull() when wrapping around RA=0:00
* Nov 1 2004 Keep wcs->rot between 0 and 360
*
* Mar 9 2005 Fix bug in wcsrotset() which set rot > 360 to 360
* Jun 27 2005 Fix ctype in calls to wcs subroutines
* Jul 21 2005 Fix bug in wcsrange() at RA ~= 0.0
*
* Apr 24 2006 Always set inverse CD matrix to 2 dimensions in wcspcset()
* May 3 2006 (void *) pointers so arguments match type, from Robert Lupton
* Jun 30 2006 Set only 2-dimensional PC matrix; that is all lin* can deal with
* Oct 30 2006 In pix2wcs(), do not limit x to between 0 and 360 if XY or LINEAR
* Oct 30 2006 In wcsc2pix(), do not precess LINEAR or XY coordinates
* Dec 21 2006 Add cpwcs() to copy WCS keywords to new suffix
*
* Jan 4 2007 Fix pointer to header in cpwcs()
* Jan 5 2007 Drop declarations of wcscon(); it is in wcs.h
* Jan 9 2006 Drop declarations of fk425e() and fk524e(); moved to wcs.h
* Jan 9 2006 Drop *pix() and *pos() external declarations; moved to wcs.h
* Jan 9 2006 Drop matinv() external declaration; it is already in wcslib.h
* Feb 15 2007 If CTYPEi contains DET, set to WCS_PIX projection
* Feb 23 2007 Fix bug when checking for "DET" in CTYPEi
* Apr 2 2007 Fix PC to CD matrix conversion
* Jul 25 2007 Compute distance between two coordinates using d2v3()
*
* Apr 7 2010 In wcstype() set number of WCS projections from NWCSTYPE
*
* Mar 11 2011 Add NOAO ZPX projection (Frank Valdes)
* Mar 14 2011 Delete j<=MAXPV PVi_j parameters (for SCAMP polynomials via Ed Los)
* Mar 17 2011 Fix WCSDEP bug found by Ed Los
* May 9 2011 Free WCS structure recursively if WCSDEP is used
* Sep 1 2011 Add TPV projection type for SCAMP TAN with PVs
*
* Oct 19 2012 Drop d1 and d2 from wcsdist(); diffi from wcsdist1()
* Oct 19 2012 Drop depwcs; it's in main wcs structure
*/