subroutine deimos (x, y)

      implicit none
      double precision x, y

*     Plate scale at 550 nm (arcmin/pixel).
      double precision scl
      parameter (scl = 0.1186474d0/60d0)

*     Camera offset from plate centre, arcmin.
      double precision xc, yc
      parameter (xc = 0d0, yc = 0d0)

*     Telescope offset from plate centre, arcmin.
      double precision xt, yt
      parameter (xt = -4.5d0, yt = 0d0)

*     Coefficients of telescope distortion.
      double precision t1, t2, t3
      parameter (t1 = +8.056121663d-2,
     :           t2 = -9.854431080d-1,
     :           t3 = +4.585744900d-2)

*     Coefficients of camera distortion.
      double precision c1, c2, c3, c4, c5
      parameter (c1 = -2.981488078d-2,
     :           c2 = +3.837548331d-2,
     :           c3 = +5.818240867d-2,
     :           c4 = +3.123005326d-3,
     :           c5 = -2.396621712d-4)

      double precision dc, dt, rc, rt, xa, ya

*     Statement functions for the distortions.
      dt(rt) = (t1 + (t2 + t3*rt)*rt)*rt
      dc(rc) = (c1 + (c2 + (c3 + (c4 + c5*rc)*rc)*rc)*rc)*rc

*     Convert (x,y) from pixels to arcmin.
      xa = x * scl
      ya = y * scl

*     Distance from camera centre, in arcmin.
      rc = sqrt((xa-xc)*(xa-xc) + (ya-yc)*(ya-yc))

*     Distance from telescope centre, in arcmin.
      rt = sqrt((xa-xt)*(xa-xt) + (ya-yt)*(ya-yt))

*     Distorted position, in pixels.
      if (rt.gt.0d0) then
         x = x + dt(rt)*(xa-xt)/rt
         y = y + dt(rt)*(ya-yt)/rt
      end if

      if (rc.gt.0d0) then
         x = x + dc(rc)*(xa-xc)/rc
         y = y + dc(rc)*(ya-yc)/rc
      end if

*     Remove first-order displacement.
      x = x + 15.413937422040

      end