diff options
Diffstat (limited to 'pkg/tbtables/fitsio/ftwldp.f')
| -rw-r--r-- | pkg/tbtables/fitsio/ftwldp.f | 289 |
1 files changed, 289 insertions, 0 deletions
diff --git a/pkg/tbtables/fitsio/ftwldp.f b/pkg/tbtables/fitsio/ftwldp.f new file mode 100644 index 00000000..69f78137 --- /dev/null +++ b/pkg/tbtables/fitsio/ftwldp.f @@ -0,0 +1,289 @@ +C------------------------------------------------------------------------------ + subroutine ftwldp(xpix,ypix,xref,yref,xrefpix,yrefpix, + & xinc,yinc,rot,type,xpos,ypos,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 position from pixel coordinates +C does: -SIN, -TAN, -ARC, -NCP, -GLS, -MER, -AIT projections +C returns 0 = good, +C 501 = angle too large for projection; +C Input: +C dbl xpix x pixel number (RA or long without rotation) +C dbl ypiy y pixel number (dec or lat without rotation) +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 xpos x (RA) coordinate (deg) +C dbl ypos y (dec) coordinate (deg) +C int status error status flag, zero + + integer status + double precision xpix,ypix,xref,yref,xrefpix,yrefpix + double precision xinc,yinc,rot,xpos,ypos + character*(*) type + integer error1,error4 + parameter (error1=501) + parameter (error4=504) + + double precision cosr,sinr,dx,dy,dz,temp + double precision sins,coss,dect,rat,dt,l,m,mg,da,dd,cos0,sin0 + double precision dec0,ra0,decout,raout + double precision geo1,geo2,geo3 + 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 *** Offset from ref pixel + dx = (xpix-xrefpix) * xinc + dy = (ypix-yrefpix) * yinc +C *** Take out rotation + cosr = dcos(rot*cond2r) + sinr = dsin(rot*cond2r) + if (rot .ne. 0.0) then + temp = dx * cosr - dy * sinr + dy = dy * cosr + dx * sinr + dx = temp + end if +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 *** default, linear result for error return + xpos = xref + dx + ypos = yref + dy +C *** Convert to radians + ra0 = xref * cond2r + dec0 = yref * cond2r + l = dx * cond2r + m = dy * cond2r + sins = l*l + m*m + decout = 0.0 + raout = 0.0 + cos0 = dcos(dec0) + sin0 = dsin(dec0) +C *** Process by case + if (itype .eq. 0) then +C *** LINEAR + rat = ra0 + l + dect = dec0 + m + else if (itype .eq. 1) then +C *** SINE from '-SIN' type + if (sins .gt. 1.0) then + status = error1 + goto 30 + end if + coss = dsqrt(1.0 - sins) + dt = sin0 * coss + cos0 * m + if ((dt .gt. 1.0) .or. (dt .lt. -1.0)) then + status = error1 + goto 30 + end if + dect = dasin(dt) + rat = cos0 * coss - sin0 * m + if ((rat .eq. 0.0) .and. (l .eq. 0.0)) then + status = error1 + goto 30 + end if + rat = datan2 (l, rat) + ra0 + else if (itype .eq. 2) then +C *** TANGENT from '-TAN' type + if (sins .gt. 1.0) then + status = error1 + goto 30 + end if + dect = cos0 - m * sin0 + if (dect .eq. 0.0) then + status = error1 + goto 30 + end if + rat = ra0 + datan2(l, dect) + dect = datan(dcos(rat-ra0) * (m * cos0 + sin0) / dect) + else if (itype .eq. 3) then +C *** Arc from '-ARC' type + if (sins .ge. twopi * twopi / 4.0) then + status = error1 + goto 30 + end if + sins = dsqrt(sins) + coss = dcos(sins) + if (sins .ne. 0.0) then + sins = dsin(sins) / sins + else + sins = 1.0 + end if + dt = m * cos0 * sins + sin0 * coss + if ((dt .gt. 1.0) .or. (dt .lt. -1.0)) then + status = error1 + goto 30 + end if + dect = dasin(dt) + da = coss - dt * sin0 + dt = l * sins * cos0 + if ((da .eq. 0.0) .and. (dt .eq. 0.0)) then + status = error1 + goto 30 + end if + rat = ra0 + datan2(dt, da) + else if (itype .eq. 4) then +C *** North Celestial Pole from '-NCP' type + dect = cos0 - m * sin0 + if (dect .eq. 0.0) then + status = error1 + goto 30 + end if + rat = ra0 + datan2(l, dect) + dt = dcos(rat-ra0) + if (dt .eq. 0.0) then + status = error1 + goto 30 + end if + dect = dect / dt + if ((dect .gt. 1.0) .or. (dect .lt. -1.0)) then + status = error1 + goto 30 + end if + dect = dacos(dect) + if (dec0 .lt. 0.0) dect = -dect + else if (itype .eq. 5) then +C *** Global Sinusoid from '-GLS' type + dect = dec0 + m + if (dabs(dect) .gt. twopi/4.0) then + status = error1 + goto 30 + end if + coss = dcos(dect) + if (dabs(l) .gt. twopi*coss/2.0) then + status = error1 + goto 30 + end if + rat = ra0 + if (coss .gt. deps) rat = rat + l / 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 = dcos(yref * cond2r) + if (geo1 .le. 0.0) geo1 = 1.0 + rat = l / geo1 + ra0 + if (dabs(rat - ra0) .gt. twopi) then + status = error1 + goto 30 + end if + dt = 0.0 + if (geo2 .ne. 0.0) dt = (m + geo3) / geo2 + dt = dexp(dt) + dect = 2.0 * datan(dt) - twopi / 4.0 + else if (itype .eq. 7) then +C *** Aitoff from '-AIT' type + 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) + rat = ra0 + dect = dec0 + if ((l .eq. 0.0) .and. (m .eq. 0.0)) goto 20 + dz = 4.0-l*l/(4.0*geo1*geo1)-((m+geo3)/geo2)*((m+geo3)/geo2) + if ((dz .gt. 4.0) .or. (dz .lt. 2.0)) then + status = error1 + goto 30 + end if + dz = 0.5 * dsqrt(dz) + dd = (m+geo3) * dz / geo2 + if (dabs(dd) .gt. 1.0) then + status = error1 + goto 30 + end if + dd = dasin(dd) + if (dabs(dcos(dd)) .lt. deps) then + status = error1 + goto 30 + end if + da = l * dz / (2.0 * geo1 * dcos(dd)) + if (dabs(da) .gt. 1.0) then + status = error1 + goto 30 + end if + da = dasin(da) + rat = ra0 + 2.0 * da + dect = dd + else if (itype .eq. 8) then +C *** Stereographic from '-STG' type + dz = (4.0 - sins) / (4.0 + sins) + if (dabs(dz) .gt. 1.0) then + status = error1 + goto 30 + end if + dect = dz * sin0 + m * cos0 * (1.0+dz) / 2.0 + if (dabs(dect) .gt. 1.0) then + status = error1 + goto 30 + end if + dect = dasin(dect) + rat = dcos(dect) + if (dabs(rat) .lt. deps) then + status = error1 + goto 30 + end if + rat = l * (1.0+dz) / (2.0 * rat) + if (dabs(rat) .gt. 1.0) then + status = error1 + goto 30 + end if + rat = dasin(rat) + mg = 1.0 + dsin(dect)*sin0 + dcos(dect)*cos0*dcos(rat) + if (dabs(mg) .lt. deps) then + status = error1 + goto 30 + end if + mg = 2.0 * (dsin(dect)*cos0 - dcos(dect)*sin0*dcos(rat)) / mg + if (dabs(mg-m) .gt. deps) rat = twopi/2.0 - rat + rat = ra0 + rat + else +C *** Unsupported Projection + status = error4 + goto 30 + end if + 20 continue +C *** Return RA in range + raout = rat + decout = dect + if (raout-ra0 .gt. twopi/2.0) raout = raout - twopi + if (raout-ra0 .lt. -twopi/2.0) raout = raout + twopi + if (raout .lt. 0.0) raout = raout + twopi +C *** Correct units back to degrees + xpos = raout / cond2r + ypos = decout / cond2r + 30 continue + end |