procedure axfchk (dcrname, axfname, rname, catname)

string      dcrname, axfname, rname, catname
string      dcrname2
real        z1, z2, zm1, zm2, magcut, xshift
struct      *tmpfile
bool        initspec

begin

   int      wcs, ipx, ipy, imx, imy, ix, iy, imin
   int      il, iu, jl, ju, m, n, icolor, nx, ny, ii, iii
   int      ixl[1000], iyl[1000], itype[1000], ishow[1000]
   int      ixu[1000], iyu[1000], ixc[1000], iyc[1000]
   int      irs, ippp[1000], islit[1000], iobj[1000]
   int      ncat, mcat, icl, iclcat[10000]
   real     rmag[1000], color[1000]
   real     xcat[10000], ycat[10000], magcat[10000]
   real     colorcat[10000]
   real     x, y, x1, y1, xc, yc, dx, dy, dx2, x2, y2
   real     xim, yim, xd, yd, xbeg, ybeg
   real     xsp, ysp, xbegsp, ybegsp
   real     d, dmin, slope, dety, col, tol
   real     mag, cenx, ceny, dright, dleft, x0, y0
   real     p0, p1, p2, dp, w1, w2, ws, disp, rad
   real     dmag, mag2, dmag2
   struct   com, key, line
   string   s1
   string   dcrimage, axffile, rimage, catfile
  
# set some parameters and initialize things
   set clobber = yes
#
# 2048x2048 display and STIS spectroscopic raster
   set stdimage=imt2048
   xd = 2048.
   yd = 2048.   
   xim = 1416.
   yim = 1368.
#   xbeg = (xd-xim)/2.
#   ybeg = (yd-yim)/2.
   xsp = 1456.
   ysp = 1560.
   xbegsp = (xd-xsp)/2.
   ybegsp = (yd-ysp)/2.
   xbeg = (xd-xsp)/2.
   ybeg = (yd-ysp)/2.
#
   w2 = 6300.
   w1 = 4400.
   ws = 5176.
   disp = 5.2
   rad = 2.
   tol = 10.
#
   m = 1000
   n = 0
   mcat = 10000
   ncat = 0
   nx = 0
   ny = 0
#
   catfile = catname
   print(" ")
   print("reading ", catfile)
   list = catfile
   key = "     "
   while (key != "   Ns") {
      if (fscan (list, line) != EOF) {
         key = substr(line, 1, 5)
#         print(line)
       }
    }
#   while (fscan(list, i, x, y, dx, dy, mag, dmag, \
#                      icl) != EOF) {
   while (fscan(list, i, x, y, dx, dy, mag, dmag, \
                      mag2, dmag2, icl, col) != EOF) {
      ncat = ncat + 1
      if (ncat <= mcat) {
         xcat[ncat] = x
         ycat[ncat] = y
         magcat[ncat] = mag
         iclcat[ncat] = icl
#         col = -99.
#          fudge for *cal.cat
#         colorcat[ncat] = col
         colorcat[ncat] = -col
#         catalog tabulated color is r-g but want g-r
#         colorcat[ncat] = 99.
#           fudge for Z5 mask
         if (ncat != i) {
            print("WARNING: catalog read mismatch")
            print(i, x, y, mag)
            print(ncat, xcat[ncat], ycat[ncat], magcat[ncat])
          }
       }
      else {
         print("WARNING: ncat > mcat: ", ncat, mcat)
         ncat = mcat
         print("resetting ncat = ", ncat)
       }
     }
#    (read catfile)
#   print(i, x, y, dx, dy, mag)
#   print(ncat, xcat[ncat], ycat[ncat], magcat[ncat])
   print("read ncat = ", ncat, " objects")
   list = ""
   print("writing objects to axfchk.tmp.cat")
   if (ncat >= 1) {
      print(ncat, tol, > "axfchk.tmp.cat")
      for (i=1; i<=ncat; i+=1) {
          print(xcat[i], ycat[i], magcat[i], iclcat[i], colorcat[i], \
                >> "axfchk.tmp.cat")
       }
    }
