/*** 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 /* strstr, NULL */ #include /* stderr */ #include #include "wcs.h" #ifndef VMS #include #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 = [/ ] */ 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 = [/ ] */ 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 = [/ ] */ 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 = [/ ] */ 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 = [/ ] */ 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 = [/ ] */ { 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 = [/ ] */ 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 = [/ ] */ 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 = [/ ] */ 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 */