procedure AptReduce ()

# task AptReduce = home$apt/AptReduce.cl

# An IRAF program to process APT images.

# This works for 800x1200 pixel images, taken with "gain lo speed lo".
# The images must have names like im?????.fits.  
#
# Usage:
#
# Put all the images from a night into a directory somewhere. cd to the 
# directory and run AptReduce. The images will be sorted into six
# categories, using the following filenames:
#
#     Images     - a list of all the im?????.fits images
#
#     Flats      - 800x1200 images that appear to be flat fields on the
#                  basis of their statistics:
#                  (mean > 16000.0 && mean < 35000.0 && stddev < 400)
#
#     Darks      - 800x1200 images that appear to be dark frames on the 
#                  basis of their statistics:
#                  (mean > 1900.0 && mean < 2500.0 && stddev < 10)
#
#     Twilights  - 800x1200 images that have IMAGETYP = 'twilight' and that
#                  satisfy:
#                  (mean > 6000.0 && mean < 35000.0)
#
#     Stars      - all other 800x1200 images
#     Unknowns   - all other images
#
# If files of the above names already exist, their contents are assumed to
# be correct, and they are not regenerated. So, if you want to redo the
# flat fields after removing some dubious frames, simply edit the Flats
# file, delete Flat.fits, and re-run AptReduce.
#
# The program begins by creating a combined flat image (called Flat.fits) from
# the Flats, and uses it to mask out stars in the twilight images. The
# twilight images are then combined to produce a master twilight image
# called Twilight. At this point the program stops and asks you to
# examine the Twilight and Flat images. You should then create an image
# called SuperFlat (perhaps by "imcopy Twilight SuperFlat", or by combining
# images from a previous night), and rerun AptReduce to process object
# files (in Stars).
#
# To process individual images, use AptProc, for example:
#
#    AptProc (dir="./", files=im23232, useflat+, bias=0, 
#      short-, flat=SuperFlat)
#
# The reduced output files have a "p" (for "processed") prepended to 
# their filenames.
#
# To start from scratch:
#
# del Darks,Flats,Images,Stars,Twilights,Unknowns
# imdel Flat,Twilight,SuperFlat
#
#
# Michael Ashley / UNSW / 27-Sep-1998

char   version = "14Jan1999"  {prompt="Program version"}

struct *imglist

begin
  real   mean, stddev
  string ImageFile, StarFile, FlatFile, TwilightFile
  string DarkFile, UnknownFile, Flat, SuperFlat
  string Twilight, Dark, img, imagetyp, t1, t2, t3
  int    n1, n2, i, NumTwilights
  real   pmin, pmax, percent

  bool   HaveImages
  bool   HaveStars
  bool   HaveFlats
  bool   HaveTwilights
  bool   HaveDark
  bool   HaveDarks
  bool   HaveUnknowns
  bool   HaveFlat
  bool   HaveSuperFlat
  bool   HaveTwilight

