Merge commit '339420dd5dd874c41f6bab5808291fb4036dd022' as 'funtools'

1 files changed, 1611 insertions, 0 deletions

diff --git a/funtools/wcs/wcsinit.c b/funtools/wcs/wcsinit.c
new file mode 100644
index 0000000..134c1bf
--- /dev/null
+++ b/funtools/wcs/wcsinit.c
@@ -0,0 +1,1611 @@
+/*** File libwcs/wcsinit.c
+ *** October 19, 2012
+ *** By Jessica Mink, jmink@cfa.harvard.edu
+ *** Harvard-Smithsonian Center for Astrophysics
+ *** Copyright (C) 1998-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:	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 wcschar();
+
+/* set up a WCS structure from a FITS image header lhstring bytes long 
+ * for a specified WCS name */
+
+struct WorldCoor *
+wcsninitn (hstring, lhstring, name)
+
+const char *hstring;	/* character string containing FITS header information
+		   	in the format <keyword>= <value> [/ <comment>] */
+int	lhstring;	/* Length of FITS header in bytes */
+const 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)
+
+const char *hstring;	/* character string containing FITS header information
+			   in the format <keyword>= <value> [/ <comment>] */
+const 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)
+
+const char *hstring;	/* character string containing FITS header information
+		   	in the format <keyword>= <value> [/ <comment>] */
+const 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)
+
+const 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)
+
+const 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[0] == ' ')
+	mchar[0] = (char) 0;
+    return (wcsinitc (hstring, mchar));
+}
+
+
+/* set up a WCS structure from a FITS image header */
+
+struct WorldCoor *
+wcsinit (hstring)
+
+const 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, wchar)
+
+const char *hstring;	/* character string containing FITS header information
+			   in the format <keyword>= <value> [/ <comment>] */
+char *wchar;		/* Suffix character for one of multiple WCS */
+{
+    struct WorldCoor *wcs, *depwcs;
+    char ctype1[32], ctype2[32], tstring[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[81];
+    char keyword[16];
+    int ieq, i, j, k, 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[80];
+    char wcsname[64];	/* Name of WCS depended on by current WCS */
+    char mchar;
+    char cspace = (char) ' ';
+    char cnull = (char) 0;
+    double mjd;
+    double rot;
+    double ut;
+    int nax;
+    int twod;
+    extern int tnxinit();
+    extern int zpxinit();
+    extern int platepos();
+    extern int dsspos();
+    void invert_wcs();
+
+    wcs = (struct WorldCoor *) calloc (1, sizeof(struct WorldCoor));
+
+    /* Set WCS character and name in structure */
+    mchar = wchar[0];
+    if (mchar == ' ')
+	mchar = cnull;
+    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;
+    wcs->wcsl.cubeface = -1;
+
+    /* 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);
+	}
+    if (naxes > 2)
+	naxes = 2;
+    wcs->naxis = naxes;
+    wcs->naxes = naxes;
+    wcs->lin.naxis = naxes;
+    wcs->nxpix = 0;
+    hgetr8 (hstring, "NAXIS1", &wcs->nxpix);
+    if (wcs->nxpix < 1)
+	hgetr8 (hstring, "IMAGEW", &wcs->nxpix);
+    if (wcs->nxpix < 1) {
+	setwcserr ("WCSINIT: No NAXIS1 or IMAGEW keyword");
+	wcsfree (wcs);
+	return (NULL);
+	}
+    wcs->nypix = 0;
+    hgetr8 (hstring, "NAXIS2", &wcs->nypix);
+    if (wcs->nypix < 1)
+	hgetr8 (hstring, "IMAGEH", &wcs->nypix);
+    if (naxes > 1 && wcs->nypix < 1) {
+	setwcserr ("WCSINIT: No NAXIS2 or IMAGEH keyword");
+	wcsfree (wcs);
+	return (NULL);
+	}
+
+    /* Reset number of axes to only those with dimension greater than one */
+    nax = 0;
+    for (i = 0; i < naxes; i++) {
+
+	/* Check for number of pixels in axis more than one */
+	strcpy (keyword, "NAXIS");
+	sprintf (temp, "%d", i+1);
+	strcat (keyword, temp);
+	if (!hgeti4 (hstring, keyword, &j)) {
+	    if (i == 0 && wcs->nxpix > 1) {
+		/* fprintf (stderr,"WCSINIT: Missing keyword %s set to %.0f from IMAGEW\n",
+			 keyword, wcs->nxpix); */
+		j = wcs->nxpix;
+		}
+	    else if (i == 1 && wcs->nypix > 1) {
+		/* fprintf (stderr,"WCSINIT: Missing keyword %s set to %.0f from IMAGEH\n",
+			 keyword, wcs->nypix); */
+		j = wcs->nypix;
+		}
+	    else
+		fprintf (stderr,"WCSINIT: Missing keyword %s assumed 1\n",keyword);
+	    }
+
+	/* Check for TAB WCS in axis */
+	strcpy (keyword, "CTYPE");
+	strcat (keyword, temp);
+	if (hgets (hstring, keyword, 16, temp)) {
+	    if (strsrch (temp, "-TAB"))
+		j = 0;
+	    }
+	if (j > 1) nax = nax + 1;
+	}
+    naxes = nax;
+    wcs->naxes = nax;
+    wcs->naxis = nax;
+
+    hgets (hstring, "INSTRUME", 16, wcs->instrument);
+    hgeti4 (hstring, "DETECTOR", &wcs->detector);
+    wcs->wcsproj = getdefwcs();
+    wcs->logwcs = 0;
+    hgeti4 (hstring, "DC-FLAG", &wcs->logwcs);
+
+    /* Initialize rotation matrices */
+    for (i = 0; i < 81; i++) wcs->pc[i] = 0.0;
+    for (i = 0; i < 81; i++) pc[i] = 0.0;
+    for (i = 0; i < naxes; i++) wcs->pc[(i*naxes)+i] = 1.0;
+    for (i = 0; i < naxes; i++) pc[(i*naxes)+i] = 1.0;
+    for (i = 0; i < 9; i++) wcs->cdelt[i] = 0.0;
+    for (i = 0; i < naxes; i++) wcs->cdelt[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;
+
+    for (i = 0; i < 10; i++) {
+	wcs->prj.p[i] = 0.0;
+	}
+
+    /* World coordinate system reference coordinate information */
+    if (hgetsc (hstring, "CTYPE1", &mchar, 16, ctype1)) {
+
+	/* Read second coordinate type */
+	strcpy (ctype2, ctype1);
+	if (!hgetsc (hstring, "CTYPE2", &mchar, 16, ctype2))
+	    twod = 0;
+	else
+	    twod = 1;
+	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 (!strcmp (wcs->units[0], "pixel"))
+		wcs->prjcode = WCS_PIX;
+	    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;
+			}
+		    }
+		if (!strcmp (wcs->units[0], "pixel"))
+		    wcs->prjcode = WCS_PIX;
+		}
+	    }
+
+	/* 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;
+	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++) {
+	    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]);
+		}
+	    }
+
+	/* 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;
+		}
+	    }
+
+	/* Initialize ZPX, defaulting to ZPN if there is a problem */
+	if (wcs->prjcode == WCS_ZPX) {
+	    if (zpxinit (hstring, wcs)) {
+		wcs->ctype[0][7] = 'N';
+		wcs->ctype[1][7] = 'N';
+		wcs->prjcode = WCS_ZPN;
+		}
+	    }
+
+	/* Set TPV to TAN as SCAMP coefficients will be added below */
+	if (wcs->prjcode == WCS_TPV) {
+	    wcs->ctype[0][6] = 'A';
+	    wcs->ctype[0][7] = 'N';
+	    wcs->ctype[1][6] = 'A';
+	    wcs->ctype[1][7] = 'N';
+	    wcs->prjcode = WCS_TAN;
+	    }
+
+	/* Coordinate reference frame, equinox, and epoch */
+	if (wcs->wcsproj > 0)
+	    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;
+
+	    /* Initialize rotation matrix */
+	    for (i = 0; i < 81; i++) {
+		pc[i] = 0.0;
+		wcs->pc[i] = 0.0;
+		}
+	    for (i = 0; i < naxes; i++)
+		pc[(i*naxes)+i] = 1.0;
+
+	    /* Read FITS WCS interim rotation matrix */
+	    if (!mchar && hgetr8 (hstring,"PC001001",&pc[0]) != 0) {
+		k = 0;
+		for (i = 0; i < naxes; i++) {
+		    for (j = 0; j < naxes; j++) {
+			if (i == j)
+			    pc[k] = 1.0;
+			else
+			    pc[k] = 0.0;
+			sprintf (keyword, "PC00%1d00%1d", i+1, j+1);
+			hgetr8 (hstring, keyword, &pc[k++]);
+			}
+		    }
+		wcspcset (wcs, cdelt1, cdelt2, pc);
+		}
+
+	    /* Read FITS WCS standard rotation matrix */
+	    else if (hgetr8c (hstring, "PC1_1", &mchar, &pc[0]) != 0) {
+		k = 0;
+		for (i = 0; i < naxes; i++) {
+		    for (j = 0; j < naxes; j++) {
+			if (i == j)
+			    pc[k] = 1.0;
+			else
+			    pc[k] = 0.0;
+			sprintf (keyword, "PC%1d_%1d", i+1, j+1);
+			hgetr8c (hstring, keyword, &mchar, &pc[k++]);
+			}
+		    }
+		wcspcset (wcs, cdelt1, cdelt2, pc);
+		}
+
+	    /* Otherwise, use CROTAn */
+	    else {
+		rot = 0.0;
+		if (ilat == 2)
+		    hgetr8c (hstring, "CROTA2", &mchar, &rot);
+		else
+		    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");
+	    }
+
+	/* SCAMP convention */
+	if (wcs->prjcode == WCS_TAN && wcs->naxis == 2) { 
+	    int n = 0;
+	    if (wcs->inv_x) {
+		poly_end(wcs->inv_x);
+		wcs->inv_x = NULL;
+		}
+	    if (wcs->inv_y) {
+		poly_end(wcs->inv_y);
+		wcs->inv_y = NULL;
+		}
+	    wcs->pvfail = 0;
+	    for (i = 0; i < (2*MAXPV); i++) {
+		wcs->projppv[i] = 0.0;
+		wcs->prj.ppv[i] = 0.0;
+		}
+	    for (k = 0; k < 2; k++) {
+		for (j = 0; j < MAXPV; j++) {
+		    sprintf(keyword, "PV%d_%d", k+1, j);
+		    if (hgetr8c(hstring, keyword,&mchar, &wcs->projppv[j+k*MAXPV]) == 0) {
+			wcs->projppv[j+k*MAXPV] = 0.0;
+			}
+		    else
+			n++;
+		    }
+		}
+
+	    /* If any PVi_j are set, add them to the structure */
+	    if (n > 0) {
+		n = 0;
+
+		for (k = MAXPV; k >= 0; k--) {
+		    /* lat comes first for compatibility reasons */
+		    wcs->prj.ppv[k] = wcs->projppv[k+wcs->wcsl.lat*MAXPV];
+		    wcs->prj.ppv[k+MAXPV] = wcs->projppv[k+wcs->wcsl.lng*MAXPV];
+		    if (!n && (wcs->prj.ppv[k] || wcs->prj.ppv[k+MAXPV]))  {
+			n = k+1;
+			}
+		    }
+		invert_wcs(wcs);
+
+		/* Need to call tanset again */
+		wcs->cel.flag = 0;
+		}
+	    }
+
+	/* If linear or pixel WCS, print "degrees" */
+	if (!strncmp (wcs->ptype,"LINEAR",6) ||
+	    !strncmp (wcs->ptype,"PIXEL",5)) {
+	    wcs->degout = -1;
+	    wcs->ndec = 5;
+	    }
+
+	/* 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;
+		}
+	    }
+
+	/* Add time of day if not part of DATE-OBS string */
+	else {
+	    hgets (hstring,"DATE-OBS",32,tstring);
+	    if (!strchr (tstring,'T')) {
+		if (hgetr8 (hstring, "UT",&ut))
+		    wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
+		else if (hgetr8 (hstring, "UTMID",&ut))
+		    wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
+		}
+	    }
+
+	wcs->wcson = 1;
+	}
+
+    else if (mchar != cnull && mchar != cspace) {
+	(void) sprintf (temp, "WCSINITC: No image scale for WCS %c", mchar);
+	setwcserr (temp);
+	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 = cnull;
+	    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 = cnull;
+	    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]);
+
+	/* 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;
+		}
+	    }
+
+	/* Add time of day if not part of DATE-OBS string */
+	else {
+	    hgets (hstring,"DATE-OBS",32,tstring);
+	    if (!strchr (tstring,'T')) {
+		if (hgetr8 (hstring, "UT",&ut))
+		    wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
+		else if (hgetr8 (hstring, "UTMID",&ut))
+		    wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
+		}
+	    }
+
+	/* 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;
+	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);
+}
+
+
+/******* invert_wcs ***********************************************************
+PROTO	void invert_wcs(wcsstruct *wcs)
+PURPOSE	Invert WCS projection mapping (using a polynomial).
+INPUT	WCS structure.
+OUTPUT	-.
+NOTES	.
+AUTHOR	E. Bertin (IAP)
+VERSION	06/11/2003
+ ***/
+
+void
+invert_wcs( struct WorldCoor *wcs)
+
+{
+    polystruct	*poly;
+    double pixin[NAXISPV],raw[NAXISPV],rawmin[NAXISPV];
+    double *outpos,*outpost, *lngpos,*lngpost;
+    double *latpos,*latpost,lngstep,latstep, rawsize, epsilon;
+    int	 group[] = {1,1};
+    /* Don't ask, this is needed by poly_init()! */
+    int	 i,j,lng,lat,deg,  maxflag;
+    char errstr[80];
+    double xmin;
+    double ymin;
+    double xmax;
+    double ymax;
+    double lngmin;
+    double latmin;
+
+    /* Check first that inversion is not straightforward */
+    lng = wcs->wcsl.lng;
+    lat = wcs->wcsl.lat;
+
+    if (wcs->naxis != NAXISPV) {
+	return;
+	}
+
+    if (strcmp(wcs->wcsl.pcode, "TAN") != 0) {
+	return;
+	}
+
+    if ((wcs->projppv[1+lng*MAXPV] == 0) &&
+	(wcs->projppv[1+lat*MAXPV] == 0)) {
+	return;
+	}
+
+    if (wcs->wcs != NULL) {
+	pix2wcs(wcs->wcs,0,0,&xmin,&ymin);
+	pix2wcs(wcs->wcs,wcs->nxpix,wcs->nypix,&xmax,&ymax);
+	}
+    else {
+	xmin = 0;
+	ymin = 0;
+	xmax = wcs->nxpix;
+	ymax = wcs->nypix;
+	}
+
+    /* We define x as "longitude" and y as "latitude" projections */
+    /* We assume that PCxx cross-terms with additional dimensions are small */
+    /* Sample the whole image with a regular grid */
+    if (lng == 0) {
+	lngstep = (xmax-xmin)/(WCS_NGRIDPOINTS-1.0);
+	lngmin  = xmin;
+	latstep = (ymax-ymin)/(WCS_NGRIDPOINTS-1.0);
+	latmin  = ymin;
+	}
+    else {
+	lngstep = (ymax-ymin)/(WCS_NGRIDPOINTS-1.0);
+	lngmin = ymin;
+	latstep = (xmax-xmin)/(WCS_NGRIDPOINTS-1.0);
+	latmin = xmin;
+	}
+
+    outpos = (double *)calloc(2*WCS_NGRIDPOINTS2,sizeof(double));
+    lngpos = (double *)calloc(WCS_NGRIDPOINTS2,sizeof(double));
+    latpos = (double *)calloc(WCS_NGRIDPOINTS2,sizeof(double));
+    raw[lat] = rawmin[lat] = 0.5+latmin;
+    raw[lng] = rawmin[lng] = 0.5+lngmin;
+    outpost = outpos;
+    lngpost = lngpos;
+    latpost = latpos;
+    for (j=WCS_NGRIDPOINTS; j--; raw[lat]+=latstep) {
+  	raw[lng] = rawmin[lng];
+	for (i=WCS_NGRIDPOINTS; i--; raw[lng]+=lngstep) {
+	    if (linrev(raw, &wcs->lin, pixin)) {
+		sprintf (errstr,"*Error*: incorrect linear conversion in %s",
+			 wcs->wcsl.pcode);
+		setwcserr (errstr);
+		}
+	    *(lngpost++) = pixin[lng];
+	    *(latpost++) = pixin[lat];
+	    raw_to_pv (&wcs->prj,pixin[lng],pixin[lat], outpost, outpost+1);
+	    outpost += 2;
+	    }
+	}
+
+    /* Invert "longitude" */
+    /* Compute the extent of the pixel in reduced projected coordinates */
+    linrev(rawmin, &wcs->lin, pixin);
+    pixin[lng] += S2D;
+    linfwd(pixin, &wcs->lin, raw);
+    rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
+                 +(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*D2S;
+    if (!rawsize) {
+	sprintf (errstr,"*Error*: incorrect linear conversion in %s",
+		 wcs->wcsl.pcode);
+	setwcserr (errstr);
+	}
+    epsilon = WCS_INVACCURACY/rawsize;
+
+    /* Find the lowest degree polynom */
+    poly = NULL;  /* to avoid gcc -Wall warnings */
+    maxflag = 1;
+    for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++) {
+	if (deg>1) {
+	    poly_end(poly);
+	    }
+	poly = poly_init(group, 2, &deg, 1);
+	poly_fit(poly, outpos, lngpos, NULL, WCS_NGRIDPOINTS2, NULL);
+	maxflag = 0;
+	outpost = outpos;
+	lngpost = lngpos;
+  	for (i=WCS_NGRIDPOINTS2; i--; outpost+=2) {
+	    if (fabs(poly_func(poly, outpost)-*(lngpost++))>epsilon) {
+		maxflag = 1;
+		break;
+		}
+	    }
+	}
+    if (maxflag) {
+	setwcserr ("WARNING: Significant inaccuracy likely to occur in projection");
+	wcs->pvfail = 1;
+	}
+
+    /* Now link the created structure */
+    wcs->prj.inv_x = wcs->inv_x = poly;
+
+    /* Invert "latitude" */
+    /* Compute the extent of the pixel in reduced projected coordinates */
+    linrev(rawmin, &wcs->lin, pixin);
+    pixin[lat] += S2D;
+    linfwd(pixin, &wcs->lin, raw);
+    rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
+                 +(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*D2S;
+    if (!rawsize) {
+	sprintf (errstr,"*Error*: incorrect linear conversion in %s",
+		 wcs->wcsl.pcode);
+	setwcserr (errstr);
+	}
+    epsilon = WCS_INVACCURACY/rawsize;
+
+    /* Find the lowest degree polynom */
+    maxflag = 1;
+    for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++) {
+	if (deg>1)
+	    poly_end(poly);
+	poly = poly_init(group, 2, &deg, 1);
+	poly_fit(poly, outpos, latpos, NULL, WCS_NGRIDPOINTS2, NULL);
+	maxflag = 0;
+	outpost = outpos;
+	latpost = latpos;
+	for (i=WCS_NGRIDPOINTS2; i--; outpost+=2) {
+	    if (fabs(poly_func(poly, outpost)-*(latpost++))>epsilon) {
+		maxflag = 1;
+		break;
+		}
+	    }
+	}
+    if (maxflag) {
+	setwcserr ("WARNING: Significant inaccuracy likely to occur in projection");
+	wcs->pvfail = 1;
+	}
+
+    /* Now link the created structure */
+    wcs->prj.inv_y = wcs->inv_y = poly;
+
+    /* Free memory */
+    free(outpos);
+    free(lngpos);
+    free(latpos);
+
+    return;
+}
+
+
+/* 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");
+    if (wcs->prjcode == WCS_PIX)
+	strcpy (wcs->radecsys, "PIXEL");
+    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];
+    char tstring[32];
+    double ut;
+
+    /* Set equinox from EQUINOX, EPOCH, or RADECSYS; default to 2000 */
+    systring[0] = 0;
+    eqstring[0] = 0;
+    if (mchar[0]) {
+	sprintf (eqkey, "EQUINOX%c", mchar[0]);
+	sprintf (radeckey, "RADECSYS%c", mchar[0]);
+	}
+    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 (!strncmp (wcs->c1type, "RA",2) || !strncmp (wcs->c1type,"DEC",3))
+	    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;
+	    }
+	}
+
+	/* Add time of day if not part of DATE-OBS string */
+    else {
+	hgets (hstring,"DATE-OBS",32,tstring);
+	if (!strchr (tstring,'T')) {
+	    if (hgetr8 (hstring, "UT",&ut))
+		wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
+	    else if (hgetr8 (hstring, "UTMID",&ut))
+		wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
+	    }
+	}
+    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]
+ * Dec  3 2003	Add back wcs->naxes for backward compatibility
+ * Dec  3 2003	Remove unused variables j,m in wcsinitc()
+ * Dec 12 2003	Fix call to setwcserr() with format in it
+ *
+ * Feb 26 2004	Add parameters for ZPX projection
+ *
+ * Jun 22 2005	Drop declaration of variable wcserrmsg which is not used
+ * Nov  9 2005	Use CROTA1 if CTYPE1 is LAT/DEC, CROTA2 if CTYPE2 is LAT/DEC
+ *
+ * Mar  9 2006	Get Epoch of observation from MJD-OBS or DATE-OBS/UT unless DSS
+ * Apr 24 2006	Initialize rotation matrices
+ * Apr 25 2006	Ignore axes with dimension of one
+ * May 19 2006	Initialize all of 9x9 PC matrix; read in loops
+ * Aug 21 2006	Limit naxes to 2 everywhere; RA and DEC should always be 1st
+ * Oct  6 2006	If units are pixels, projection type is PIXEL
+ * Oct 30 2006	Initialize cube face to -1, not a cube projection
+ *
+ * Jan  4 2007	Drop declarations of wcsinitc() and wcsinitn() already in wcs.h
+ * Jan  8 2007	Change WCS letter from char to char*
+ * Feb  1 2007	Read IRAF log wavelength flag DC-FLAG to wcs.logwcs
+ * Feb 15 2007	Check for wcs->wcsproj > 0 instead of CTYPEi != LINEAR or PIXEL
+ * Mar 13 2007	Try for RA, DEC, SECPIX if WCS character is space or null
+ * Apr 27 2007	Ignore axes with TAB WCS for now
+ * Oct 17 2007	Fix bug testing &mchar instead of mchar in if statement
+ *
+ * May  9 2008	Initialize TNX projection when projection types first set
+ * Jun 27 2008	If NAXIS1 and NAXIS2 not present, check for IMAGEW and IMAGEH
+ *
+ * Mar 24 2009	Fix dimension bug if NAXISi not present (fix from John Burns)
+ *
+ * Mar 11 2011	Add NOAO ZPX projection (Frank Valdes)
+ * Mar 18 2011	Add invert_wcs() by Emmanuel Bertin for SCAMP
+ * Mar 18 2011	Change Bertin's ARCSEC/DEG to S2D and DEG/ARCSEC to D2S
+ * Sep  1 2011	Add TPV as TAN with SCAMP PVs
+ *
+ * Oct 19 2012	Drop unused variable iszpx; fix bug in latmin assignment
+ */