#
   rimage = rname
   print(" ")
   print("copying r image ", rimage, " to rtemp")
   imdel("rtemp", go_ahead=yes, verify=no)
   imexpr("repl(1150,1456)", "rtemp", dim="1456,1560")
   imcopy(rimage, "rtemp[1:1416,1:1560]")
   print("displaying copied r image rtemp")
   display(image="rtemp", frame=3,  \
           zscale=no, zrange=no, z1=zm1, z2=zm2, erase = no)
   display(image="rtemp", frame=2,  \
           zscale=no, zrange=no, z1=zm1, z2=zm2, erase = no)
# , > "dev$null")
#
   print(" ")   
   print("making shifted r image: rtemp2")
   imdel("rtemp2", go_ahead=yes, verify=no)
   imcopy("rtemp", "rtemp2")
   imshift("rtemp2", "rtemp2", xshift=1.0, yshift=0.)
#
   print(" ")
   if (initspec == yes) {
      dcrimage = dcrname
      display(image=dcrimage, frame=4, zscale=no, \
              zrange=no, z1=z1, z2=z2, ztrans="log")
      print("displaying dcr image ", dcrimage)
      display(image=dcrimage, frame=1, zscale=no, \
              zrange=no, z1=z1, z2=z2, ztrans="log")
    }
   else {
      dcrimage = dcrname2
      display(image=dcrimage, frame=4, zscale=no, \
              zrange=no, z1=z1, z2=z2, ztrans="log")
      dcrimage = dcrname
      print("displaying dcr image ", dcrimage)
      display(image=dcrimage, frame=1, ztrans="none")
    }
