Repatch (from linux) of OSX IRAF

1 files changed, 230 insertions, 0 deletions

diff --git a/pkg/tbtables/fitsio/ftxypx.f b/pkg/tbtables/fitsio/ftxypx.f
new file mode 100644
index 00000000..4a21e55f
--- /dev/null
+++ b/pkg/tbtables/fitsio/ftxypx.f
@@ -0,0 +1,230 @@
+C------------------------------------------------------------------------------
+        subroutine ftxypx(xpos,ypos,xref,yref,xrefpix,yrefpix,
+     &                    xinc,yinc,rot,type,xpix,ypix,status)
+
+C       Fortran version of worldpos.c -- WCS Algorithms from Classic AIPS
+C       Translated by James Kent Blackburn -- HEASARC/GSFC/NASA -- November 1994
+C       routine to determine accurate pixel coordinates from an RA and Dec
+C       does: -SIN, -TAN, -ARC, -NCP, -GLS, -MER, -AIT projections
+C       returns   0 = good, 
+C               501 = angle too large for projection;
+C               502 = bad values
+C               503 = ???undocumented error - looks like an underflow???
+C Input:
+C dbl   xpos   x (RA) coordinate (deg)
+C dbl   ypos   y (dec) coordinate (deg)
+C dbl   xref    x reference coordinate value (deg)
+C dbl   yref    y reference coordinate value (deg)
+C dbl   xrefpix x reference pixel
+C dbl   yrefpix y reference pixel
+C dbl   xinc    x coordinate increment (deg)
+C dbl   yinc    y coordinate increment (deg)
+C dbl   rot     rotation (deg)  (from N through E)
+C chr   type    projection type code e.g. "-SIN"
+C Output:
+C dbl   xpix    x pixel number  (RA or long without rotation)
+C dbl   ypiy    y pixel number  (dec or lat without rotation)
+C int   status  error status flag, zero
+
+        integer           status
+        double precision  xpos,ypos,xref,yref,xrefpix,yrefpix
+        double precision  xinc,yinc,rot,xpix,ypix
+        character*(*)     type
+        integer           error1,error2,error3,error4
+        parameter         (error1=501)
+        parameter         (error2=502)
+        parameter         (error3=503)
+        parameter         (error4=504)
+        double precision  dx,dy,dz,r,ra0,dec0,ra,dec
+        double precision  coss,sins,dt,da,dd,sint,oldxpos
+        double precision  l,m,geo1,geo2,geo3,sinr,cosr
+        double precision  cond2r
+        parameter         (cond2r=1.745329252d-2)
+        double precision  twopi,deps
+        parameter         (twopi = 6.28318530717959)
+        parameter         (deps = 1.0d-5)
+        integer           i,itype
+        character*4       ctypes(8)
+        data ctypes/ '-SIN', '-TAN', '-ARC', '-NCP',
+     &               '-GLS', '-MER', '-AIT', '-STG' /
+
+        if (status .gt. 0) return
+C   *** 0 hour wrap around test
+        oldxpos = xpos
+        dt = (xpos - xref)
+        if (dt .gt. +180) xpos = xpos - 360
+        if (dt .lt. -180) xpos = xpos + 360
+C   *** Default values - Linear
+        dx = xpos - xref
+        dy = ypos - yref
+        dz = 0.0
+C   *** Correct for rotation
+        r = rot * cond2r
+        cosr = dcos(r)
+        sinr = dsin(r)
+        dz = dx * cosr + dy * sinr
+        dy = dy * cosr - dx * sinr
+        dx = dz
+C   *** Check axis increments - bail out if either 0
+        if ((xinc .eq. 0.0) .or. (yinc .eq. 0.0)) then
+          xpix = 0.0
+          ypix = 0.0
+          status = error2
+          goto 30
+        end if
+        xpix = dx / xinc + xrefpix
+        ypix = dy / yinc + yrefpix
+C   *** Find type of coordinate transformation (0 is linear)
+        itype = 0
+        do 10 i = 1, 8
+          if (ctypes(i) .eq. type) itype = i
+  10    continue
+C   *** Done if linear
+        if (itype .eq. 0) goto 30
+C   *** Non-Linear position
+        ra0 = xref * cond2r
+        dec0 = yref * cond2r
+        ra = xpos * cond2r
+        dec = ypos * cond2r
+C   *** Compute directional cosine
+        coss = dcos(dec)
+        sins = dsin(dec)
+        l = dsin(ra-ra0) * coss
+        sint = sins * dsin(dec0) + coss * dcos(dec0) * dcos(ra-ra0)
+C   *** Process by case
+        if (itype .eq. 1) then
+C   *** SINE from '-SIN' type
+          if (sint .lt. 0.0) then 
+            status = error1
+            goto 30
+          end if
+          m = sins * dcos(dec0) - coss * dsin(dec0) * dcos(ra-ra0)
+        else if (itype .eq. 2) then
+C   *** TANGENT from '-TAN' type
+          if (sint .le. 0.0) then 
+            status = error1
+            goto 30
+          end if
+          m = sins * dsin(dec0) + coss * dcos(dec0) * dcos(ra-ra0)
+          l = l / m
+          m = (sins*dcos(dec0) - coss*dsin(dec0)*dcos(ra-ra0)) / m
+        else if (itype .eq. 3) then
+C   *** Arc from '-ARC' type
+          m = sins*dsin(dec0) + coss*dcos(dec0)*dcos(ra-ra0)
+          if (m .lt. -1.0) m = -1.0
+          if (m .gt. 1.0) m = 1.0
+          m = dacos(m)
+          if (m .ne. 0) then
+            m = m / dsin(m)
+          else
+            m = 1.0
+          end if
+          l = l * m
+          m = (sins*dcos(dec0) - coss*dsin(dec0)*dcos(ra-ra0)) * m
+        else if (itype .eq. 4) then
+C   *** North Celestial Pole from '-NCP' type
+          if (dec0 .eq. 0.0) then
+            status = error1
+            goto 30
+          else
+            m = (dcos(dec0) - coss * dcos(ra-ra0)) / dsin(dec0)
+          end if
+        else if (itype .eq. 5) then
+C   *** Global Sinusoid from '-GLS' type
+          dt = ra - ra0
+          if (dabs(dec) .gt. twopi/4.0) then 
+            status = error1
+            goto 30
+          end if
+          if (dabs(dec0) .gt. twopi/4.0) then 
+            status = error1
+            goto 30
+          end if
+          m = dec - dec0
+          l = dt * coss
+        else if (itype .eq. 6) then
+C   *** Mercator from '-MER' type
+          dt = yinc * cosr + xinc * sinr
+          if (dt .eq. 0.0) dt = 1.0
+          dy = (yref/2.0 + 45.0) * cond2r
+          dx = dy + dt / 2.0 * cond2r
+          dy = dlog(dtan(dy))
+          dx = dlog(dtan (dx))
+          geo2 = dt * cond2r / (dx - dy)
+          geo3 = geo2 * dy
+          geo1 = cos (yref * cond2r)
+          if (geo1 .le. 0.0) geo1 = 1.0
+          dt = ra - ra0
+          l = geo1 * dt
+          dt = dec / 2.0 + twopi / 8.0
+          dt = dtan(dt)
+          if (dt .lt. deps) then 
+            status = error2
+            goto 30
+          end if
+          m = geo2 * dlog(dt) - geo3
+        else if (itype .eq. 7) then
+C   *** Aitoff from '-AIT' type
+          l = 0.0
+          m = 0.0
+          da = (ra - ra0) / 2.0
+          if (dabs(da) .gt. twopi/4.0) then
+            status = error1
+            goto 30
+          end if
+          dt = yinc * cosr + xinc * sinr
+          if (dt .eq. 0.0) dt = 1.0
+          dt = dt * cond2r
+          dy = yref * cond2r
+          dx = dsin(dy+dt)/dsqrt((1.0+dcos(dy+dt))/2.0) -
+     &         dsin(dy)/dsqrt((1.0+dcos(dy))/2.0)
+          if (dx .eq. 0.0) dx = 1.0
+          geo2 = dt / dx
+          dt = xinc * cosr - yinc * sinr
+          if (dt .eq. 0.0) dt = 1.0
+          dt = dt * cond2r
+          dx = 2.0 * dcos(dy) * dsin(dt/2.0)
+          if (dx .eq. 0.0) dx = 1.0
+          geo1 = dt*dsqrt((1.0+dcos(dy)*dcos(dt/2.0))/2.0)/dx
+          geo3 = geo2 * dsin(dy) / dsqrt((1.0+dcos(dy))/2.0)
+          dt = dsqrt ((1.0 + dcos(dec) * dcos(da))/2.0)
+          if (dabs(dt) .lt. deps) then
+            status = error3
+            goto 30
+          end if
+          l = 2.0 * geo1 * dcos(dec) * dsin(da) / dt
+          m = geo2 * dsin(dec) / dt - geo3
+        else if (itype .eq. 8) then
+C   *** Stereographic from '-STG' type
+          da = ra - ra0
+          if (dabs(dec) .gt. twopi/4.0) then
+            status = error1
+            goto 30
+          end if
+          dd = 1.0 + sins*dsin(dec0) + coss*dcos(dec0)*dcos(da)
+          if (dabs(dd) .lt. deps) then 
+            status = error1
+            goto 30
+          end if
+          dd = 2.0 / dd
+          l = l * dd
+          m = dd * (sins*dcos(dec0) - coss*dsin(dec0)*dcos(da))
+        else
+C   *** Unsupported Projection
+          status = error4
+          goto 30
+        end if
+C   *** Convert back to degrees
+        dx = l / cond2r
+        dy = m / cond2r
+C   *** Correct for rotation
+        dz = dx * cosr + dy * sinr
+        dy = dy * cosr - dx * sinr
+        dx = dz
+C   *** Convert to PIXELS ... yeah!
+        xpix = dx / xinc + xrefpix
+        ypix = dy / yinc + yrefpix
+  30    continue
+C   *** reset xpos to correct for in place modification
+        xpos = oldxpos
+        end