# Restore defaults for all the IRAF routines we use.

  unlearn AptProc
  unlearn sections
  unlearn imcopy
  unlearn imcombine
  unlearn imdelete 
  unlearn imstatistics 
  unlearn imarith
  unlearn imsurfit

  ImageFile     = "Images"
  StarFile      = "Stars"
  FlatFile      = "Flats"
  TwilightFile  = "Twilights"
  DarkFile      = "Darks"
  UnknownFile   = "Unknowns"
  Flat          = "Flat"
  SuperFlat     = "SuperFlat"
  Twilight      = "Twilight"
  Dark          = "Dark"

  HaveImages    = access(ImageFile)
  HaveStars     = access(StarFile)
  HaveFlats     = access(FlatFile)
  HaveTwilights = access(TwilightFile)
  HaveDarks     = access(DarkFile)
  HaveUnknowns  = access(UnknownFile)

  HaveFlat      = access(Flat//".imh")      || access(Flat//".fits")
  HaveSuperFlat = access(SuperFlat//".imh") || access(SuperFlat//".fits")
  HaveTwilight  = access(Twilight//".imh")  || access(Twilight//".fits")
  HaveDark      = access(Dark//".imh")      || access(Dark//".fits")

  if (!HaveImages) {
    printf ("finding images...")
    sections ("im?????.fits", option="fullname", > ImageFile)
  }

  if (!HaveStars) {
     printf ("", > StarFile)
  }

  if (!HaveFlats) {
     printf ("", > FlatFile)
  }

  if (!HaveTwilights) {
     printf ("", > TwilightFile)
  }

  if (!HaveDarks) {
     printf ("", > DarkFile)
  }

  if (!HaveUnknowns) {
     printf ("", > UnknownFile)
  }

  imglist = ImageFile
  if (!HaveStars || !HaveFlats || !HaveTwilights || !HaveDarks || !HaveUnknowns) {
    printf ("identifying images...")
    while (fscan (imglist, img) != EOF) {
      i = strlen (img)
      if (substr (img, i-4, i) == ".fits") {
        img = substr (img, 1, i-5)
      }
      imgets (image = img, param = "i_naxis1")
      n1 = int (imgets.value)
      imgets (image = img, param = "i_naxis2")
      n2 = int (imgets.value)
      if (n1 == 800 && n2 == 1200) {
        imgets (image = img, param = "IMAGETYP")
        imagetyp = substr (imgets.value, 1, 8)
        imstatistics (images=img//"[400:500,400:500]",
                      fields="mean,stddev",format-) | scan (mean, stddev)
        if (imagetyp == "twilight" && mean > 6000.0 && mean < 35000.0) {
          if (!HaveTwilights) {
            printf ("%s\n", img, >> TwilightFile)
          }
        } else if (!HaveStars || !HaveFlats || !HaveDarks) {
          if (mean > 16000.0 && mean < 35000.0 && stddev < 400) {
            if (!HaveFlats) {
              printf ("%s\n", img, >> FlatFile)
            }
          } else if (mean > 1900.0 && mean < 2500.0 && stddev < 10) {
            if (!HaveDarks) {
              printf ("%s\n", img, >> DarkFile)
            }
          } else {
            if (!HaveStars) {
              printf ("%s\n", img, >> StarFile)
            }
          }
        }
      } else {
        if (!HaveUnknowns) {
          printf ("%s\n", img, >> UnknownFile)
        }
      }
    }
  }

  if (!HaveFlat) {

# Create some temporary files.
  
    t1 = mktemp ("tmp")
    t2 = mktemp ("tmp")
    t3 = mktemp ("tmp")
    printf ("processing flat field images...\n")

# Make sure there are at least two flat field images:

    imglist = FlatFile
    i = 0
    while (fscan (imglist, img) != EOF) {
      i += 1
    }
    if (i < 2) {
      printf ("The file \"Flats\" contains fewer than two images, this is\n")
      printf ("not sufficient to continue. Try copying a \"Flat.fits\"\n")
      printf ("file from another night.\n")
      goto cleanup
    }

    AptProc (dir="./", files="@"//FlatFile, useflat-, bias=0, short-)
    imcombine ("@"//FlatFile//"//p", output=t1, combine="average",
        reject="avsigclip",scale="median",weight="median",mclip+,lsigma=6,
        hsigma=6)

    maxpercent = 0.0
    imglist = FlatFile
    while (fscan (imglist, img) != EOF) {
      i = strlen (img)
      if (substr (img, i-4, i) == ".fits") {
        img = substr (img, 1, i-5)
      }
      imdelete (t2//","//t3, ver-, >& "dev$null")
      img = img//"p"
      imarith (img, "/", t1, t2)
      imsurfit (input=t2, output=t3, xorder=3, yorder=3, 
          type_ou="fit", function="spline3", cross+, regions="all")
      imdelete (t2, ver-)   
      imarith (img, "/", t3, img)
      imstatistics (images=t3, fields="min,max", format-) | scan (pmin,pmax)
      percent = 200.0 * (pmax - pmin) / (pmax + pmin)
      printf ("%s has %5.2f%s residual slope\n", img, percent, "%")
      if (percent > maxpercent) {
        maxpercent = percent
      }
    }
    imcombine ("@"//FlatFile//"//p", output=Flat, combine="average",
        reject="avsigclip",scale="none",weight="none",mclip+,lsigma=3,
        hsigma=2)
    imstatistics (images=Flat, fields="mean", format-) | scan (mean)
    imarith (Flat, "/", mean, Flat)
    if (maxpercent > 0.5) {
      printf ("Warning, at least one flat field has a large variation\n")
      goto cleanup
    }
  }

  if (!HaveTwilight) { 
 
     printf ("processing twilight images...\n")

# Make sure there are at least five twilight images:

    imglist = TwilightFile
    i = 0
    while (fscan (imglist, img) != EOF) {
      i += 1
    }
    if (i < 5) {
      printf ("The file \"Twilights\" contains fewer than 5 images, this is\n")
      printf ("not sufficient to continue. Try copying a \"Twilight.fits\"\n")
      printf ("file from another night.\n")
      goto cleanup
    }

# Flatten the twilights using the Flat image. This is used purely for
# the initial stage of masking stars.

    AptProc (dir="./", files="@"//TwilightFile, useflat+, bias=0, 
       short-, flat=Flat)

    imglist = TwilightFile

# Create some temporary files.
  
    t1 = mktemp ("tmp")
    t2 = mktemp ("tmp")
    t3 = mktemp ("tmp")

# For each twilight frame, create a mask to eliminate stars:

    NumTwilights = 0;
    while (fscan (imglist, img) != EOF) {
      NumTwilights += 1;
      printf ("%s ", img)
      i = strlen (img)
      if (substr (img, i-4, i) == ".fits") {
        img = substr (img, 1, i-5)
      }
      imdelete (t1//","//t2//","//t3, ver-, >& "dev$null")
      img = img//"p"

# Fit a smooth surface, and then create a mask showing where the
# twilight image is more than 1.5 stddev above the flat value.

      imsurfit (input=img, output=t1, xorder=4, yorder=4, type_ou="fit",
          function="spline3", cross+, regions="all")
      imarith (img, "/", t1, t2)
      imstatistics (images=t2, fields="mean,stddev", format-) | scan (mean,stddev)
      imexpr ("a < "//str(mean+1.5*stddev)//" ? 0 : 1", t3, t2)
      imdelete (img//"m.pl", ver-, >& "dev$null")
      imcopy (input=t3, output=img//"m.pl", >& "dev$null")
      hedit (images=img, fields="BPM", value=img//"m.pl", add+, delete-, verify-, >& "dev$null")
    }

    imdelete (t1, ver-, >& "dev$null")
    imdelete (t1//".pl", ver-, >& "dev$null")

# Start again with the raw twilight images. This time, don't flat field.

    AptProc (dir="./", files="@"//TwilightFile, useflat-, bias=0, 
             short-)
    imdelete (Twilight//"m.pl", ver-, >& "dev$null")

# Produce a super twilight, called Twilight.fits.

    imcombine ("@"//TwilightFile//"//p", output=Twilight, 
       plfile=Twilight//"m.pl",
       combine="average", masktype="goodval", maskval=0, 
       reject="minmax",scale="median",weight="median", nlow=0, nhigh=2) 
#         reject="avsigclip",scale="median",weight="median",mclip+,lsigma=6,
#        hsigma=3)
    imcopy (Twilight//"m.pl", t1//".pl")
    imdelete (Twilight//"m.pl", ver-, >& "dev$null")
    imexpr ("a != "//str(NumTwilights)//" ? 0 : 1", Twilight//"m.pl",
         t1//".pl")

# Normalize Twilight to 1.00.

    imstatistics (images=Twilight, fields="mean", format-) | scan (mean)
    imarith (Twilight, "/", mean, Twilight)

# Reprocess the raw twilight images using Twilight.fits as the flat.

    AptProc (dir="./", files="@"//TwilightFile, useflat+, bias=0, 
       short-, flat=Twilight)

# For each twilight image...

    imglist = TwilightFile
    while (fscan (imglist, img) != EOF) {
      i = strlen (img)
      if (substr (img, i-4, i) == ".fits") {
        img = substr (img, 1, i-5)
      }
      img = img//"p[*:4,*:4]"

# Fit a smooth surface, and measure the maximum deviation from flatness.
# This can be useful in detecting cloud.

      imdelete (t1, ver-, >& "dev$null")
      imdelete (t1//".pl", ver-, >& "dev$null")

      imsurfit (input=img, output=t1, xorder=5, yorder=5, type_ou="fit",
          function="spline3", cross+, regions="all", upper=20, niter=2)
      imstatistics (images=t1, fields="min,max", format-) | scan (pmin,pmax)
      percent = 200.0 * (pmax - pmin) / (pmax + pmin)
      printf ("%s has %5.2f%s residual variation\n", img, percent, "%")
    }
    printf ("If an image has a variation >0.5%%, perhaps it contains cloud.\n")
    printf ("Please now examine Twilight.fits, compare it with Flat.fits\n")
    printf ("and create SuperFlat.fits (perhaps just using imcopy), then\n")
    printf ("rerun AptReduce to process the files in Stars.\n")
    goto cleanup
  }

  if (!HaveSuperFlat) {
    printf ("The SuperFlat image does not exist.\n")
    printf ("You must create this yourself. A first-order attempt at this\n")
    printf ("would be to use \"imcopy Twilight SuperFlat\".\n")
    goto cleanup
  }

  AptProc (dir="./", files="@"//StarFile, useflat+, bias=0, 
      short-, flat=SuperFlat)

cleanup:
  imdelete (t1, ver-, >& "dev$null")
  imdelete (t1//".pl", ver-, >& "dev$null")
  imdelete (t2, ver-, >& "dev$null")
  imdelete (t2//".pl", ver-, >& "dev$null")
  imdelete (t3, ver-, >& "dev$null")
  imdelete (t3//".pl", ver-, >& "dev$null")
  beep
end