# , > "dev$null")
#   print(" ")   
#   print("making shifted dcr image: temp")
#   imdel("temp", go_ahead=yes, verify=no)
#   imcopy(dcrimage, "temp")
#   imshift("temp", "temp", xshift=1.0, yshift=0.)
#
   axffile = axfname
   print(" ")
   print("reading ", axffile)
   print("1 1 1 : radii ", rad, > "axfchk.tmp")
   print("1 1 1 : txtsize 1", >> "axfchk.tmp")
   print(" ", > "axfchk.tmp2")
   print(" ", > "axfchk.tmp3")
   list = axffile
   if (fscan(list, line) != EOF) print(line)
   if (fscan(list, line) != EOF) print(" ")
   while (fscan(list, i, cenx, ceny, dright, dleft, \
                j, mag, irs, dety, col) != EOF) {
             n = n + 1
             if (n <= m) {
                cenx = cenx + xshift

                x1 = cenx - abs(dright)
                x2 = cenx + abs(dleft)
                x0 = (x1+x2) / 2.
                y0 = ceny - dety
                dx = abs(dright) + abs(dleft)
                dy = 6.
                y1 = y0 - dy/2.
                y2 = y0 + dy/2. 

                p2 = y0 + (w2-ws)/disp 
                p1 = y0 + (w1-ws)/disp 
#                p2 = y0 + (w2-ws)/disp - 200.
#                p1 = y0 + (w1-ws)/disp - 200.
#                  subtract 200 as fudge for Z5 mask
                p0 = (p1+p2)/2. + 1638.-1560. 
                dp = p2-p1
                p1 = p0 - dp/2.
                p2 = p0 + dp/2.
                        
                ixl[n] = nint(x1)
                iyl[n] = nint(p1)
                ixu[n] = nint(x2)
                iyu[n] = nint(p2)
                ixc[n] = nint(cenx-xshift)
                iyc[n] = nint(ceny)

                ippp[n] = j
                rmag[n] = mag
                color[n] = col
                iobj[n] = irs
                ishow[n] = 1

                if (irs > 0) {
                   itype[n] = 1
                   islit[n] = irs

                   print("1 1 1 : color 204 ", >> "axfchk.tmp2")
                   print(x1, p1, " 1 b ", >> "axfchk.tmp2")
                   print(x2, p2, " 1 b ", >> "axfchk.tmp2")
                   print((x1-1), (p1-1), " 1 b ", >> "axfchk.tmp2")
                   print((x2+1), (p2+1), " 1 b ", >> "axfchk.tmp2")

                   print("1 1 1 : color 207 ", >> "axfchk.tmp")
                   print((x1+1), (y1-1), " 1 b ", >> "axfchk.tmp")
                   print((x2-1), (y2+1), " 1 b ", >> "axfchk.tmp")
                   print((x1-4), (y1-10), " 1 : text ", islit[n], \
                         >> "axfchk.tmp")

                   print("1 1 1 : color 204 ", >> "axfchk.tmp")
                   print((cenx-dy/2.), (ceny-dy/2.), " 1 b ",  \
                         >> "axfchk.tmp")
                   print((cenx+dy/2.), (ceny+dy/2.), " 1 b ",  \
                         >> "axfchk.tmp")
                   print((x1-4), (ceny+7), " 1 : text ", iobj[n], \
                         >> "axfchk.tmp")

                   print("1 1 1 : color 207 ", >> "axfchk.tmp3")
                   print((x1+1-xshift), (y1-1), " 1 b ", \
                         >> "axfchk.tmp3")
                   print((x2-1-xshift), (y2+1), " 1 b ", \
                         >> "axfchk.tmp3")
                   print((x1-4-xshift), (y1-10), " 1 : text ", \
                         islit[n], >> "axfchk.tmp3")

                   print("1 1 1 : color 204 ", >> "axfchk.tmp3")
                   print((cenx-dy/2.-xshift), (ceny-dy/2.), " 1 b ",  \
                         >> "axfchk.tmp3")
                   print((cenx+dy/2.-xshift), (ceny+dy/2.), " 1 b ",  \
                         >> "axfchk.tmp3")
                   print((x1-4-xshift), (ceny+7), " 1 : text ", \
                         iobj[n], >> "axfchk.tmp3")

#                   print("1 1 1 : mark circle ", >> "axfchk.tmp")
#                   print(nint(cenx), nint(ceny), " 1 c ", \
#                      >> "axfchk.tmp")

                 }
                else if (irs <= -1000) {
                   itype[n] = 3 
                   islit[n] = abs(irs) - 1000

                   print("1 1 1 : color 205 ", >> "axfchk.tmp")
                   print("1 1 1 : mark plus ", >> "axfchk.tmp")
                   print(nint(cenx), nint(ceny), " 1 + ", \
                      >> "axfchk.tmp")
                   print(x0, (ceny+7), " 1 : text ", iobj[n], \
                         >> "axfchk.tmp")

                   print("1 1 1 : color 205 ", >> "axfchk.tmp3")
                   print("1 1 1 : mark plus ", >> "axfchk.tmp3")
                   print(nint(cenx-xshift), nint(ceny), " 1 + ", \
                      >> "axfchk.tmp3")
                   print((x0-xshift), (ceny+7), " 1 : text ", \
                         iobj[n], >> "axfchk.tmp3")
                 }
                else {
                   itype[n] = 2
                   islit[n] = abs(irs) 

                   print("1 1 1 : color 205 ", >> "axfchk.tmp")
                   print((cenx-dy/2.), (ceny-dy/2.), " 1 b ",  \
                         >> "axfchk.tmp")
                   print((cenx+dy/2.), (ceny+dy/2.), " 1 b ",  \
                         >> "axfchk.tmp")
                   print(x0, (ceny+7), " 1 : text ", iobj[n], \
                         >> "axfchk.tmp")

                   print("1 1 1 : color 205 ", >> "axfchk.tmp3")
                   print((cenx-dy/2.-xshift), (ceny-dy/2.), " 1 b ",  \
                         >> "axfchk.tmp3")
                   print((cenx+dy/2.-xshift), (ceny+dy/2.), " 1 b ",  \
                         >> "axfchk.tmp3")
                   print((x0-xshift), (ceny+7), " 1 : text ", \
                         iobj[n], >> "axfchk.tmp3")
                 }
              }
             else {
               print("WARNING: n > m: ", n, m)
               n = m
               print("resetting n = ", n)
              }
      }
