/*** File libwcs/wcsinit.c
*** December 1, 2003
*** By Doug Mink, dmink@cfa.harvard.edu
*** Harvard-Smithsonian Center for Astrophysics
*** Copyright (C) 1998-2003
*** 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: dmink@cfa.harvard.edu
Postal address: Doug Mink
Smithsonian Astrophysical Observatory
60 Garden St.
Cambridge, MA 02138 USA
* Module: wcsinit.c (World Coordinate Systems)
* Purpose: Convert FITS WCS to pixels and vice versa:
* Subroutine: wcsinit (hstring) sets a WCS structure from an image header
* Subroutine: wcsninit (hstring,lh) sets a WCS structure from fixed-length header
* Subroutine: wcsinitn (hstring, name) sets a WCS structure for specified WCS
* Subroutine: wcsninitn (hstring,lh, name) sets a WCS structure for specified WCS
* Subroutine: wcsinitc (hstring, mchar) sets a WCS structure if multiple
* Subroutine: wcsninitc (hstring,lh,mchar) sets a WCS structure if multiple
* Subroutine: wcschar (hstring, name) returns suffix for specifed WCS
* Subroutine: wcseq (hstring, wcs) set radecsys and equinox from image header
* Subroutine: wcseqm (hstring, wcs, mchar) set radecsys and equinox if multiple
*/
#include <string.h> /* strstr, NULL */
#include <stdio.h> /* stderr */
#include <math.h>
#include "wcs.h"
#ifndef VMS
#include <stdlib.h>
#endif
static void wcseq();
static void wcseqm();
static void wcsioset();
void wcsrotset();
char wcserrmsg[80];
char wcschar();
struct WorldCoor *wcsinitn();
struct WorldCoor *wcsinitc();
/* set up a WCS structure from a FITS image header lhstring bytes long
* for a specified WCS name */
struct WorldCoor *
wcsninitn (hstring, lhstring, name)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
int lhstring; /* Length of FITS header in bytes */
char *name; /* character string with identifying name of WCS */
{
hlength (hstring, lhstring);
return (wcsinitn (hstring, name));
}
/* set up a WCS structure from a FITS image header for specified WCSNAME */
struct WorldCoor *
wcsinitn (hstring, name)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
char *name; /* character string with identifying name of WCS */
{
char mchar; /* Suffix character for one of multiple WCS */
mchar = wcschar (hstring, name);
if (mchar == '_') {
fprintf (stderr, "WCSINITN: WCS name %s not matched in FITS header\n",
name);
return (NULL);
}
return (wcsinitc (hstring, mchar));
}
/* WCSCHAR -- Find the letter for a specific WCS conversion */
char
wcschar (hstring, name)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
char *name; /* Name of WCS conversion to be matched
(case-independent) */
{
char *upname, *uppercase();
char cwcs, charwcs;
int iwcs;
char keyword[12];
char *upval, value[72];
/* If no WCS character, return 0 */
if (name == NULL)
return ((char) 0);
/* Convert input name to upper case */
upname = uppercase (name);
/* If single character name, return that character */
if (strlen (upname) == 1)
return (upname[0]);
/* Try to match input name to available WCSNAME names in header */
strcpy (keyword, "WCSNAME");
keyword[8] = (char) 0;
charwcs = '_';
for (iwcs = 0; iwcs < 27; iwcs++) {
if (iwcs > 0)
cwcs = (char) (64 + iwcs);
else
cwcs = (char) 0;
keyword[7] = cwcs;
if (hgets (hstring, keyword, 72, value)) {
upval = uppercase (value);
if (!strcmp (upval, upname))
charwcs = cwcs;
free (upval);
}
}
free (upname);
return (charwcs);
}
/* Make string of arbitrary case all uppercase */
char *
uppercase (string)
char *string;
{
int lstring, i;
char *upstring;
lstring = strlen (string);
upstring = (char *) calloc (1,lstring+1);
for (i = 0; i < lstring; i++) {
if (string[i] > 96 && string[i] < 123)
upstring[i] = string[i] - 32;
else
upstring[i] = string[i];
}
upstring[lstring] = (char) 0;
return (upstring);
}
/* set up a WCS structure from a FITS image header lhstring bytes long */
struct WorldCoor *
wcsninit (hstring, lhstring)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
int lhstring; /* Length of FITS header in bytes */
{
char mchar; /* Suffix character for one of multiple WCS */
mchar = (char) 0;
hlength (hstring, lhstring);
return (wcsinitc (hstring, mchar));
}
/* set up a WCS structure from a FITS image header lhstring bytes long */
struct WorldCoor *
wcsninitc (hstring, lhstring, mchar)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
int lhstring; /* Length of FITS header in bytes */
char mchar; /* Suffix character for one of multiple WCS */
{
hlength (hstring, lhstring);
if (mchar == ' ')
mchar = (char) 0;
return (wcsinitc (hstring, mchar));
}
/* set up a WCS structure from a FITS image header */
struct WorldCoor *
wcsinit (hstring)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
{
char mchar; /* Suffix character for one of multiple WCS */
mchar = (char) 0;
return (wcsinitc (hstring, mchar));
}
/* set up a WCS structure from a FITS image header for specified suffix */
struct WorldCoor *
wcsinitc (hstring, mchar)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
char mchar; /* Suffix character for one of multiple WCS */
{
struct WorldCoor *wcs, *depwcs;
char ctype1[32], ctype2[32];
char pvkey1[8],pvkey2[8],pvkey3[8];
char *hcoeff; /* pointer to first coeff's in header */
char decsign;
double rah,ram,ras, dsign,decd,decm,decs;
double dec_deg,ra_hours, secpix, ra0, ra1, dec0, dec1, cvel;
double cdelt1, cdelt2, cd[4], pc[16];
char keyword[16];
int ieq, i, j, m, naxes, cd11p, cd12p, cd21p, cd22p;
int ilat; /* coordinate for latitude or declination */
/*
int ix1, ix2, iy1, iy2, idx1, idx2, idy1, idy2;
double dxrefpix, dyrefpix;
char temp[48];
*/
char wcsname[64]; /* Name of WCS depended on by current WCS */
double mjd;
double rot;
int twod;
int iszpx = 0;
extern int tnxinit();
extern int platepos();
extern int dsspos();
wcs = (struct WorldCoor *) calloc (1, sizeof(struct WorldCoor));
/* Set WCS character and name in structure */
if (mchar == ' ')
mchar = (char) 0;
wcs->wcschar = mchar;
if (hgetsc (hstring, "WCSNAME",mchar, 63, wcsname)) {
wcs->wcsname = (char *) calloc (strlen (wcsname)+2, 1);
strcpy (wcs->wcsname, wcsname);
}
/* Set WCSLIB flags so that structures will be reinitialized */
wcs->cel.flag = 0;
wcs->lin.flag = 0;
wcs->wcsl.flag = 0;
/* Initialize to no plate fit */
wcs->ncoeff1 = 0;
wcs->ncoeff2 = 0;
/* Initialize to no CD matrix */
cdelt1 = 0.0;
cdelt2 = 0.0;
cd[0] = 0.0;
cd[1] = 0.0;
cd[2] = 0.0;
cd[3] = 0.0;
pc[0] = 0.0;
wcs->rotmat = 0;
wcs->rot = 0.0;
/* Header parameters independent of projection */
naxes = 0;
hgeti4c (hstring, "WCSAXES", mchar, &naxes);
if (naxes == 0)
hgeti4 (hstring, "WCSAXES", &naxes);
if (naxes == 0)
hgeti4 (hstring, "NAXIS", &naxes);
if (naxes == 0)
hgeti4 (hstring, "WCSDIM", &naxes);
if (naxes < 1) {
setwcserr ("WCSINIT: No WCSAXES, NAXIS, or WCSDIM keyword");
wcsfree (wcs);
return (NULL);
}
wcs->naxis = naxes;
wcs->lin.naxis = naxes;
wcs->nxpix = 0;
hgetr8 (hstring, "NAXIS1", &wcs->nxpix);
if (wcs->nxpix < 1) {
setwcserr ("WCSINIT: No NAXIS1 keyword");
wcsfree (wcs);
return (NULL);
}
wcs->nypix = 0;
hgetr8 (hstring, "NAXIS2", &wcs->nypix);
hgets (hstring, "INSTRUME", 16, wcs->instrument);
hgeti4 (hstring, "DETECTOR", &wcs->detector);
wcs->wcsproj = getdefwcs();
for (i = 0; i < 16; i++) wcs->pc[i] = 0.0;
for (i = 0; i < naxes; i++) wcs->pc[(i*naxes)+i] = 1.0;
/* If the current world coordinate system depends on another, set it now */
if (hgetsc (hstring, "WCSDEP",mchar, 63, wcsname)) {
if ((wcs->wcs = wcsinitn (hstring, wcsname)) == NULL) {
setwcserr ("WCSINIT: depended on WCS could not be set");
wcsfree (wcs);
return (NULL);
}
depwcs = wcs->wcs;
depwcs->wcsdep = wcs;
}
else
wcs->wcs = NULL;
/* Read radial velocity from image header */
wcs->radvel = 0.0;
wcs->zvel = 0.0;
cvel = 299792.5;
if (hgetr8c (hstring, "VSOURCE", mchar, &wcs->radvel))
wcs->zvel = wcs->radvel / cvel;
else if (hgetr8c (hstring, "ZSOURCE", mchar, &wcs->zvel))
wcs->radvel = wcs->zvel * cvel;
else if (hgetr8 (hstring, "VELOCITY", &wcs->radvel))
wcs->zvel = wcs->radvel / cvel;
/* World coordinate system reference coordinate information */
if (hgetsc (hstring, "CTYPE1", mchar, 16, ctype1)) {
if (!strncmp (ctype1+5,"ZPX", 3)) {
iszpx = 1;
ctype1[7] = 'N';
}
else
iszpx = 0;
/* Read second coordinate type */
strcpy (ctype2, ctype1);
if (!hgetsc (hstring, "CTYPE2", mchar, 16, ctype2))
twod = 0;
else {
twod = 1;
if (!strncmp (ctype2+5,"ZPX", 3)) {
iszpx = 1;
ctype2[7] = 'N';
}
}
strcpy (wcs->ctype[0], ctype1);
strcpy (wcs->ctype[1], ctype2);
if (strsrch (ctype2, "LAT") || strsrch (ctype2, "DEC"))
ilat = 2;
else
ilat = 1;
/* Read third and fourth coordinate types, if present */
strcpy (wcs->ctype[2], "");
hgetsc (hstring, "CTYPE3", mchar, 9, wcs->ctype[2]);
strcpy (wcs->ctype[3], "");
hgetsc (hstring, "CTYPE4", mchar, 9, wcs->ctype[3]);
/* Set projection type in WCS data structure */
if (wcstype (wcs, ctype1, ctype2)) {
wcsfree (wcs);
return (NULL);
}
/* Get units, if present, for linear coordinates */
if (wcs->prjcode == WCS_LIN) {
if (!hgetsc (hstring, "CUNIT1", mchar, 16, wcs->units[0])) {
if (!mgetstr (hstring, "WAT1", "units", 16, wcs->units[0])) {
wcs->units[0][0] = 0;
}
}
if (twod) {
if (!hgetsc (hstring, "CUNIT2", mchar, 16, wcs->units[1])) {
if (!mgetstr (hstring, "WAT2", "units", 16, wcs->units[1])) {
wcs->units[1][0] = 0;
}
}
}
}
/* Reference pixel coordinates and WCS value */
wcs->crpix[0] = 1.0;
hgetr8c (hstring, "CRPIX1", mchar, &wcs->crpix[0]);
wcs->crpix[1] = 1.0;
hgetr8c (hstring, "CRPIX2", mchar, &wcs->crpix[1]);
wcs->xrefpix = wcs->crpix[0];
wcs->yrefpix = wcs->crpix[1];
wcs->crval[0] = 0.0;
sprintf (keyword,"CRVAL1%c",mchar);
hgetr8c (hstring, "CRVAL1", mchar, &wcs->crval[0]);
wcs->crval[1] = 0.0;
hgetr8c (hstring, "CRVAL2", mchar, &wcs->crval[1]);
if (wcs->syswcs == WCS_NPOLE)
wcs->crval[1] = 90.0 - wcs->crval[1];
if (wcs->syswcs == WCS_SPA)
wcs->crval[1] = wcs->crval[1] - 90.0;
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
if (wcs->coorflip) {
wcs->cel.ref[0] = wcs->crval[1];
wcs->cel.ref[1] = wcs->crval[0];
}
else {
wcs->cel.ref[0] = wcs->crval[0];
wcs->cel.ref[1] = wcs->crval[1];
}
wcs->longpole = 999.0;
hgetr8c (hstring, "LONPOLE", mchar, &wcs->longpole);
wcs->cel.ref[2] = wcs->longpole;
wcs->latpole = 999.0;
hgetr8c (hstring, "LATPOLE", mchar, &wcs->latpole);
wcs->cel.ref[3] = wcs->latpole;
wcs->lin.crpix = wcs->crpix;
wcs->lin.cdelt = wcs->cdelt;
wcs->lin.pc = wcs->pc;
/* Projection constants (this should be projection-dependent */
wcs->prj.r0 = 0.0;
hgetr8c (hstring, "PROJR0", mchar, &wcs->prj.r0);
/* FITS WCS interim proposal projection constants */
for (i = 0; i < 10; i++) {
wcs->prj.p[i] = 0.0;
sprintf (keyword,"PROJP%d",i);
hgetr8c (hstring, keyword, mchar, &wcs->prj.p[i]);
}
sprintf (pvkey1, "PV%d_1", ilat);
sprintf (pvkey2, "PV%d_2", ilat);
sprintf (pvkey3, "PV%d_3", ilat);
/* FITS WCS standard projection constants (projection-dependent) */
if (wcs->prjcode == WCS_AZP || wcs->prjcode == WCS_SIN ||
wcs->prjcode == WCS_COP || wcs->prjcode == WCS_COE ||
wcs->prjcode == WCS_COD || wcs->prjcode == WCS_COO) {
hgetr8c (hstring, pvkey1, mchar, &wcs->prj.p[1]);
hgetr8c (hstring, pvkey2, mchar, &wcs->prj.p[2]);
}
else if (wcs->prjcode == WCS_SZP) {
hgetr8c (hstring, pvkey1, mchar, &wcs->prj.p[1]);
hgetr8c (hstring, pvkey2, mchar, &wcs->prj.p[2]);
if (wcs->prj.p[3] == 0.0)
wcs->prj.p[3] = 90.0;
hgetr8c (hstring, pvkey3, mchar, &wcs->prj.p[3]);
}
else if (wcs->prjcode == WCS_CEA) {
if (wcs->prj.p[1] == 0.0)
wcs->prj.p[1] = 1.0;
hgetr8c (hstring, pvkey1, mchar, &wcs->prj.p[1]);
}
else if (wcs->prjcode == WCS_CYP) {
if (wcs->prj.p[1] == 0.0)
wcs->prj.p[1] = 1.0;
hgetr8c (hstring, pvkey1, mchar, &wcs->prj.p[1]);
if (wcs->prj.p[2] == 0.0)
wcs->prj.p[2] = 1.0;
hgetr8c (hstring, pvkey2, mchar, &wcs->prj.p[2]);
}
else if (wcs->prjcode == WCS_AIR) {
if (wcs->prj.p[1] == 0.0)
wcs->prj.p[1] = 90.0;
hgetr8c (hstring, pvkey1, mchar, &wcs->prj.p[1]);
}
else if (wcs->prjcode == WCS_BON) {
hgetr8c (hstring, pvkey1, mchar, &wcs->prj.p[1]);
}
else if (wcs->prjcode == WCS_ZPN) {
for (i = 0; i < 10; i++) {
sprintf (keyword,"PV%d_%d", ilat, i);
hgetr8c (hstring, keyword, mchar, &wcs->prj.p[i]);
}
}
/* If ZPX, read coefficients from WATi keyword */
if (iszpx) {
char mkey[8];
sprintf (mkey,"WAT%d", ilat);
for (i = 0; i < 10; i++) {
wcs->prj.p[i] = 0.0;
sprintf (keyword,"projp%d",i);
mgetr8 (hstring, mkey, keyword, &wcs->prj.p[i]);
}
}
/* Coordinate reference frame, equinox, and epoch */
if (strncmp (wcs->ptype,"LINEAR",6) &&
strncmp (wcs->ptype,"PIXEL",5))
wcseqm (hstring, wcs, mchar);
wcsioset (wcs);
/* Read distortion coefficients, if present */
distortinit (wcs, hstring);
/* Use polynomial fit instead of projection, if present */
wcs->ncoeff1 = 0;
wcs->ncoeff2 = 0;
cd11p = hgetr8c (hstring, "CD1_1", mchar, &cd[0]);
cd12p = hgetr8c (hstring, "CD1_2", mchar, &cd[1]);
cd21p = hgetr8c (hstring, "CD2_1", mchar, &cd[2]);
cd22p = hgetr8c (hstring, "CD2_2", mchar, &cd[3]);
if (wcs->wcsproj != WCS_OLD &&
(hcoeff = ksearch (hstring,"CO1_1")) != NULL) {
wcs->prjcode = WCS_PLT;
(void)strcpy (wcs->ptype, "PLATE");
for (i = 0; i < 20; i++) {
sprintf (keyword,"CO1_%d", i+1);
wcs->x_coeff[i] = 0.0;
if (hgetr8 (hcoeff, keyword, &wcs->x_coeff[i]))
wcs->ncoeff1 = i + 1;
}
hcoeff = ksearch (hstring,"CO2_1");
for (i = 0; i < 20; i++) {
sprintf (keyword,"CO2_%d",i+1);
wcs->y_coeff[i] = 0.0;
if (hgetr8 (hcoeff, keyword, &wcs->y_coeff[i]))
wcs->ncoeff2 = i + 1;
}
/* Compute a nominal scale factor */
platepos (wcs->crpix[0], wcs->crpix[1], wcs, &ra0, &dec0);
platepos (wcs->crpix[0], wcs->crpix[1]+1.0, wcs, &ra1, &dec1);
wcs->yinc = dec1 - dec0;
wcs->xinc = -wcs->yinc;
/* Compute image rotation angle */
wcs->wcson = 1;
wcsrotset (wcs);
rot = degrad (wcs->rot);
/* Compute scale at reference pixel */
platepos (wcs->crpix[0], wcs->crpix[1], wcs, &ra0, &dec0);
platepos (wcs->crpix[0]+cos(rot),
wcs->crpix[1]+sin(rot), wcs, &ra1, &dec1);
wcs->cdelt[0] = -wcsdist (ra0, dec0, ra1, dec1);
wcs->xinc = wcs->cdelt[0];
platepos (wcs->crpix[0]+sin(rot),
wcs->crpix[1]+cos(rot), wcs, &ra1, &dec1);
wcs->cdelt[1] = wcsdist (ra0, dec0, ra1, dec1);
wcs->yinc = wcs->cdelt[1];
/* Set CD matrix from header */
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);
}
/* Else use CD matrix, if present */
else if (cd11p || cd12p || cd21p || cd22p) {
wcs->rotmat = 1;
wcscdset (wcs, cd);
}
/* Else get scaling from CDELT1 and CDELT2 */
else if (hgetr8c (hstring, "CDELT1", mchar, &cdelt1) != 0) {
hgetr8c (hstring, "CDELT2", mchar, &cdelt2);
/* If CDELT1 or CDELT2 is 0 or missing */
if (cdelt1 == 0.0 || (wcs->nypix > 1 && cdelt2 == 0.0)) {
if (ksearch (hstring,"SECPIX") != NULL ||
ksearch (hstring,"PIXSCALE") != NULL ||
ksearch (hstring,"PIXSCAL1") != NULL ||
ksearch (hstring,"XPIXSIZE") != NULL ||
ksearch (hstring,"SECPIX1") != NULL) {
secpix = 0.0;
hgetr8 (hstring,"SECPIX",&secpix);
if (secpix == 0.0)
hgetr8 (hstring,"PIXSCALE",&secpix);
if (secpix == 0.0) {
hgetr8 (hstring,"SECPIX1",&secpix);
if (secpix != 0.0) {
if (cdelt1 == 0.0)
cdelt1 = -secpix / 3600.0;
if (cdelt2 == 0.0) {
hgetr8 (hstring,"SECPIX2",&secpix);
cdelt2 = secpix / 3600.0;
}
}
else {
hgetr8 (hstring,"XPIXSIZE",&secpix);
if (secpix != 0.0) {
if (cdelt1 == 0.0)
cdelt1 = -secpix / 3600.0;
if (cdelt2 == 0.0) {
hgetr8 (hstring,"YPIXSIZE",&secpix);
cdelt2 = secpix / 3600.0;
}
}
else {
hgetr8 (hstring,"PIXSCAL1",&secpix);
if (secpix != 0.0 && cdelt1 == 0.0)
cdelt1 = -secpix / 3600.0;
if (cdelt2 == 0.0) {
hgetr8 (hstring,"PIXSCAL2",&secpix);
cdelt2 = secpix / 3600.0;
}
}
}
}
else {
if (cdelt1 == 0.0)
cdelt1 = -secpix / 3600.0;
if (cdelt2 == 0.0)
cdelt2 = secpix / 3600.0;
}
}
}
if (cdelt2 == 0.0 && wcs->nypix > 1)
cdelt2 = -cdelt1;
wcs->cdelt[2] = 1.0;
wcs->cdelt[3] = 1.0;
/* Use rotation matrix, if present */
for (i = 0; i < 16; i++)
wcs->pc[i] = 0.0;
/* FITS WCS interim rotation matrix */
if (!mchar && hgetr8 (hstring,"PC001001",&pc[0]) != 0) {
hgetr8 (hstring,"PC001002",&pc[1]);
if (naxes < 3) {
hgetr8 (hstring,"PC002001",&pc[2]);
pc[3] = wcs->pc[0];
hgetr8 (hstring,"PC002002",&pc[3]);
}
if (naxes == 3) {
hgetr8 (hstring,"PC001003",&pc[2]);
hgetr8 (hstring,"PC002001",&pc[3]);
pc[4] = wcs->pc[0];
hgetr8 (hstring,"PC002002",&pc[4]);
hgetr8 (hstring,"PC002003",&pc[5]);
hgetr8 (hstring,"PC003001",&pc[6]);
hgetr8 (hstring,"PC003002",&pc[7]);
pc[8] = 1.0;
hgetr8 (hstring,"PC003003",&pc[8]);
}
if (naxes > 3) {
hgetr8 (hstring,"PC001003",&pc[2]);
hgetr8 (hstring,"PC001004",&pc[3]);
hgetr8 (hstring,"PC002001",&pc[4]);
pc[5] = wcs->pc[0];
hgetr8 (hstring,"PC002002",&pc[5]);
hgetr8 (hstring,"PC002003",&pc[6]);
hgetr8 (hstring,"PC002004",&pc[7]);
hgetr8 (hstring,"PC003001",&pc[8]);
hgetr8 (hstring,"PC003002",&pc[9]);
pc[10] = 1.0;
hgetr8 (hstring,"PC003003",&pc[10]);
hgetr8 (hstring,"PC003004",&pc[11]);
hgetr8 (hstring,"PC004001",&pc[12]);
hgetr8 (hstring,"PC004002",&pc[13]);
hgetr8 (hstring,"PC004003",&pc[14]);
pc[15] = 1.0;
hgetr8 (hstring,"PC004004",&pc[15]);
}
wcspcset (wcs, cdelt1, cdelt2, pc);
}
/* FITS WCS standard rotation matrix */
else if (hgetr8c (hstring, "PC1_1", mchar, &pc[0]) != 0) {
hgetr8c (hstring, "PC1_2", mchar, &pc[1]);
if (naxes < 3) {
hgetr8c (hstring, "PC2_1", mchar, &pc[2]);
pc[3] = wcs->pc[0];
hgetr8c (hstring, "PC2_2", mchar, &pc[3]);
}
if (naxes == 3) {
hgetr8c (hstring, "PC1_3", mchar, &pc[2]);
hgetr8c (hstring, "PC2_1", mchar, &pc[3]);
pc[4] = wcs->pc[0];
hgetr8c (hstring, "PC2_2", mchar, &pc[4]);
hgetr8c (hstring, "PC2_3", mchar, &pc[5]);
hgetr8c (hstring, "PC3_1", mchar, &pc[6]);
hgetr8c (hstring, "PC3_2", mchar, &pc[7]);
pc[8] = 1.0;
hgetr8c (hstring, "PC3_3", mchar, &pc[8]);
}
if (naxes > 3) {
hgetr8c (hstring, "PC1_3", mchar, &pc[2]);
hgetr8c (hstring, "PC1_4", mchar, &pc[3]);
hgetr8c (hstring, "PC2_1", mchar, &pc[4]);
pc[5] = wcs->pc[0];
hgetr8c (hstring, "PC2_2", mchar, &pc[5]);
hgetr8c (hstring, "PC2_3", mchar, &pc[6]);
hgetr8c (hstring, "PC2_4", mchar, &pc[7]);
hgetr8c (hstring, "PC3_1", mchar, &pc[8]);
hgetr8c (hstring, "PC3_2", mchar, &pc[9]);
pc[10] = 1.0;
hgetr8c (hstring, "PC3_3", mchar, &pc[10]);
hgetr8c (hstring, "PC3_4", mchar, &pc[11]);
hgetr8c (hstring, "PC4_1", mchar, &pc[12]);
hgetr8c (hstring, "PC4_2", mchar, &pc[13]);
hgetr8c (hstring, "PC4_3", mchar, &pc[14]);
pc[15] = 1.0;
hgetr8c (hstring, "PC4_4", mchar, &pc[15]);
}
wcspcset (wcs, cdelt1, cdelt2, pc);
}
/* Otherwise, use CROTAn */
else {
rot = 0.0;
hgetr8c (hstring, "CROTA2", mchar, &rot);
if (rot == 0.) {
hgetr8c (hstring,"CROTA1", mchar, &rot);
}
wcsdeltset (wcs, cdelt1, cdelt2, rot);
}
}
/* If no scaling is present, set to 1 per pixel, no rotation */
else {
wcs->xinc = 1.0;
wcs->yinc = 1.0;
wcs->cdelt[0] = 1.0;
wcs->cdelt[1] = 1.0;
wcs->rot = 0.0;
wcs->rotmat = 0;
setwcserr ("WCSINIT: setting CDELT to 1");
}
/* Initialize TNX, defaulting to TAN if there is a problem */
if (wcs->prjcode == WCS_TNX) {
if (tnxinit (hstring, wcs)) {
wcs->ctype[0][6] = 'A';
wcs->ctype[0][7] = 'N';
wcs->ctype[1][6] = 'A';
wcs->ctype[1][7] = 'N';
wcs->prjcode = WCS_TAN;
}
}
/* If linear or pixel WCS, print "degrees" */
if (!strncmp (wcs->ptype,"LINEAR",6) ||
!strncmp (wcs->ptype,"PIXEL",5)) {
wcs->degout = -1;
wcs->ndec = 5;
}
wcs->wcson = 1;
}
else if (mchar != (char) 0) {
setwcserr ("WCSINITC: No image scale for WCS %c", mchar);
wcsfree (wcs);
return (NULL);
}
/* Plate solution coefficients */
else if (ksearch (hstring,"PLTRAH") != NULL) {
wcs->prjcode = WCS_DSS;
hcoeff = ksearch (hstring,"PLTRAH");
hgetr8 (hcoeff,"PLTRAH",&rah);
hgetr8 (hcoeff,"PLTRAM",&ram);
hgetr8 (hcoeff,"PLTRAS",&ras);
ra_hours = rah + (ram / (double)60.0) + (ras / (double)3600.0);
wcs->plate_ra = hrrad (ra_hours);
decsign = '+';
hgets (hcoeff,"PLTDECSN", 1, &decsign);
if (decsign == '-')
dsign = -1.;
else
dsign = 1.;
hgetr8 (hcoeff,"PLTDECD",&decd);
hgetr8 (hcoeff,"PLTDECM",&decm);
hgetr8 (hcoeff,"PLTDECS",&decs);
dec_deg = dsign * (decd+(decm/(double)60.0)+(decs/(double)3600.0));
wcs->plate_dec = degrad (dec_deg);
hgetr8 (hstring,"EQUINOX",&wcs->equinox);
hgeti4 (hstring,"EQUINOX",&ieq);
if (ieq == 1950)
strcpy (wcs->radecsys,"FK4");
else
strcpy (wcs->radecsys,"FK5");
wcs->epoch = wcs->equinox;
hgetr8 (hstring,"EPOCH",&wcs->epoch);
(void)sprintf (wcs->center,"%2.0f:%2.0f:%5.3f %c%2.0f:%2.0f:%5.3f %s",
rah,ram,ras,decsign,decd,decm,decs,wcs->radecsys);
hgetr8 (hstring,"PLTSCALE",&wcs->plate_scale);
hgetr8 (hstring,"XPIXELSZ",&wcs->x_pixel_size);
hgetr8 (hstring,"YPIXELSZ",&wcs->y_pixel_size);
hgetr8 (hstring,"CNPIX1",&wcs->x_pixel_offset);
hgetr8 (hstring,"CNPIX2",&wcs->y_pixel_offset);
hcoeff = ksearch (hstring,"PPO1");
for (i = 0; i < 6; i++) {
sprintf (keyword,"PPO%d", i+1);
wcs->ppo_coeff[i] = 0.0;
hgetr8 (hcoeff,keyword,&wcs->ppo_coeff[i]);
}
hcoeff = ksearch (hstring,"AMDX1");
for (i = 0; i < 20; i++) {
sprintf (keyword,"AMDX%d", i+1);
wcs->x_coeff[i] = 0.0;
hgetr8 (hcoeff, keyword, &wcs->x_coeff[i]);
}
hcoeff = ksearch (hstring,"AMDY1");
for (i = 0; i < 20; i++) {
sprintf (keyword,"AMDY%d",i+1);
wcs->y_coeff[i] = 0.0;
hgetr8 (hcoeff, keyword, &wcs->y_coeff[i]);
}
wcs->wcson = 1;
(void)strcpy (wcs->c1type, "RA");
(void)strcpy (wcs->c2type, "DEC");
(void)strcpy (wcs->ptype, "DSS");
wcs->degout = 0;
wcs->ndec = 3;
/* Compute a nominal reference pixel at the image center */
strcpy (wcs->ctype[0], "RA---DSS");
strcpy (wcs->ctype[1], "DEC--DSS");
wcs->crpix[0] = 0.5 * wcs->nxpix;
wcs->crpix[1] = 0.5 * wcs->nypix;
wcs->xrefpix = wcs->crpix[0];
wcs->yrefpix = wcs->crpix[1];
dsspos (wcs->crpix[0], wcs->crpix[1], wcs, &ra0, &dec0);
wcs->crval[0] = ra0;
wcs->crval[1] = dec0;
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
/* Compute a nominal scale factor */
dsspos (wcs->crpix[0], wcs->crpix[1]+1.0, wcs, &ra1, &dec1);
wcs->yinc = dec1 - dec0;
wcs->xinc = -wcs->yinc;
wcsioset (wcs);
/* Compute image rotation angle */
wcs->wcson = 1;
wcsrotset (wcs);
rot = degrad (wcs->rot);
/* Compute image scale at center */
dsspos (wcs->crpix[0]+cos(rot),
wcs->crpix[1]+sin(rot), wcs, &ra1, &dec1);
wcs->cdelt[0] = -wcsdist (ra0, dec0, ra1, dec1);
dsspos (wcs->crpix[0]+sin(rot),
wcs->crpix[1]+cos(rot), wcs, &ra1, &dec1);
wcs->cdelt[1] = wcsdist (ra0, dec0, ra1, dec1);
/* Set all other image scale parameters */
wcsdeltset (wcs, wcs->cdelt[0], wcs->cdelt[1], wcs->rot);
}
/* Approximate world coordinate system if plate scale is known */
else if ((ksearch (hstring,"SECPIX") != NULL ||
ksearch (hstring,"PIXSCALE") != NULL ||
ksearch (hstring,"PIXSCAL1") != NULL ||
ksearch (hstring,"XPIXSIZE") != NULL ||
ksearch (hstring,"SECPIX1") != NULL)) {
secpix = 0.0;
hgetr8 (hstring,"SECPIX",&secpix);
if (secpix == 0.0)
hgetr8 (hstring,"PIXSCALE",&secpix);
if (secpix == 0.0) {
hgetr8 (hstring,"SECPIX1",&secpix);
if (secpix != 0.0) {
cdelt1 = -secpix / 3600.0;
hgetr8 (hstring,"SECPIX2",&secpix);
cdelt2 = secpix / 3600.0;
}
else {
hgetr8 (hstring,"XPIXSIZE",&secpix);
if (secpix != 0.0) {
cdelt1 = -secpix / 3600.0;
hgetr8 (hstring,"YPIXSIZE",&secpix);
cdelt2 = secpix / 3600.0;
}
else {
hgetr8 (hstring,"PIXSCAL1",&secpix);
cdelt1 = -secpix / 3600.0;
hgetr8 (hstring,"PIXSCAL2",&secpix);
cdelt2 = secpix / 3600.0;
}
}
}
else {
cdelt2 = secpix / 3600.0;
cdelt1 = -cdelt2;
}
/* Get rotation angle from the header, if it's there */
rot = 0.0;
hgetr8 (hstring,"CROTA1", &rot);
if (wcs->rot == 0.)
hgetr8 (hstring,"CROTA2", &rot);
/* Set CD and PC matrices */
wcsdeltset (wcs, cdelt1, cdelt2, rot);
/* By default, set reference pixel to center of image */
wcs->crpix[0] = 0.5 + (wcs->nxpix * 0.5);
wcs->crpix[1] = 0.5 + (wcs->nypix * 0.5);
/* Get reference pixel from the header, if it's there */
if (ksearch (hstring,"CRPIX1") != NULL) {
hgetr8 (hstring,"CRPIX1",&wcs->crpix[0]);
hgetr8 (hstring,"CRPIX2",&wcs->crpix[1]);
}
/* Use center of detector array as reference pixel
else if (ksearch (hstring,"DETSIZE") != NULL ||
ksearch (hstring,"DETSEC") != NULL) {
char *ic;
hgets (hstring, "DETSIZE", 32, temp);
ic = strchr (temp, ':');
if (ic != NULL)
*ic = ' ';
ic = strchr (temp, ',');
if (ic != NULL)
*ic = ' ';
ic = strchr (temp, ':');
if (ic != NULL)
*ic = ' ';
ic = strchr (temp, ']');
if (ic != NULL)
*ic = (char) 0;
sscanf (temp, "%d %d %d %d", &idx1, &idx2, &idy1, &idy2);
dxrefpix = 0.5 * (double) (idx1 + idx2 - 1);
dyrefpix = 0.5 * (double) (idy1 + idy2 - 1);
hgets (hstring, "DETSEC", 32, temp);
ic = strchr (temp, ':');
if (ic != NULL)
*ic = ' ';
ic = strchr (temp, ',');
if (ic != NULL)
*ic = ' ';
ic = strchr (temp, ':');
if (ic != NULL)
*ic = ' ';
ic = strchr (temp, ']');
if (ic != NULL)
*ic = (char) 0;
sscanf (temp, "%d %d %d %d", &ix1, &ix2, &iy1, &iy2);
wcs->crpix[0] = dxrefpix - (double) (ix1 - 1);
wcs->crpix[1] = dyrefpix - (double) (iy1 - 1);
} */
wcs->xrefpix = wcs->crpix[0];
wcs->yrefpix = wcs->crpix[1];
wcs->crval[0] = -999.0;
if (!hgetra (hstring,"RA",&wcs->crval[0])) {
setwcserr ("WCSINIT: No RA with SECPIX, no WCS");
wcsfree (wcs);
return (NULL);
}
wcs->crval[1] = -999.0;
if (!hgetdec (hstring,"DEC",&wcs->crval[1])) {
setwcserr ("WCSINIT No DEC with SECPIX, no WCS");
wcsfree (wcs);
return (NULL);
}
wcs->xref = wcs->crval[0];
wcs->yref = wcs->crval[1];
wcs->coorflip = 0;
wcs->cel.ref[0] = wcs->crval[0];
wcs->cel.ref[1] = wcs->crval[1];
wcs->cel.ref[2] = 999.0;
if (!hgetr8 (hstring,"LONPOLE",&wcs->cel.ref[2]))
hgetr8 (hstring,"LONGPOLE",&wcs->cel.ref[2]);
wcs->cel.ref[3] = 999.0;
hgetr8 (hstring,"LATPOLE",&wcs->cel.ref[3]);
/* Coordinate reference frame and equinox */
(void) wcstype (wcs, "RA---TAN", "DEC--TAN");
wcs->coorflip = 0;
wcseq (hstring,wcs);
wcsioset (wcs);
wcs->degout = 0;
wcs->ndec = 3;
/* Epoch of image (from observation date, if possible) */
if (hgetr8 (hstring, "MJD-OBS", &mjd))
wcs->epoch = 1900.0 + (mjd - 15019.81352) / 365.242198781;
else if (!hgetdate (hstring,"DATE-OBS",&wcs->epoch)) {
if (!hgetdate (hstring,"DATE",&wcs->epoch)) {
if (!hgetr8 (hstring,"EPOCH",&wcs->epoch))
wcs->epoch = wcs->equinox;
}
}
wcs->wcson = 1;
}
else {
setwcserr ("WCSINIT: No image scale");
wcsfree (wcs);
return (NULL);
}
wcs->lin.crpix = wcs->crpix;
wcs->lin.cdelt = wcs->cdelt;
wcs->lin.pc = wcs->pc;
wcs->printsys = 1;
wcs->tabsys = 0;
wcs->linmode = 0;
/* Initialize special WCS commands */
setwcscom (wcs);
return (wcs);
}
/* Set coordinate system of image, input, and output */
static void
wcsioset (wcs)
struct WorldCoor *wcs;
{
if (strlen (wcs->radecsys) == 0 || wcs->prjcode == WCS_LIN)
strcpy (wcs->radecsys, "LINEAR");
wcs->syswcs = wcscsys (wcs->radecsys);
if (wcs->syswcs == WCS_B1950)
strcpy (wcs->radecout, "FK4");
else if (wcs->syswcs == WCS_J2000)
strcpy (wcs->radecout, "FK5");
else
strcpy (wcs->radecout, wcs->radecsys);
wcs->sysout = wcscsys (wcs->radecout);
wcs->eqout = wcs->equinox;
strcpy (wcs->radecin, wcs->radecsys);
wcs->sysin = wcscsys (wcs->radecin);
wcs->eqin = wcs->equinox;
return;
}
static void
wcseq (hstring, wcs)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
struct WorldCoor *wcs; /* World coordinate system data structure */
{
char mchar; /* Suffix character for one of multiple WCS */
mchar = (char) 0;
wcseqm (hstring, wcs, mchar);
return;
}
static void
wcseqm (hstring, wcs, mchar)
char *hstring; /* character string containing FITS header information
in the format <keyword>= <value> [/ <comment>] */
struct WorldCoor *wcs; /* World coordinate system data structure */
char mchar; /* Suffix character for one of multiple WCS */
{
int ieq = 0;
int eqhead = 0;
char systring[32], eqstring[32];
char radeckey[16], eqkey[16];
/* Set equinox from EQUINOX, EPOCH, or RADECSYS; default to 2000 */
systring[0] = 0;
eqstring[0] = 0;
if (mchar) {
sprintf (eqkey, "EQUINOX%c", mchar);
sprintf (radeckey, "RADECSYS%c", mchar);
}
else {
strcpy (eqkey, "EQUINOX");
sprintf (radeckey, "RADECSYS");
}
if (!hgets (hstring, eqkey, 31, eqstring)) {
if (hgets (hstring, "EQUINOX", 31, eqstring))
strcpy (eqkey, "EQUINOX");
}
if (!hgets (hstring, radeckey, 31, systring)) {
if (hgets (hstring, "RADECSYS", 31, systring))
sprintf (radeckey, "RADECSYS");
}
if (eqstring[0] == 'J') {
wcs->equinox = atof (eqstring+1);
ieq = atoi (eqstring+1);
strcpy (systring, "FK5");
}
else if (eqstring[0] == 'B') {
wcs->equinox = atof (eqstring+1);
ieq = (int) atof (eqstring+1);
strcpy (systring, "FK4");
}
else if (hgeti4 (hstring, eqkey, &ieq)) {
hgetr8 (hstring, eqkey, &wcs->equinox);
eqhead = 1;
}
else if (hgeti4 (hstring,"EPOCH",&ieq)) {
if (ieq == 0) {
ieq = 1950;
wcs->equinox = 1950.0;
}
else {
hgetr8 (hstring,"EPOCH",&wcs->equinox);
eqhead = 1;
}
}
else if (systring[0] != (char)0) {
if (!strncmp (systring,"FK4",3)) {
wcs->equinox = 1950.0;
ieq = 1950;
}
else if (!strncmp (systring,"ICRS",4)) {
wcs->equinox = 2000.0;
ieq = 2000;
}
else if (!strncmp (systring,"FK5",3)) {
wcs->equinox = 2000.0;
ieq = 2000;
}
else if (!strncmp (systring,"GAL",3)) {
wcs->equinox = 2000.0;
ieq = 2000;
}
else if (!strncmp (systring,"ECL",3)) {
wcs->equinox = 2000.0;
ieq = 2000;
}
}
if (ieq == 0) {
wcs->equinox = 2000.0;
ieq = 2000;
if (wcs->c1type[0] == 'R' || wcs->c1type[0] == 'D')
strcpy (systring,"FK5");
}
/* Epoch of image (from observation date, if possible) */
if (!hgetdate (hstring,"DATE-OBS",&wcs->epoch)) {
if (!hgetdate (hstring,"DATE",&wcs->epoch)) {
if (!hgetr8 (hstring,"EPOCH",&wcs->epoch))
wcs->epoch = wcs->equinox;
}
}
if (wcs->epoch == 0.0)
wcs->epoch = wcs->equinox;
/* Set coordinate system from keyword, if it is present */
if (systring[0] == (char) 0)
hgets (hstring, radeckey, 31, systring);
if (systring[0] != (char) 0) {
strcpy (wcs->radecsys,systring);
if (!eqhead) {
if (!strncmp (wcs->radecsys,"FK4",3))
wcs->equinox = 1950.0;
else if (!strncmp (wcs->radecsys,"FK5",3))
wcs->equinox = 2000.0;
else if (!strncmp (wcs->radecsys,"ICRS",4))
wcs->equinox = 2000.0;
else if (!strncmp (wcs->radecsys,"GAL",3) && ieq == 0)
wcs->equinox = 2000.0;
}
}
/* Otherwise set coordinate system from equinox */
/* Systemless coordinates cannot be translated using b, j, or g commands */
else if (wcs->syswcs != WCS_NPOLE) {
if (ieq > 1980)
strcpy (wcs->radecsys,"FK5");
else
strcpy (wcs->radecsys,"FK4");
}
/* Set galactic coordinates if GLON or GLAT are in C1TYPE */
if (wcs->c1type[0] == 'G')
strcpy (wcs->radecsys,"GALACTIC");
else if (wcs->c1type[0] == 'E')
strcpy (wcs->radecsys,"ECLIPTIC");
else if (wcs->c1type[0] == 'S')
strcpy (wcs->radecsys,"SGALACTC");
else if (wcs->c1type[0] == 'H')
strcpy (wcs->radecsys,"HELIOECL");
else if (wcs->c1type[0] == 'A')
strcpy (wcs->radecsys,"ALTAZ");
else if (wcs->c1type[0] == 'L')
strcpy (wcs->radecsys,"LINEAR");
wcs->syswcs = wcscsys (wcs->radecsys);
return;
}
/* Jun 11 1998 Split off header-dependent WCS initialization from other subs
* Jun 15 1998 Fix major bug in wcsinit() when synthesizing WCS from header
* Jun 18 1998 Fix bug in CD initialization; split PC initialization off
* Jun 18 1998 Split PC initialization off into subroutine wcspcset()
* Jun 24 1998 Set equinox from RADECSYS only if EQUINOX and EPOCH not present
* Jul 6 1998 Read third and fourth axis CTYPEs
* Jul 7 1998 Initialize eqin and eqout to equinox,
* Jul 9 1998 Initialize rotation matrices correctly
* Jul 13 1998 Initialize rotation, scale for polynomial and DSS projections
* Aug 6 1998 Fix CROTA computation for DSS projection
* Sep 4 1998 Fix CROTA, CDELT computation for DSS and polynomial projections
* Sep 14 1998 If DATE-OBS not found, check for DATE
* Sep 14 1998 If B or J present in EQUINOX, use that info to set system
* Sep 29 1998 Initialize additional WCS commands from the environment
* Sep 29 1998 Fix bug which read DATE as number rather than formatted date
* Dec 2 1998 Read projection constants from header (bug fix)
*
* Feb 9 1999 Set rotation angle correctly when using DSS projection
* Feb 19 1999 Fill in CDELTs from scale keyword if absent or zero
* Feb 19 1999 Add PIXSCALE as possible default arcseconds per pixel
* Apr 7 1999 Add error checking for NAXIS and NAXIS1 keywords
* Apr 7 1999 Do not set systring if epoch is 0 and not RA/Dec
* Jul 8 1999 In RADECSYS, use FK5 and FK4 instead of J2000 and B1950
* Oct 15 1999 Free wcs using wcsfree()
* Oct 20 1999 Add multiple WCS support using new subroutine names
* Oct 21 1999 Delete unused variables after lint; declare dsspos()
* Nov 9 1999 Add wcschar() to check WCSNAME keywords for desired WCS
* Nov 9 1999 Check WCSPREx keyword to find out if chained WCS's
*
* Jan 6 1999 Add wcsinitn() to initialize from specific WCSNAME
* Jan 24 2000 Set CD matrix from header even if using polynomial
* Jan 27 2000 Fix MJD to epoch conversion for when MJD-OBS is the only date
* Jan 28 2000 Set CD matrix for DSS projection, too
* Jan 28 2000 Use wcsproj instead of oldwcs
* Dec 18 2000 Fix error in hgets() call in wcschar()
* Dec 29 2000 Compute inverse CD matrix even if polynomial solution
* Dec 29 2000 Add PROJR0 keyword for WCSLIB projections
* Dec 29 2000 Use CDi_j matrix if any elements are present
*
* Jan 31 2001 Fix to allow 1D WCS
* Jan 31 2001 Treat single character WCS name as WCS character
* Feb 20 2001 Implement WCSDEPx nested WCS's
* Feb 23 2001 Initialize all 4 terms of CD matrix
* Feb 28 2001 Fix bug which read CRPIX1 into CRPIX2
* Mar 20 2001 Compare mchar to (char)0, not null
* Mar 21 2001 Move ic declaration into commented out code
* Jul 12 2001 Read PROJPn constants into proj.p array instead of PVn
* Sep 7 2001 Set system to galactic or ecliptic based on CTYPE, not RADECSYS
* Oct 11 2001 Set ctype[0] as well as ctype[1] to TAN for TNX projections
* Oct 19 2001 WCSDIM keyword overrides zero value of NAXIS
*
* Feb 19 2002 Add XPIXSIZE/YPIXSIZE (KPNO) as default image scale keywords
* Mar 12 2002 Add LONPOLE as well as LONGPOLE for WCSLIB 2.8
* Apr 3 2002 Implement hget8c() and hgetsc() to simplify code
* Apr 3 2002 Add PVj_n projection constants in addition to PROJPn
* Apr 19 2002 Increase numeric keyword value length from 16 to 31
* Apr 19 2002 Fix bug which didn't set radecsys keyword name
* Apr 24 2002 If no WCS present for specified letter, return null
* Apr 26 2002 Implement WCSAXESa keyword as first choice for number of axes
* Apr 26 2002 Add wcschar and wcsname to WCS structure
* May 9 2002 Add radvel and zvel to WCS structure
* May 13 2002 Free everything which is allocated
* May 28 2002 Read 10 prj.p instead of maximum of 100
* May 31 2002 Fix bugs with PV reading
* May 31 2002 Initialize syswcs, sysin, sysout in wcsioset()
* Sep 25 2002 Fix subroutine calls for radvel and latpole
* Dec 6 2002 Correctly compute pixel at center of image for default CRPIX
*
* Jan 2 2002 Do not reinitialize projection vector for PV input
* Jan 3 2002 For ZPN, read PVi_0 to PVi_9, not PVi_1 to PVi_10
* Mar 27 2003 Clean up default center computation
* Apr 3 2003 Add input for SIRTF distortion coefficients
* May 8 2003 Change PROJP reading to start with 0 instead of 1
* May 22 2003 Add ZPX approximation, reading projpn from WATi
* May 28 2003 Avoid reinitializing coefficients set by PROJP
* Jun 26 2003 Initialize xref and yref to -999.0
* Sep 23 2003 Change mgets() to mgetstr() to avoid name collision at UCO Lick
* Oct 1 2003 Rename wcs->naxes to wcs->naxis to match WCSLIB 3.2
* Nov 3 2003 Initialize distortion coefficients in distortinit() in distort.c
* Dec 1 2003 Change p[0,1,2] initializations to p[1,2,3]
*/
syntax highlighted by Code2HTML, v. 0.9.1