diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2016-10-25 20:57:49 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2016-10-25 20:57:49 (GMT) |
commit | d1c4bf158203c4e8ec29fdeb83fd311e36320885 (patch) | |
tree | 15874534e282f67505ce4af5ba805a1ff70ec43e /funtools/wcs/worldpos.c | |
parent | e19a18e035dc4d0e8e215f9b452bb9ef6f58b9d7 (diff) | |
parent | 339420dd5dd874c41f6bab5808291fb4036dd022 (diff) | |
download | blt-d1c4bf158203c4e8ec29fdeb83fd311e36320885.zip blt-d1c4bf158203c4e8ec29fdeb83fd311e36320885.tar.gz blt-d1c4bf158203c4e8ec29fdeb83fd311e36320885.tar.bz2 |
Merge commit '339420dd5dd874c41f6bab5808291fb4036dd022' as 'funtools'
Diffstat (limited to 'funtools/wcs/worldpos.c')
-rw-r--r-- | funtools/wcs/worldpos.c | 693 |
1 files changed, 693 insertions, 0 deletions
diff --git a/funtools/wcs/worldpos.c b/funtools/wcs/worldpos.c new file mode 100644 index 0000000..8693789 --- /dev/null +++ b/funtools/wcs/worldpos.c @@ -0,0 +1,693 @@ +/* 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 + */ |