rm funtools for update

1 files changed, 0 insertions, 693 deletions

diff --git a/funtools/wcs/worldpos.c b/funtools/wcs/worldpos.c
deleted file mode 100644
index 8693789..0000000
--- a/funtools/wcs/worldpos.c
+++ /dev/null
@@ -1,693 +0,0 @@
-/*  worldpos.c -- WCS Algorithms from Classic AIPS.
- *  September 1, 2011
- *  Copyright (C) 1994-2011
- *  Associated Universities, Inc. Washington DC, USA.
- *  With code added by Jessica Mink, Smithsonian Astrophysical Observatory
- *                 and Allan Brighton and Andreas Wicenec, ESO
- *                 and Frank Valdes, NOAO
-
- * Module:	worldpos.c
- * Purpose:	Perform forward and reverse WCS computations for 8 projections
- * Subroutine:	worldpos() converts from pixel location to RA,Dec 
- * Subroutine:	worldpix() converts from RA,Dec         to pixel location   
-
-    -=-=-=-=-=-=-
-
-    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 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 worldpix(), perform
-    forward and reverse WCS computations for 8 types of projective
-    geometries ("-SIN", "-TAN", "-ARC", "-NCP", "-GLS" or "-SFL", "-MER",
-     "-AIT", "-STG", "CAR", and "COE"):
-
-	worldpos() converts from pixel location to RA,Dec 
-	worldpix() 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 worldpix() and worldpos().
-    6) E.Greisen recommended removal of various redundant if-statements,
-       and addition of a 360d difference test to MER case of worldpos(). 
-    7) D.Mink changed input to data structure and implemented rotation matrix.
-*/
-#include <math.h>
-#include <string.h>
-#include <stdio.h>
-#include "wcs.h"
-
-int
-worldpos (xpix, ypix, wcs, xpos, ypos)
-
-/* Routine to determine accurate position for pixel coordinates */
-/* returns 0 if successful otherwise 1 = angle too large for projection; */
-/* does: -SIN, -TAN, -ARC, -NCP, -GLS or -SFL, -MER, -AIT projections */
-/* anything else is linear */
-
-/* Input: */
-double	xpix;		/* x pixel number  (RA or long without rotation) */
-double	ypix;		/* y pixel number  (Dec or lat without rotation) */
-struct WorldCoor *wcs;		/* WCS parameter structure */
-
-/* Output: */
-double	*xpos;		/* x (RA) coordinate (deg) */
-double	*ypos;		/* y (dec) coordinate (deg) */
-
-{
-  double cosr, sinr, dx, dy, dz, tx;
-  double sins, coss, dt, l, m, mg, da, dd, cos0, sin0;
-  double rat = 0.0;
-  double dect = 0.0;
-  double mt, a, y0, td, r2;  /* allan: for COE */
-  double dec0, ra0, decout, raout;
-  double geo1, geo2, geo3;
-  double cond2r=1.745329252e-2;
-  double twopi = 6.28318530717959;
-  double deps = 1.0e-5;
-
-  /* Structure elements */
-  double xref;		/* X reference coordinate value (deg) */
-  double yref;		/* Y reference coordinate value (deg) */
-  double xrefpix;	/* X reference pixel */
-  double yrefpix;	/* Y reference pixel */
-  double xinc;		/* X coordinate increment (deg) */
-  double yinc;		/* Y coordinate increment (deg) */
-  double rot;		/* Optical axis rotation (deg)  (N through E) */
-  int itype = wcs->prjcode;
-
-  /* Set local projection parameters */
-  xref = wcs->xref;
-  yref = wcs->yref;
-  xrefpix = wcs->xrefpix;
-  yrefpix = wcs->yrefpix;
-  xinc = wcs->xinc;
-  yinc = wcs->yinc;
-  rot = degrad (wcs->rot);
-  cosr = cos (rot);
-  sinr = sin (rot);
-
-  /* Offset from ref pixel */
-  dx = xpix - xrefpix;
-  dy = ypix - yrefpix;
-
-  /* Scale and rotate using CD matrix */
-  if (wcs->rotmat) {
-    tx = dx * wcs->cd[0] + dy * wcs->cd[1];
-    dy = dx * wcs->cd[2] + dy * wcs->cd[3];
-    dx = tx;
-    }
-
-  /* Scale and rotate using CDELTn and CROTA2 */
-  else {
-
-    /* Check axis increments - bail out if either 0 */
-    if ((xinc==0.0) || (yinc==0.0)) {
-      *xpos=0.0;
-      *ypos=0.0;
-      return 2;
-      }
-
-    /* Scale using CDELT */
-    dx = dx * xinc;
-    dy = dy * yinc;
-
-    /* Take out rotation from CROTA */
-    if (rot != 0.0) {
-      tx = dx * cosr - dy * sinr;
-      dy = dx * sinr + dy * cosr;
-      dx = tx;
-      }
-    }
-
-  /* Flip coordinates if necessary */
-  if (wcs->coorflip) {
-    tx = dx;
-    dx = dy;
-    dy = tx;
-    }
-
-  /* Default, linear result for error or pixel return  */
-  *xpos = xref + dx;
-  *ypos = yref + dy;
-  if (itype <= 0)
-    return 0;
-
-  /* Convert to radians  */
-  if (wcs->coorflip) {
-    dec0 = degrad (xref);
-    ra0 = degrad (yref);
-    }
-  else {
-    ra0 = degrad (xref);
-    dec0 = degrad (yref);
-    }
-  l = degrad (dx);
-  m = degrad (dy);
-  sins = l*l + m*m;
-  decout = 0.0;
-  raout = 0.0;
-  cos0 = cos (dec0);
-  sin0 = sin (dec0);
-
-  /* Process by case  */
-  switch (itype) {
-
-    case WCS_CAR:   /* -CAR Cartesian (was WCS_PIX pixel and WCS_LIN linear) */
-      rat =  ra0 + l;
-      dect = dec0 + m;
-      break;
-
-    case WCS_SIN: /* -SIN sin*/ 
-      if (sins>1.0) return 1;
-      coss = sqrt (1.0 - sins);
-      dt = sin0 * coss + cos0 * m;
-      if ((dt>1.0) || (dt<-1.0)) return 1;
-      dect = asin (dt);
-      rat = cos0 * coss - sin0 * m;
-      if ((rat==0.0) && (l==0.0)) return 1;
-      rat = atan2 (l, rat) + ra0;
-      break;
-
-    case WCS_TAN:   /* -TAN tan */
-    case WCS_TNX:   /* -TNX tan with polynomial correction */
-    case WCS_TPV:   /* -TPV tan with polynomial correction */
-    case WCS_ZPX:   /* -ZPX zpn with polynomial correction */
-      if (sins>1.0) return 1;
-      dect = cos0 - m * sin0;
-      if (dect==0.0) return 1;
-      rat = ra0 + atan2 (l, dect);
-      dect = atan (cos(rat-ra0) * (m * cos0 + sin0) / dect);
-      break;
-
-    case WCS_ARC:   /* -ARC Arc*/
-      if (sins>=twopi*twopi/4.0) return 1;
-      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 1;
-      dect = asin (dt);
-      da = coss - dt * sin0;
-      dt = l * sins * cos0;
-      if ((da==0.0) && (dt==0.0)) return 1;
-      rat = ra0 + atan2 (dt, da);
-      break;
-
-    case WCS_NCP:   /* -NCP North celestial pole*/
-      dect = cos0 - m * sin0;
-      if (dect==0.0) return 1;
-      rat = ra0 + atan2 (l, dect);
-      dt = cos (rat-ra0);
-      if (dt==0.0) return 1;
-      dect = dect / dt;
-      if ((dect>1.0) || (dect<-1.0)) return 1;
-      dect = acos (dect);
-      if (dec0<0.0) dect = -dect;
-      break;
-
-    case WCS_GLS:   /* -GLS global sinusoid */
-    case WCS_SFL:   /* -SFL Samson-Flamsteed */
-      dect = dec0 + m;
-      if (fabs(dect)>twopi/4.0) return 1;
-      coss = cos (dect);
-      if (fabs(l)>twopi*coss/2.0) return 1;
-      rat = ra0;
-      if (coss>deps) rat = rat + l / coss;
-      break;
-
-    case WCS_MER:   /* -MER mercator*/
-      dt = yinc * cosr + xinc * sinr;
-      if (dt==0.0) dt = 1.0;
-      dy = degrad (yref/2.0 + 45.0);
-      dx = dy + dt / 2.0 * cond2r;
-      dy = log (tan (dy));
-      dx = log (tan (dx));
-      geo2 = degrad (dt) / (dx - dy);
-      geo3 = geo2 * dy;
-      geo1 = cos (degrad (yref));
-      if (geo1<=0.0) geo1 = 1.0;
-      rat = l / geo1 + ra0;
-      if (fabs(rat - ra0) > twopi) return 1; /* 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 WCS_AIT:   /* -AIT Aitoff*/
-      dt = yinc*cosr + xinc*sinr;
-      if (dt==0.0) dt = 1.0;
-      dt = degrad (dt);
-      dy = degrad (yref);
-      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 = degrad (dt);
-      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 1;;
-      dz = 0.5 * sqrt (dz);
-      dd = (m+geo3) * dz / geo2;
-      if (fabs(dd)>1.0) return 1;;
-      dd = asin (dd);
-      if (fabs(cos(dd))<deps) return 1;;
-      da = l * dz / (2.0 * geo1 * cos(dd));
-      if (fabs(da)>1.0) return 1;;
-      da = asin (da);
-      rat = ra0 + 2.0 * da;
-      dect = dd;
-      break;
-
-    case WCS_STG:   /* -STG Sterographic*/
-      dz = (4.0 - sins) / (4.0 + sins);
-      if (fabs(dz)>1.0) return 1;
-      dect = dz * sin0 + m * cos0 * (1.0+dz) / 2.0;
-      if (fabs(dect)>1.0) return 1;
-      dect = asin (dect);
-      rat = cos(dect);
-      if (fabs(rat)<deps) return 1;
-      rat = l * (1.0+dz) / (2.0 * rat);
-      if (fabs(rat)>1.0) return 1;
-      rat = asin (rat);
-      mg = 1.0 + sin(dect) * sin0 + cos(dect) * cos0 * cos(rat);
-      if (fabs(mg)<deps) return 1;
-      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;
-
-    case WCS_COE:    /* COE projection code from Andreas Wicenic, ESO */
-      td = tan (dec0);
-      y0 = 1.0 / td;
-      mt = y0 - m;
-      if (dec0 < 0.)
-	a = atan2 (l,-mt);
-      else
-	a = atan2 (l, mt);
-      rat = ra0 - (a / sin0);
-      r2 = (l * l) + (mt * mt);
-      dect = asin (1.0 / (sin0 * 2.0) * (1.0 + sin0*sin0 * (1.0 - r2)));
-      break;
-  }
-
-  /* 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 */
-
-  /* Convert units back to degrees  */
-  *xpos = raddeg (raout);
-  *ypos = raddeg (decout);
-
-  return 0;
-}  /* End of worldpos */
-
-
-int
-worldpix (xpos, ypos, wcs, xpix, ypix)
-
-/*-----------------------------------------------------------------------*/
-/* 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                                                       */
-/* does: SIN, TAN, ARC, NCP, GLS or SFL, MER, AIT, STG, CAR, COE projections    */
-/* anything else is linear                                               */
-
-/* Input: */
-double	xpos;		/* x (RA) coordinate (deg) */
-double	ypos;		/* y (dec) coordinate (deg) */
-struct WorldCoor *wcs;	/* WCS parameter structure */
-
-/* Output: */
-double	*xpix;		/* x pixel number  (RA or long without rotation) */
-double	*ypix;		/* y pixel number  (dec or lat without rotation) */
-{
-  double dx, dy, ra0, dec0, ra, dec, coss, sins, dt, da, dd, sint;
-  double l, m, geo1, geo2, geo3, sinr, cosr, tx, x, a2, a3, a4;
-  double rthea,gamby2,a,b,c,phi,an,rap,v,tthea,co1,co2,co3,co4,ansq; /* COE */
-  double cond2r=1.745329252e-2, deps=1.0e-5, twopi=6.28318530717959;
-
-/* Structure elements */
-  double xref;		/* x reference coordinate value (deg) */
-  double yref;		/* y reference coordinate value (deg) */
-  double xrefpix;	/* x reference pixel */
-  double yrefpix;	/* y reference pixel */
-  double xinc;		/* x coordinate increment (deg) */
-  double yinc;		/* y coordinate increment (deg) */
-  double rot;		/* Optical axis rotation (deg)  (from N through E) */
-  int itype;
-
-  /* Set local projection parameters */
-  xref = wcs->xref;
-  yref = wcs->yref;
-  xrefpix = wcs->xrefpix;
-  yrefpix = wcs->yrefpix;
-  xinc = wcs->xinc;
-  yinc = wcs->yinc;
-  rot = degrad (wcs->rot);
-  cosr = cos (rot);
-  sinr = sin (rot);
-
-  /* Projection type */
-  itype = wcs->prjcode;
-
-  /* Nonlinear position */
-  if (itype > 0) {
-    if (wcs->coorflip) {
-      dec0 = degrad (xref);
-      ra0 = degrad (yref);
-      dt = xpos - yref;
-      }
-    else {
-      ra0 = degrad (xref);
-      dec0 = degrad (yref);
-      dt = xpos - xref;
-      }
-
-    /* 0h wrap-around tests added by D.Wells 10/12/1994: */
-    /* Modified to exclude weird reference pixels by D.Mink 2/3/2004 */
-    if (xrefpix*xinc > 180.0 || xrefpix*xinc < -180.0) {
-	if (dt > 360.0) xpos -= 360.0;
-	if (dt < 0.0) xpos += 360.0;
-	}
-    else {
-	if (dt > 180.0) xpos -= 360.0;
-	if (dt < -180.0) xpos += 360.0;
-	}
-    /* NOTE: changing input argument xpos is OK (call-by-value in C!) */
-
-    ra = degrad (xpos);
-    dec = degrad (ypos);
-
-    /* Compute direction cosine */
-    coss = cos (dec);
-    sins = sin (dec);
-    l = sin(ra-ra0) * coss;
-    sint = sins * sin(dec0) + coss * cos(dec0) * cos(ra-ra0);
-    }
-  else {
-    l = 0.0;
-    sint = 0.0;
-    sins = 0.0;
-    coss = 0.0;
-    ra = 0.0;
-    dec = 0.0;
-    ra0 = 0.0;
-    dec0 = 0.0;
-    m = 0.0;
-    }
-
-  /* Process by case  */
-  switch (itype) {
-
-    case WCS_CAR:   /* -CAR Cartesian */
-      l = ra - ra0;
-      m = dec - dec0;
-      break;
-
-    case WCS_SIN:   /* -SIN sin*/ 
-	if (sint<0.0) return 1;
-	m = sins * cos(dec0) - coss * sin(dec0) * cos(ra-ra0);
-	break;
-
-    case WCS_TNX:   /* -TNX tan with polynomial correction */
-    case WCS_TPV:   /* -TPV tan with polynomial correction */
-    case WCS_ZPX:   /* -ZPX zpn with polynomial correction */
-    case WCS_TAN:   /* -TAN tan */
-	if (sint<=0.0) return 1;
- 	m = sins * sin(dec0) + coss * cos(dec0) * cos(ra-ra0);
-	l = l / m;
-	m = (sins * cos(dec0) - coss * sin(dec0) * cos(ra-ra0)) / m;
-	break;
-
-    case WCS_ARC:   /* -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 WCS_NCP:   /* -NCP North celestial pole*/
-	if (dec0==0.0) 
-	    return 1;  /* can't stand the equator */
-	else
-	    m = (cos(dec0) - coss * cos(ra-ra0)) / sin(dec0);
-	break;
-
-    case WCS_GLS:   /* -GLS global sinusoid */
-    case WCS_SFL:   /* -SFL Samson-Flamsteed */
-	dt = ra - ra0;
-	if (fabs(dec)>twopi/4.0) return 1;
-	if (fabs(dec0)>twopi/4.0) return 1;
-	m = dec - dec0;
-	l = dt * coss;
-	break;
-
-    case WCS_MER:   /* -MER mercator*/
-	dt = yinc * cosr + xinc * sinr;
-	if (dt==0.0) dt = 1.0;
-	dy = degrad (yref/2.0 + 45.0);
-	dx = dy + dt / 2.0 * cond2r;
-	dy = log (tan (dy));
-	dx = log (tan (dx));
-	geo2 = degrad (dt) / (dx - dy);
-	geo3 = geo2 * dy;
-	geo1 = cos (degrad (yref));
-	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 2;
-	m = geo2 * log (dt) - geo3;
-	break;
-
-    case WCS_AIT:   /* -AIT Aitoff*/
-	l = 0.0;
-	m = 0.0;
-	da = (ra - ra0) / 2.0;
-	if (fabs(da)>twopi/4.0) return 1;
-	dt = yinc*cosr + xinc*sinr;
-	if (dt==0.0) dt = 1.0;
-	dt = degrad (dt);
-	dy = degrad (yref);
-	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 = degrad (dt);
-	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 3;
-	l = 2.0 * geo1 * cos(dec) * sin(da) / dt;
-	m = geo2 * sin(dec) / dt - geo3;
-	break;
-
-    case WCS_STG:   /* -STG Sterographic*/
-	da = ra - ra0;
-	if (fabs(dec)>twopi/4.0) return 1;
-	dd = 1.0 + sins * sin(dec0) + coss * cos(dec0) * cos(da);
-	if (fabs(dd)<deps) return 1;
-	dd = 2.0 / dd;
-	l = l * dd;
-	m = dd * (sins * cos(dec0) - coss * sin(dec0) * cos(da));
-	break;
-
-    case WCS_COE:    /* allan: -COE projection added, AW, ESO*/
-	gamby2 = sin (dec0);
-	tthea = tan (dec0);
-	rthea = 1. / tthea;
-	a = -2. * tthea;
-	b = tthea * tthea;
-	c = tthea / 3.;
-	a2 = a * a;
-	a3 = a2 * a;
-	a4 = a2 * a2;
-	co1 = a/2.;
-	co2 = -0.125 * a2 + b/2.;
-	co3 = -0.25 * a*b + 0.0625 * a3 + c/2.0;
-	co4 = -0.125 * b*b - 0.25 * a*c + 0.1875 * b*a2 - (5.0/128.0)*a4;
-	phi = ra0 - ra;
-	an = phi * gamby2;
-	v = dec - dec0;
-	rap = rthea * (1.0 + v * (co1+v * (co2+v * (co3+v * co4))));
-	ansq = an * an;
-	if (wcs->rotmat)
-	    l = rap * an * (1.0 - ansq/6.0) * (wcs->cd[0] / fabs(wcs->cd[0]));
-	else
-	    l = rap * an * (1.0 - ansq/6.0) * (xinc / fabs(xinc));
-	m = rthea - (rap * (1.0 - ansq/2.0));
-	break;
-
-    }  /* end of itype switch */
-
-  /* Convert back to degrees  */
-  if (itype > 0) {
-    dx = raddeg (l);
-    dy = raddeg (m);
-    }
-
-  /* For linear or pixel projection */
-  else {
-    dx = xpos - xref;
-    dy = ypos - yref;
-    }
-
-  if (wcs->coorflip) {
-    tx = dx;
-    dx = dy;
-    dy = tx;
-    }
-
-  /* Scale and rotate using CD matrix */
-  if (wcs->rotmat) {
-    tx = dx * wcs->dc[0] + dy * wcs->dc[1];
-    dy = dx * wcs->dc[2] + dy * wcs->dc[3];
-    dx = tx;
-    }
-
-  /* Scale and rotate using CDELTn and CROTA2 */
-  else {
-
-    /* Correct for rotation */
-    if (rot!=0.0) {
-      tx = dx*cosr + dy*sinr;
-      dy = dy*cosr - dx*sinr;
-      dx = tx;
-      }
-
-    /* Scale using CDELT */
-    if (xinc != 0.)
-      dx = dx / xinc;
-    if (yinc != 0.)
-      dy = dy / yinc;
-    }
-
-  /* Convert to pixels  */
-  *xpix = dx + xrefpix;
-  if (itype == WCS_CAR) {
-    if (*xpix > wcs->nxpix) {
-      x = *xpix - (360.0 / xinc);
-      if (x > 0.0) *xpix = x;
-      }
-    else if (*xpix < 0) {
-      x = *xpix + (360.0 / xinc);
-      if (x <= wcs->nxpix) *xpix = x;
-      }
-    }
-  *ypix = dy + yrefpix;
-
-  return 0;
-}  /* end worldpix */
-
- 
-/* Oct 26 1995	Fix bug which interchanged RA and Dec twice when coorflip
- *
- * Oct 31 1996	Fix CD matrix use in WORLDPIX
- * Nov  4 1996	Eliminate extra code for linear projection in WORLDPIX
- * Nov  5 1996	Add coordinate flip in WORLDPIX
- *
- * May 22 1997	Avoid angle wraparound when CTYPE is pixel
- * Jun  4 1997	Return without angle conversion from worldpos if type is PIXEL
- *
- * Oct 20 1997	Add chip rotation; compute rotation angle trig functions
- * Jan 23 1998	Change PCODE to PRJCODE
- * Jan 26 1998	Remove chip rotation code
- * Feb  5 1998	Make cd[] and dc[] vectors; use xinc, yinc, rot from init
- * Feb 23 1998	Add NOAO TNX projection as TAN
- * Apr 28 1998  Change projection flags to WCS_*
- * May 27 1998	Skip limit checking for linear projection
- * Jun 25 1998	Fix inverse for CAR projection
- * Aug  5 1998	Allan Brighton: Added COE projection (code from A. Wicenec, ESO)
- * Sep 30 1998	Fix bug in COE inverse code to get sign correct
- *
- * Oct 21 1999	Drop unused y from worldpix()
- *
- * Apr  3 2002	Use GLS and SFL interchangeably
- *
- * Feb  3 2004	Let ra be >180 in worldpix() if ref pixel is >180 deg away
- *
- * Jun 20 2006	Initialize uninitialized variables
- *
- * Mar 11 2011	Initialize ZPX
- * Sep  1 2011	Add TPV projection as TAN
- */