Repatch (from linux) of OSX IRAF

1 files changed, 786 insertions, 0 deletions

diff --git a/pkg/tbtables/cfitsio/wcsutil.c.OLD b/pkg/tbtables/cfitsio/wcsutil.c.OLD
new file mode 100644
index 00000000..c5b87fa8
--- /dev/null
+++ b/pkg/tbtables/cfitsio/wcsutil.c.OLD
@@ -0,0 +1,786 @@
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "fitsio2.h"
+/*--------------------------------------------------------------------------*/
+int ffgics(fitsfile *fptr,    /* I - FITS file pointer           */
+           double *xrval,     /* O - X reference value           */
+           double *yrval,     /* O - Y reference value           */
+           double *xrpix,     /* O - X reference pixel           */
+           double *yrpix,     /* O - Y reference pixel           */
+           double *xinc,      /* O - X increment per pixel       */
+           double *yinc,      /* O - Y increment per pixel       */
+           double *rot,       /* O - rotation angle (degrees)    */
+           char *type,        /* O - type of projection ('-tan') */
+           int *status)       /* IO - error status               */
+/*
+       read the values of the celestial coordinate system keywords.
+       These values may be used as input to the subroutines that
+       calculate celestial coordinates. (ffxypx, ffwldp)
+
+       Modified in Nov 1999 to convert the CD matrix keywords back
+       to the old CDELTn form, and to swap the axes if the dec-like
+       axis is given first, and to assume default values if any of the
+       keywords are not present.
+*/
+{
+    int tstat = 0, cd_exists = 0, pc_exists = 0;
+    char ctype[FLEN_VALUE];
+    double cd11 = 0.0, cd21 = 0.0, cd22 = 0.0, cd12 = 0.0;
+    double pc11 = 1.0, pc21 = 0.0, pc22 = 1.0, pc12 = 0.0;
+    double pi =  3.1415926535897932;
+    double phia, phib, temp;
+    double toler = .0002;  /* tolerance for angles to agree (radians) */
+                           /*   (= approximately 0.01 degrees) */
+
+    if (*status > 0)
+       return(*status);
+
+    tstat = 0;
+    if (ffgkyd(fptr, "CRVAL1", xrval, NULL, &tstat))
+       *xrval = 0.;
+
+    tstat = 0;
+    if (ffgkyd(fptr, "CRVAL2", yrval, NULL, &tstat))
+       *yrval = 0.;
+
+    tstat = 0;
+    if (ffgkyd(fptr, "CRPIX1", xrpix, NULL, &tstat))
+        *xrpix = 0.;
+
+    tstat = 0;
+    if (ffgkyd(fptr, "CRPIX2", yrpix, NULL, &tstat))
+        *yrpix = 0.;
+
+    /* look for CDELTn first, then CDi_j keywords */
+    tstat = 0;
+    if (ffgkyd(fptr, "CDELT1", xinc, NULL, &tstat))
+    {
+        /* CASE 1: no CDELTn keyword, so look for the CD matrix */
+        tstat = 0;
+        if (ffgkyd(fptr, "CD1_1", &cd11, NULL, &tstat))
+            tstat = 0;  /* reset keyword not found error */
+        else
+            cd_exists = 1;  /* found at least 1 CD_ keyword */
+
+        if (ffgkyd(fptr, "CD2_1", &cd21, NULL, &tstat))
+            tstat = 0;  /* reset keyword not found error */
+        else
+            cd_exists = 1;  /* found at least 1 CD_ keyword */
+
+        if (ffgkyd(fptr, "CD1_2", &cd12, NULL, &tstat))
+            tstat = 0;  /* reset keyword not found error */
+        else
+            cd_exists = 1;  /* found at least 1 CD_ keyword */
+
+        if (ffgkyd(fptr, "CD2_2", &cd22, NULL, &tstat))
+            tstat = 0;  /* reset keyword not found error */
+        else
+            cd_exists = 1;  /* found at least 1 CD_ keyword */
+
+        if (cd_exists)  /* convert CDi_j back to CDELTn */
+        {
+            /* there are 2 ways to compute the angle: */
+            phia = atan2( cd21, cd11);
+            phib = atan2(-cd12, cd22);
+
+            /* ensure that phia <= phib */
+            temp = minvalue(phia, phib);
+            phib = maxvalue(phia, phib);
+            phia = temp;
+
+            /* there is a possible 180 degree ambiguity in the angles */
+            /* so add 180 degress to the smaller value if the values  */
+            /* differ by more than 90 degrees = pi/2 radians.         */
+            /* (Later, we may decide to take the other solution by    */
+            /* subtracting 180 degrees from the larger value).        */
+
+            if ((phib - phia) > (pi / 2.))
+               phia += pi;
+
+            if (fabs(phia - phib) > toler) 
+            {
+               /* angles don't agree, so looks like there is some skewness */
+               /* between the axes.  Return with an error to be safe. */
+               *status = APPROX_WCS_KEY;
+            }
+      
+            phia = (phia + phib) /2.;  /* use the average of the 2 values */
+            *xinc = cd11 / cos(phia);
+            *yinc = cd22 / cos(phia);
+            *rot = phia * 180. / pi;
+
+            /* common usage is to have a positive yinc value.  If it is */
+            /* negative, then subtract 180 degrees from rot and negate  */
+            /* both xinc and yinc.  */
+
+            if (*yinc < 0)
+            {
+                *xinc = -(*xinc);
+                *yinc = -(*yinc);
+                *rot = *rot - 180.;
+            }
+        }
+        else   /* no CD matrix keywords either */
+        {
+            *xinc = 1.;
+
+            /* there was no CDELT1 keyword, but check for CDELT2 just in case */
+            tstat = 0;
+            if (ffgkyd(fptr, "CDELT2", yinc, NULL, &tstat))
+                *yinc = 1.;
+
+            tstat = 0;
+            if (ffgkyd(fptr, "CROTA2", rot, NULL, &tstat))
+                *rot=0.;
+        }
+    }
+    else  /* Case 2: CDELTn + optional PC matrix */
+    {
+        if (ffgkyd(fptr, "CDELT2", yinc, NULL, &tstat))
+            *yinc = 1.;
+
+        tstat = 0;
+        if (ffgkyd(fptr, "CROTA2", rot, NULL, &tstat))
+        {
+            *rot=0.;
+
+            /* no CROTA2 keyword, so look for the PC matrix */
+            tstat = 0;
+            if (ffgkyd(fptr, "PC1_1", &pc11, NULL, &tstat))
+                tstat = 0;  /* reset keyword not found error */
+            else
+                pc_exists = 1;  /* found at least 1 PC_ keyword */
+
+            if (ffgkyd(fptr, "PC2_1", &pc21, NULL, &tstat))
+                tstat = 0;  /* reset keyword not found error */
+            else
+                pc_exists = 1;  /* found at least 1 PC_ keyword */
+
+            if (ffgkyd(fptr, "PC1_2", &pc12, NULL, &tstat))
+                tstat = 0;  /* reset keyword not found error */
+            else
+                pc_exists = 1;  /* found at least 1 PC_ keyword */
+
+            if (ffgkyd(fptr, "PC2_2", &pc22, NULL, &tstat))
+                tstat = 0;  /* reset keyword not found error */
+            else
+                pc_exists = 1;  /* found at least 1 PC_ keyword */
+
+            if (pc_exists)  /* convert PCi_j back to CDELTn */
+            {
+                /* there are 2 ways to compute the angle: */
+                phia = atan2( pc21, pc11);
+                phib = atan2(-pc12, pc22);
+
+                /* ensure that phia <= phib */
+                temp = minvalue(phia, phib);
+                phib = maxvalue(phia, phib);
+                phia = temp;
+
+                /* there is a possible 180 degree ambiguity in the angles */
+                /* so add 180 degress to the smaller value if the values  */
+                /* differ by more than 90 degrees = pi/2 radians.         */
+                /* (Later, we may decide to take the other solution by    */
+                /* subtracting 180 degrees from the larger value).        */
+
+                if ((phib - phia) > (pi / 2.))
+                   phia += pi;
+
+                if (fabs(phia - phib) > toler) 
+                {
+                  /* angles don't agree, so looks like there is some skewness */
+                  /* between the axes.  Return with an error to be safe. */
+                  *status = APPROX_WCS_KEY;
+                }
+      
+                phia = (phia + phib) /2.;  /* use the average of the 2 values */
+                *rot = phia * 180. / pi;
+            }
+        }
+    }
+
+    /* get the type of projection, if any */
+    tstat = 0;
+    if (ffgkys(fptr, "CTYPE1", ctype, NULL, &tstat))
+         type[0] = '\0';
+    else
+    {
+        /* copy the projection type string */
+        strncpy(type, &ctype[4], 4);
+        type[4] = '\0';
+
+        /* check if RA and DEC are inverted */
+        if (!strncmp(ctype, "DEC-", 4) || !strncmp(ctype+1, "LAT", 3))
+        {
+            /* the latitudinal axis is given first, so swap them */
+
+/*
+ this case was removed on 12/9.  Apparently not correct.
+
+            if ((*xinc / *yinc) < 0. )  
+                *rot = -90. - (*rot);
+            else
+*/
+            *rot = 90. - (*rot);
+
+            /* Empirical tests with ds9 show the y-axis sign must be negated */
+            /* and the xinc and yinc values must NOT be swapped. */
+            *yinc = -(*yinc);
+
+            temp = *xrval;
+            *xrval = *yrval;
+            *yrval = temp;
+        }   
+    }
+
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int ffgtcs(fitsfile *fptr,    /* I - FITS file pointer           */
+           int xcol,          /* I - column containing the RA coordinate  */
+           int ycol,          /* I - column containing the DEC coordinate */
+           double *xrval,     /* O - X reference value           */
+           double *yrval,     /* O - Y reference value           */
+           double *xrpix,     /* O - X reference pixel           */
+           double *yrpix,     /* O - Y reference pixel           */
+           double *xinc,      /* O - X increment per pixel       */
+           double *yinc,      /* O - Y increment per pixel       */
+           double *rot,       /* O - rotation angle (degrees)    */
+           char *type,        /* O - type of projection ('-sin') */
+           int *status)       /* IO - error status               */
+/*
+       read the values of the celestial coordinate system keywords
+       from a FITS table where the X and Y or RA and DEC coordinates
+       are stored in separate column.  
+       These values may be used as input to the subroutines that
+       calculate celestial coordinates. (ffxypx, ffwldp)
+*/
+{
+    char comm[FLEN_COMMENT],ctype[FLEN_VALUE],keynam[FLEN_KEYWORD];
+    int tstatus = 0;
+
+    if (*status > 0)
+       return(*status);
+
+    ffkeyn("TCRVL",xcol,keynam,status);
+    ffgkyd(fptr,keynam,xrval,comm,status);
+
+    ffkeyn("TCRVL",ycol,keynam,status);
+    ffgkyd(fptr,keynam,yrval,comm,status);
+
+    ffkeyn("TCRPX",xcol,keynam,status);
+    ffgkyd(fptr,keynam,xrpix,comm,status);
+
+    ffkeyn("TCRPX",ycol,keynam,status);
+    ffgkyd(fptr,keynam,yrpix,comm,status);
+
+    ffkeyn("TCDLT",xcol,keynam,status);
+    ffgkyd(fptr,keynam,xinc,comm,status);
+
+    ffkeyn("TCDLT",ycol,keynam,status);
+    ffgkyd(fptr,keynam,yinc,comm,status);
+
+    ffkeyn("TCTYP",xcol,keynam,status);
+    ffgkys(fptr,keynam,ctype,comm,status);
+
+    if (*status > 0)
+    {
+      ffpmsg
+      ("ffgtcs could not find all the celestial coordinate keywords");
+      return(*status = NO_WCS_KEY); 
+    }
+
+    /* copy the projection type string */
+    strncpy(type, &ctype[4], 4);
+    type[4] = '\0';
+
+    *rot=0.;   /* default rotation is 0  */
+    ffkeyn("TCROT",ycol,keynam,status);
+    ffgkyd(fptr,keynam,rot,comm,&tstatus);  /* keyword may not exist */
+
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int ffwldp(double xpix, double ypix, double xref, double yref,
+      double xrefpix, double yrefpix, double xinc, double yinc, double rot,
+      char *type, double *xpos, double *ypos, int *status)
+
+/* WDP 1/97: change the name of the routine from 'worldpos' to 'ffwldp' */
+
+/*  worldpos.c -- WCS Algorithms from Classic AIPS.
+    Copyright (C) 1994
+    Associated Universities, Inc. Washington DC, USA.
+   
+    This library is free software; you can redistribute it and/or modify it
+    under the terms of the GNU Library 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 Library General Public
+    License for more details.
+   
+    You should have received a copy of the GNU Library General Public License
+    along with this library; if not, write to the Free Software Foundation,
+    Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
+   
+    Correspondence concerning AIPS should be addressed as follows:
+           Internet email: aipsmail@nrao.edu
+           Postal address: AIPS Group
+                           National Radio Astronomy Observatory
+                           520 Edgemont Road
+                           Charlottesville, VA 22903-2475 USA
+
+                 -=-=-=-=-=-=-
+
+    These two ANSI C functions, worldpos() and xypix(), perform
+    forward and reverse WCS computations for 8 types of projective
+    geometries ("-SIN", "-TAN", "-ARC", "-NCP", "-GLS", "-MER", "-AIT"
+    and "-STG"):
+
+        worldpos() converts from pixel location to RA,Dec 
+        xypix()    converts from RA,Dec         to pixel location   
+
+    where "(RA,Dec)" are more generically (long,lat). These functions
+    are based on the WCS implementation of Classic AIPS, an
+    implementation which has been in production use for more than ten
+    years. See the two memos by Eric Greisen
+
+        ftp://fits.cv.nrao.edu/fits/documents/wcs/aips27.ps.Z
+	ftp://fits.cv.nrao.edu/fits/documents/wcs/aips46.ps.Z
+
+    for descriptions of the 8 projective geometries and the
+    algorithms.  Footnotes in these two documents describe the
+    differences between these algorithms and the 1993-94 WCS draft
+    proposal (see URL below). In particular, these algorithms support
+    ordinary field rotation, but not skew geometries (CD or PC matrix
+    cases). Also, the MER and AIT algorithms work correctly only for
+    CRVALi=(0,0). Users should note that GLS projections with yref!=0
+    will behave differently in this code than in the draft WCS
+    proposal.  The NCP projection is now obsolete (it is a special
+    case of SIN).  WCS syntax and semantics for various advanced
+    features is discussed in the draft WCS proposal by Greisen and
+    Calabretta at:
+    
+        ftp://fits.cv.nrao.edu/fits/documents/wcs/wcs.all.ps.Z
+    
+                -=-=-=-
+
+    The original version of this code was Emailed to D.Wells on
+    Friday, 23 September by Bill Cotton <bcotton@gorilla.cv.nrao.edu>,
+    who described it as a "..more or less.. exact translation from the
+    AIPSish..". Changes were made by Don Wells <dwells@nrao.edu>
+    during the period October 11-13, 1994:
+    1) added GNU license and header comments
+    2) added testpos.c program to perform extensive circularity tests
+    3) changed float-->double to get more than 7 significant figures
+    4) testpos.c circularity test failed on MER and AIT. B.Cotton
+       found that "..there were a couple of lines of code [in] the wrong
+       place as a result of merging several Fortran routines." 
+    5) testpos.c found 0h wraparound in xypix() and worldpos().
+    6) E.Greisen recommended removal of various redundant if-statements,
+       and addition of a 360d difference test to MER case of worldpos(). 
+*/
+
+/*-----------------------------------------------------------------------*/
+/* routine to determine accurate position for pixel coordinates          */
+/* returns 0 if successful otherwise:                                    */
+/* 1 = angle too large for projection;                                   */
+/* (WDP 1/97: changed the return value to 501 instead of 1)              */
+/* does: -SIN, -TAN, -ARC, -NCP, -GLS, -MER, -AIT projections            */
+/* anything else is linear (== -CAR)                                     */
+/* Input:                                                                */
+/*   f   xpix    x pixel number  (RA or long without rotation)           */
+/*   f   ypiy    y pixel number  (dec or lat without rotation)           */
+/*   d   xref    x reference coordinate value (deg)                      */
+/*   d   yref    y reference coordinate value (deg)                      */
+/*   f   xrefpix x reference pixel                                       */
+/*   f   yrefpix y reference pixel                                       */
+/*   f   xinc    x coordinate increment (deg)                            */
+/*   f   yinc    y coordinate increment (deg)                            */
+/*   f   rot     rotation (deg)  (from N through E)                      */
+/*   c  *type    projection type code e.g. "-SIN";                       */
+/* Output:                                                               */
+/*   d   *xpos   x (RA) coordinate (deg)                                 */
+/*   d   *ypos   y (dec) coordinate (deg)                                */
+/*-----------------------------------------------------------------------*/
+ {double cosr, sinr, dx, dy, dz, temp, x, y, z;
+  double sins, coss, dect, rat, dt, l, m, mg, da, dd, cos0, sin0;
+  double dec0, ra0, decout, raout;
+  double geo1, geo2, geo3;
+  double cond2r=1.745329252e-2;
+  double twopi = 6.28318530717959, deps = 1.0e-5;
+  int   i, itype;
+  char ctypes[9][5] ={"-CAR","-SIN","-TAN","-ARC","-NCP", "-GLS", "-MER",
+     "-AIT", "-STG"};
+
+  if (*status > 0)
+     return(*status);
+
+/*   Offset from ref pixel  */
+  dx = (xpix-xrefpix) * xinc;
+  dy = (ypix-yrefpix) * yinc;
+/*   Take out rotation  */
+  cosr = cos(rot*cond2r);
+  sinr = sin(rot*cond2r);
+  if (rot!=0.0)
+    {temp = dx * cosr - dy * sinr;
+     dy = dy * cosr + dx * sinr;
+     dx = temp;}
+/*  find type  */
+/* WDP 1/97: removed support for default type for better error checking */
+/*  itype = 0;   default type is linear */
+  itype = -1;  /* no default type */
+  for (i=0;i<9;i++) if (!strncmp(type, ctypes[i], 4)) itype = i;
+/* default, linear result for error return  */
+  *xpos = xref + dx;
+  *ypos = yref + dy;
+/* convert to radians  */
+  ra0 = xref * cond2r;
+  dec0 = yref * cond2r;
+  l = dx * cond2r;
+  m = dy * cond2r;
+  sins = l*l + m*m;
+  cos0 = cos(dec0);
+  sin0 = sin(dec0);
+
+/* process by case  */
+  switch (itype) {
+    case 0:   /* linear -CAR */
+      rat =  ra0 + l;
+      dect = dec0 + m;
+      break;
+    case 1:   /* -SIN sin*/ 
+      if (sins>1.0) return(*status = 501);
+      coss = sqrt (1.0 - sins);
+      dt = sin0 * coss + cos0 * m;
+      if ((dt>1.0) || (dt<-1.0)) return(*status = 501);
+      dect = asin (dt);
+      rat = cos0 * coss - sin0 * m;
+      if ((rat==0.0) && (l==0.0)) return(*status = 501);
+      rat = atan2 (l, rat) + ra0;
+      break;
+    case 2:   /* -TAN tan */
+      x = cos0*cos(ra0) - l*sin(ra0) - m*cos(ra0)*sin0;
+      y = cos0*sin(ra0) + l*cos(ra0) - m*sin(ra0)*sin0;
+      z = sin0                       + m*         cos0;
+      rat  = atan2( y, x );
+      dect = atan ( z / sqrt(x*x+y*y) );
+      break;
+    case 3:   /* -ARC Arc*/
+      if (sins>=twopi*twopi/4.0) return(*status = 501);
+      sins = sqrt(sins);
+      coss = cos (sins);
+      if (sins!=0.0) sins = sin (sins) / sins;
+      else
+	sins = 1.0;
+      dt = m * cos0 * sins + sin0 * coss;
+      if ((dt>1.0) || (dt<-1.0)) return(*status = 501);
+      dect = asin (dt);
+      da = coss - dt * sin0;
+      dt = l * sins * cos0;
+      if ((da==0.0) && (dt==0.0)) return(*status = 501);
+      rat = ra0 + atan2 (dt, da);
+      break;
+    case 4:   /* -NCP North celestial pole*/
+      dect = cos0 - m * sin0;
+      if (dect==0.0) return(*status = 501);
+      rat = ra0 + atan2 (l, dect);
+      dt = cos (rat-ra0);
+      if (dt==0.0) return(*status = 501);
+      dect = dect / dt;
+      if ((dect>1.0) || (dect<-1.0)) return(*status = 501);
+      dect = acos (dect);
+      if (dec0<0.0) dect = -dect;
+      break;
+    case 5:   /* -GLS global sinusoid */
+      dect = dec0 + m;
+      if (fabs(dect)>twopi/4.0) return(*status = 501);
+      coss = cos (dect);
+      if (fabs(l)>twopi*coss/2.0) return(*status = 501);
+      rat = ra0;
+      if (coss>deps) rat = rat + l / coss;
+      break;
+    case 6:   /* -MER mercator*/
+      dt = yinc * cosr + xinc * sinr;
+      if (dt==0.0) dt = 1.0;
+      dy = (yref/2.0 + 45.0) * cond2r;
+      dx = dy + dt / 2.0 * cond2r;
+      dy = log (tan (dy));
+      dx = log (tan (dx));
+      geo2 = dt * cond2r / (dx - dy);
+      geo3 = geo2 * dy;
+      geo1 = cos (yref*cond2r);
+      if (geo1<=0.0) geo1 = 1.0;
+      rat = l / geo1 + ra0;
+      if (fabs(rat - ra0) > twopi) return(*status = 501); /* added 10/13/94 DCW/EWG */
+      dt = 0.0;
+      if (geo2!=0.0) dt = (m + geo3) / geo2;
+      dt = exp (dt);
+      dect = 2.0 * atan (dt) - twopi / 4.0;
+      break;
+    case 7:   /* -AIT Aitoff*/
+      dt = yinc*cosr + xinc*sinr;
+      if (dt==0.0) dt = 1.0;
+      dt = dt * cond2r;
+      dy = yref * cond2r;
+      dx = sin(dy+dt)/sqrt((1.0+cos(dy+dt))/2.0) -
+	  sin(dy)/sqrt((1.0+cos(dy))/2.0);
+      if (dx==0.0) dx = 1.0;
+      geo2 = dt / dx;
+      dt = xinc*cosr - yinc* sinr;
+      if (dt==0.0) dt = 1.0;
+      dt = dt * cond2r;
+      dx = 2.0 * cos(dy) * sin(dt/2.0);
+      if (dx==0.0) dx = 1.0;
+      geo1 = dt * sqrt((1.0+cos(dy)*cos(dt/2.0))/2.0) / dx;
+      geo3 = geo2 * sin(dy) / sqrt((1.0+cos(dy))/2.0);
+      rat = ra0;
+      dect = dec0;
+      if ((l==0.0) && (m==0.0)) break;
+      dz = 4.0 - l*l/(4.0*geo1*geo1) - ((m+geo3)/geo2)*((m+geo3)/geo2) ;
+      if ((dz>4.0) || (dz<2.0)) return(*status = 501);;
+      dz = 0.5 * sqrt (dz);
+      dd = (m+geo3) * dz / geo2;
+      if (fabs(dd)>1.0) return(*status = 501);;
+      dd = asin (dd);
+      if (fabs(cos(dd))<deps) return(*status = 501);;
+      da = l * dz / (2.0 * geo1 * cos(dd));
+      if (fabs(da)>1.0) return(*status = 501);;
+      da = asin (da);
+      rat = ra0 + 2.0 * da;
+      dect = dd;
+      break;
+    case 8:   /* -STG Sterographic*/
+      dz = (4.0 - sins) / (4.0 + sins);
+      if (fabs(dz)>1.0) return(*status = 501);
+      dect = dz * sin0 + m * cos0 * (1.0+dz) / 2.0;
+      if (fabs(dect)>1.0) return(*status = 501);
+      dect = asin (dect);
+      rat = cos(dect);
+      if (fabs(rat)<deps) return(*status = 501);
+      rat = l * (1.0+dz) / (2.0 * rat);
+      if (fabs(rat)>1.0) return(*status = 501);
+      rat = asin (rat);
+      mg = 1.0 + sin(dect) * sin0 + cos(dect) * cos0 * cos(rat);
+      if (fabs(mg)<deps) return(*status = 501);
+      mg = 2.0 * (sin(dect) * cos0 - cos(dect) * sin0 * cos(rat)) / mg;
+      if (fabs(mg-m)>deps) rat = twopi/2.0 - rat;
+      rat = ra0 + rat;
+      break;
+
+    default:
+      /* fall through to here on error */
+      return(*status = 504);
+  }
+
+/*  return ra in range  */
+  raout = rat;
+  decout = dect;
+  if (raout-ra0>twopi/2.0) raout = raout - twopi;
+  if (raout-ra0<-twopi/2.0) raout = raout + twopi;
+  if (raout < 0.0) raout += twopi; /* added by DCW 10/12/94 */
+
+/*  correct units back to degrees  */
+  *xpos  = raout  / cond2r;
+  *ypos  = decout  / cond2r;
+  return(*status);
+}  /* End of worldpos */
+/*--------------------------------------------------------------------------*/
+int ffxypx(double xpos, double ypos, double xref, double yref, 
+      double xrefpix, double yrefpix, double xinc, double yinc, double rot,
+      char *type, double *xpix, double *ypix, int *status)
+/* WDP  1/97: changed name of routine from xypix to ffxypx    */
+/*-----------------------------------------------------------------------*/
+/* routine to determine accurate pixel coordinates for an RA and Dec     */
+/* returns 0 if successful otherwise:                                    */
+/* 1 = angle too large for projection;                                   */
+/* 2 = bad values                                                        */
+/* WDP 1/97: changed the return values to 501 and 502 instead of 1 and 2 */
+/* does: -SIN, -TAN, -ARC, -NCP, -GLS, -MER, -AIT projections            */
+/* anything else is linear                                               */
+/* Input:                                                                */
+/*   d   xpos    x (RA) coordinate (deg)                                 */
+/*   d   ypos    y (dec) coordinate (deg)                                */
+/*   d   xref    x reference coordinate value (deg)                      */
+/*   d   yref    y reference coordinate value (deg)                      */
+/*   f   xrefpix x reference pixel                                       */
+/*   f   yrefpix y reference pixel                                       */
+/*   f   xinc    x coordinate increment (deg)                            */
+/*   f   yinc    y coordinate increment (deg)                            */
+/*   f   rot     rotation (deg)  (from N through E)                      */
+/*   c  *type    projection type code e.g. "-SIN";                       */
+/* Output:                                                               */
+/*   f  *xpix    x pixel number  (RA or long without rotation)           */
+/*   f  *ypiy    y pixel number  (dec or lat without rotation)           */
+/*-----------------------------------------------------------------------*/
+ {double dx, dy, dz, r, ra0, dec0, ra, dec, coss, sins, dt, da, dd, sint;
+  double l, m, geo1, geo2, geo3, sinr, cosr, cos0, sin0;
+  double cond2r=1.745329252e-2, deps=1.0e-5, twopi=6.28318530717959;
+  int   i, itype;
+  char ctypes[9][5] ={"-CAR","-SIN","-TAN","-ARC","-NCP", "-GLS", "-MER",
+     "-AIT", "-STG"};
+
+  /* 0h wrap-around tests added by D.Wells 10/12/94: */
+  dt = (xpos - xref);
+  if (dt >  180) xpos -= 360;
+  if (dt < -180) xpos += 360;
+  /* NOTE: changing input argument xpos is OK (call-by-value in C!) */
+
+/* default values - linear */
+  dx = xpos - xref;
+  dy = ypos - yref;
+/*  dz = 0.0; */
+/*  Correct for rotation */
+  r = rot * cond2r;
+  cosr = cos (r);
+  sinr = sin (r);
+  dz = dx*cosr + dy*sinr;
+  dy = dy*cosr - dx*sinr;
+  dx = dz;
+/*     check axis increments - bail out if either 0 */
+  if ((xinc==0.0) || (yinc==0.0)) {*xpix=0.0; *ypix=0.0; return(*status = 502);}
+/*     convert to pixels  */
+  *xpix = dx / xinc + xrefpix;
+  *ypix = dy / yinc + yrefpix;
+
+/*  find type  */
+/* WDP 1/97: removed support for default type for better error checking */
+/*  itype = 0;   default type is linear */
+  itype = -1;  /* no default type */
+  for (i=0;i<9;i++) if (!strncmp(type, ctypes[i], 4)) itype = i;
+  if (itype==0) return(*status);  /* done if linear */
+
+/* Non linear position */
+  ra0 = xref * cond2r;
+  dec0 = yref * cond2r;
+  ra = xpos * cond2r;
+  dec = ypos * cond2r;
+
+/* compute direction cosine */
+  coss = cos (dec);
+  sins = sin (dec);
+  cos0 = cos (dec0);
+  sin0 = sin (dec0);
+  l = sin(ra-ra0) * coss;
+  sint = sins * sin0 + coss * cos0 * cos(ra-ra0);
+
+/* process by case  */
+  switch (itype) {
+    case 1:   /* -SIN sin*/ 
+         if (sint<0.0) return(*status = 501);
+         m = sins * cos(dec0) - coss * sin(dec0) * cos(ra-ra0);
+      break;
+    case 2:   /* -TAN tan */
+         if (sint<=0.0) return(*status = 501);
+         if( cos0<0.001 ) {
+            /* Do a first order expansion around pole */
+            m = (coss * cos(ra-ra0)) / (sins * sin0);
+            m = (-m + cos0 * (1.0 + m*m)) / sin0;
+         } else {
+            m = ( sins/sint - sin0 ) / cos0;
+         }
+	 if( fabs(sin(ra0)) < 0.3 ) {
+	    l  = coss*sin(ra)/sint - cos0*sin(ra0) + m*sin(ra0)*sin0;
+	    l /= cos(ra0);
+	 } else {
+	    l  = coss*cos(ra)/sint - cos0*cos(ra0) + m*cos(ra0)*sin0;
+	    l /= -sin(ra0);
+	 }
+      break;
+    case 3:   /* -ARC Arc*/
+         m = sins * sin(dec0) + coss * cos(dec0) * cos(ra-ra0);
+         if (m<-1.0) m = -1.0;
+         if (m>1.0) m = 1.0;
+         m = acos (m);
+         if (m!=0) 
+            m = m / sin(m);
+         else
+            m = 1.0;
+         l = l * m;
+         m = (sins * cos(dec0) - coss * sin(dec0) * cos(ra-ra0)) * m;
+      break;
+    case 4:   /* -NCP North celestial pole*/
+         if (dec0==0.0) 
+	     return(*status = 501);  /* can't stand the equator */
+         else
+	   m = (cos(dec0) - coss * cos(ra-ra0)) / sin(dec0);
+      break;
+    case 5:   /* -GLS global sinusoid */
+         dt = ra - ra0;
+         if (fabs(dec)>twopi/4.0) return(*status = 501);
+         if (fabs(dec0)>twopi/4.0) return(*status = 501);
+         m = dec - dec0;
+         l = dt * coss;
+      break;
+    case 6:   /* -MER mercator*/
+         dt = yinc * cosr + xinc * sinr;
+         if (dt==0.0) dt = 1.0;
+         dy = (yref/2.0 + 45.0) * cond2r;
+         dx = dy + dt / 2.0 * cond2r;
+         dy = log (tan (dy));
+         dx = log (tan (dx));
+         geo2 = dt * cond2r / (dx - dy);
+         geo3 = geo2 * dy;
+         geo1 = cos (yref*cond2r);
+         if (geo1<=0.0) geo1 = 1.0;
+         dt = ra - ra0;
+         l = geo1 * dt;
+         dt = dec / 2.0 + twopi / 8.0;
+         dt = tan (dt);
+         if (dt<deps) return(*status = 502);
+         m = geo2 * log (dt) - geo3;
+         break;
+    case 7:   /* -AIT Aitoff*/
+         da = (ra - ra0) / 2.0;
+         if (fabs(da)>twopi/4.0) return(*status = 501);
+         dt = yinc*cosr + xinc*sinr;
+         if (dt==0.0) dt = 1.0;
+         dt = dt * cond2r;
+         dy = yref * cond2r;
+         dx = sin(dy+dt)/sqrt((1.0+cos(dy+dt))/2.0) -
+             sin(dy)/sqrt((1.0+cos(dy))/2.0);
+         if (dx==0.0) dx = 1.0;
+         geo2 = dt / dx;
+         dt = xinc*cosr - yinc* sinr;
+         if (dt==0.0) dt = 1.0;
+         dt = dt * cond2r;
+         dx = 2.0 * cos(dy) * sin(dt/2.0);
+         if (dx==0.0) dx = 1.0;
+         geo1 = dt * sqrt((1.0+cos(dy)*cos(dt/2.0))/2.0) / dx;
+         geo3 = geo2 * sin(dy) / sqrt((1.0+cos(dy))/2.0);
+         dt = sqrt ((1.0 + cos(dec) * cos(da))/2.0);
+         if (fabs(dt)<deps) return(*status = 503);
+         l = 2.0 * geo1 * cos(dec) * sin(da) / dt;
+         m = geo2 * sin(dec) / dt - geo3;
+      break;
+    case 8:   /* -STG Sterographic*/
+         da = ra - ra0;
+         if (fabs(dec)>twopi/4.0) return(*status = 501);
+         dd = 1.0 + sins * sin(dec0) + coss * cos(dec0) * cos(da);
+         if (fabs(dd)<deps) return(*status = 501);
+         dd = 2.0 / dd;
+         l = l * dd;
+         m = dd * (sins * cos(dec0) - coss * sin(dec0) * cos(da));
+      break;
+
+    default:
+      /* fall through to here on error */
+      return(*status = 504);
+
+  }  /* end of itype switch */
+
+/*   back to degrees  */
+  dx = l / cond2r;
+  dy = m / cond2r;
+/*  Correct for rotation */
+  dz = dx*cosr + dy*sinr;
+  dy = dy*cosr - dx*sinr;
+  dx = dz;
+/*     convert to pixels  */
+  *xpix = dx / xinc + xrefpix;
+  *ypix = dy / yinc + yrefpix;
+  return(*status);
+}  /* end xypix */