#      (read axffile)
     list = ""
#
   print("displaying n = ", n, " objects ")
   tvmark(2, "", command="axfchk.tmp3", > "dev$null")     
##   if (initspec == yes)
##      tvmark(1, "", command="axfchk.tmp2", > "dev$null")     
   tvmark(1, "", command="axfchk.tmp", > "dev$null")     
#     print("displaying tempbox")
#     display("tempbox", 3, ztrans="none")
#
#     print("making and displaying tempsum")
#     imarith("tempdcr", "+", "tempbox", "tempsum"
#     display("tempsum", 4, ztrans="none")
#
   while (fscan (imcur, x, y, wcs, com) != EOF) {
       key = substr(com, 1, 1)

       if (key == "q")
          break
#       else if (key == "r") {
#          print(" ")
#          print("redisplaying ", dcrimage)
#          display(image=dcrimage, frame=1, zscale=no, zrange=no, 
#                  z1=z1, z2=z2)
#         }
       else if (key == "p") {
          print(" ")
#          print(x, y, wcs, com)
          if (fscan (imcur, x1, y1, wcs, com) != EOF) {
#             print(x1, y1, wcs, com)
            imx = nint(x)
            ipx = nint(x1)
            imy = nint(y)
            ipy = nint(y1)
            if (x < x1 && y < y1) {
               print("1 1 1 : l ", imy, ipy, > "axfchk.tmp")
               print("1 1 1 : x ", imx, ipx, >> "axfchk.tmp")
               print("implotting average of y = ", imy, " to ", ipy)
               print("  plot x range = ", imx, " to ", ipx)
               if (initspec == yes) 
                  implot(dcrimage, 1000, coords="axfchk.tmp")
               else
                  implot(dcrname2, 1000, coords="axfchk.tmp")
              }
            else {
               beep()
               print("improper x or y range: ", imx, ipx, imy, ipy)
              }
           }
         }
#         (key == "p")
       else if (key == "x" || key == "y") {
         n = n + 1
         if (n <= m) {
            print(" ")
            print("adding box ", n)
            print(x, y, wcs, com)
            ixl[n] = nint(x)
            iyl[n] = nint(y)
            print(ixl[n], iyl[n])

            if (fscan (imcur, x, y, wcs, com) != EOF) {
               print(x, y, wcs, com)
               ixu[n] = nint(x)
               iyu[n] = nint(y)
               print(ixu[n], iyu[n])
               ixc[n] = (ixl[n] + ixu[n]) / 2
               iyc[n] = (iyl[n] + iyu[n]) / 2
               if (key == "x") {
                  itype[n] = 1
                  print("1 1 1 : color 205 ", > "axfchk.tmp")
                  nx = nx + 1
                }
               else {
                  itype[n] = 2
                  print("1 1 1 : color 206 ", > "axfchk.tmp")
                  ny = ny + 1
                }
               print(ixl[n], iyl[n], wcs, " b ", >> "axfchk.tmp")
               print(ixu[n], iyu[n], wcs, " b ", >> "axfchk.tmp")
               tvmark(1, "", command="axfchk.tmp", > "dev$null") 
             }
         else {
               print("WARNING: n > m: ", n, m)
               n = m
               print("resetting n = ", n)
           }
          }
#          (n <= m)
        }
#        (key == "x" || "y")
       else if (key == "i" || key == "t") {
         imin = 0
         if (wcs == 201) 
            x = x + xshift
         ix = nint(x)
         iy = nint(y)
         print(" ")
         for (i=1; i<=n; i+=1) {
             if (ix >= ixl[i] && ix <= ixu[i] &&
                 iy >= iyl[i] && iy <= iyu[i]) {
                printf("slit %3d  obj %5d  mag %6.2f  color %6.2f"// \
                       "  ppp # %5d  x %4d  y %5d \n", \
                       islit[i], iobj[i], rmag[i], color[i], ippp[i], \
                       nint(ixc[i]), nint(iyc[i]))
                imin = i
              }
          }
         if (imin == 0) 
            print("x, y not in any spectrum")
         else if (key == "t" && wcs == 101) {
            print(" ", > "axfchk.tmp2")
            i = imin
            if (ishow[i] == 0) {
               print("1 1 1 : color 204 ", >> "axfchk.tmp2")
               ishow[i] = 1
             }
            else {
               print("1 1 1 : color 203 ", >> "axfchk.tmp2")
               ishow[i] = 0
             }            
            print(ixl[i], iyl[i], " 1 b ", >> "axfchk.tmp2")
            print(ixu[i], iyu[i], " 1 b ", >> "axfchk.tmp2")
            print((ixl[i]-1), (iyl[i]-1), " 1 b ", >> "axfchk.tmp2")
            print((ixu[i]+1), (iyu[i]+1), " 1 b ", >> "axfchk.tmp2")
            tvmark(1, "", command="axfchk.tmp2", > "dev$null") 
          }
#         (key == "t")
        }
#        (key == "i" || "t")
       else if (key == "m") {
         if (wcs != 201) {
            beep()
            print(" ")
            print("not in r image frame 2 or 3")
          }
         else {
            print(" ")
#            print("x = ", x, " y = ", y)
            dmin = 1.e6
            imin = 0
            for (i=1; i<=ncat; i+=1) {
                d = sqrt( (x-xcat[i])**2. +  \
                          (y-ycat[i])**2. )
                if (d <= dmin) {
                   dmin = d
                   imin = i
                 }
              }
            if (imin != 0 && dmin <= tol) {
                printf("ppp # %5d  cl %2d  mag %6.2f  color %6.2f"// \
                       "  x %7.2f  y %7.2f \n", \
                       imin, iclcat[imin], magcat[imin], \
                        colorcat[imin], xcat[imin], ycat[imin])
             }
            else {
               print("no object within ", tol, \
                     " pixels of x, y = ",  x, y)
             }
          }
        }
#        (key == "m")
       else if (key == "n" || key == "d" || key == "o") {
         if (wcs != 201 && wcs != 301) {
            beep()
            print(" ")
            print("not in r image frame 2 or frame 3")
          }
         else 
            print(x, y, > "axfchk.tmp")
            readcat
            list = "axfchk.tmp"
            if (fscan(list, i) != EOF) {
               if (i != 0) {
                  if (fscan(list, j, icl, mag, col, x, y) != EOF) {
                   if (key == "n") {
                     print(" ", > "axfchk.tmp3")
                     if (icl == 1 || icl == 2) {
                        print(" 1 1 1 : color 206", >> "axfchk.tmp3") 
                        if (mag <= magcut) {
                           print((x-2), (y-2), " 1 b ", \
                                 >> "axfchk.tmp3")   
                           print((x+2), (y+2), " 1 b ", \
                                 >> "axfchk.tmp3")   
                         }
                        else {
                           print(" 1 1 1 : pointsize 3",
                                 >> "axfchk.tmp3")   
                           print(" 1 1 1 : mark point", \
                                 >> "axfchk.tmp3")   
                           print(x, y, " 1 . ", >> "axfchk.tmp3")   
                         }
                      }
#                      (if galaxy)
                     else if (icl == 3 || icl == 4) {
                        print(" 1 1 1 : color 210", >> "axfchk.tmp3") 
                        if (mag <= magcut) {
                           print(" 1 1 1 : txtsize 1",
                                 >> "axfchk.tmp3")   
                           print(" 1 1 1 : mark plus", \
                                 >> "axfchk.tmp3")   
                           print(x, y, " 1 + ", >> "axfchk.tmp3")   
                         }
                        else {
                           print(" 1 1 1 : pointsize 3",
                                 >> "axfchk.tmp3")   
                           print(" 1 1 1 : mark point", \
                                 >> "axfchk.tmp3")   
                           print(x, y, " 1 . ", >> "axfchk.tmp3")   
                         }
                      }
#                      (if star)
                     if (wcs == 201) 
                        tvmark(2, "", command="axfchk.tmp3", \
                               > "dev$null")             
                     else if (wcs == 301)
                        tvmark(3, "", command="axfchk.tmp3", \
                               > "dev$null")             
                    }
#                    (key == "n")
                   else if (key == "d") {
                     if (icl == 1 || icl == 2) 
                        imx = nint(x) - 1
                     else if (icl == 3 || icl == 4) 
                        imx = nint(x) - 2
                      
                     ipx = nint(x) + 3
                     imy = nint(y) - 2
                     ipy = nint(y) + 2
                     s1 = "rtemp2[" + imx // ":" + ipx // "," \
                                    + imy // ":" + ipy // "]" 
#                     print("s1 = ", s1)
                      xc = (xbeg + (real(imx) + real(ipx) ) / 2. - 1.) / xd
#                      note one pixel shift in x only, and use of shifted
#                       image "rtemp2" to account for funny iraf display bug 
                      yc = (ybeg + (real(imy) + real(ipy) ) / 2. ) / yd
#                     print("xc, yc = ", xc, yc)
                      if (wcs == 201) 
                       display(image=s1, frame=2, xcenter=xc, ycenter=yc, \
                               zscale=no, zrange=no, z1=zm1, z2=zm2, 
                               erase = no, > "dev$null")
                      else if (wcs == 301)
                       display(image=s1, frame=3, xcenter=xc, ycenter=yc, \
                               zscale=no, zrange=no, z1=zm1, z2=zm2, 
                               erase = no, > "dev$null")
                    }
#                    (key == "d")              
                   }
#                   (read "axfchk.tmp")
                }
             }
            list = ""
        }
#        (key == "n" || "d")
       else if (key == "s") {
         print(" ")
         print("enter slit number:")
         if (scan(j) != EOF) {
            ix = 0
            for (i=1; i<=n; i+=1) {
             if (islit[i] == j) {
                printf("slit %3d  obj %5d  mag %6.2f  color %6.2f"// \
                       "  ppp # %5d  x %4d  y %5d \n", \
                       islit[i], iobj[i], rmag[i], color[i], ippp[i], \
                       nint(ixc[i]), nint(iyc[i]))
                ix = 1
              }
             }
            if (ix == 0) 
               print("not found")
          }
        }
#        (key == "s")
       else if (key == "l") {
         print(" ")
         print(dcrimage, > "axfchk.tmp.list")
         for (i=1; i<=n; i+=1) {
            if (ishow[i] == 1) {
               printf("slit %3d  obj %5d  mag %6.2f  color %6.2f"// \
                      "  ppp # %5d  x %4d  y %5d \n", \
                      islit[i], iobj[i], rmag[i], color[i], ippp[i], \
                      nint(ixc[i]), nint(iyc[i]))
               printf("slit %3d  obj %5d  mag %6.2f  color %6.2f"// \
                      "  ppp # %5d  x %4d  y %5d \n", \
                      islit[i], iobj[i], rmag[i], color[i], ippp[i], \
                      nint(ixc[i]), nint(iyc[i]), >> "axfchk.tmp.list")
             }
          }
        }
#        (key == "l")
       else {
         beep()
         print(" ")
         print("WARNING: unknown command ", key)
        }
   }
#   (read image cursor)
   set clobber = no
end