diff options
Diffstat (limited to 'funtools/wcs/wcsinit.c')
| -rw-r--r-- | funtools/wcs/wcsinit.c | 1611 |
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, °, 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, °, 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 + */ |