# *****************************************************
procedure satbgon (input, x1, y1, x2, y2, width)
#
# Replace input FITS image with a new version 
# in which all pixels affected by a satellite have been
# replaced with a constant value.
#
# By Matthew Ashby, IRAC IT, 2006 May 26.
# 
# USAGE: satbgon image 
#
# INPUTS: 
#       image:       input image
#       x1:          logical x-pixel value of start of plane trail
#       y1:          logical y-pixel value of start of trail
#       x2:          logical x-pixel value of end of trail
#       y2:          logical y-pixel value of end of trail
#    width:          width in pixels of trail
#    value:          replace pixel with this constant
#
# MODIFIES: The input file is replaced with a modified version.
#
# ASSUMES: 
# *****************************************************

string  input   {prompt="input image ",mode = "q"}
int	x1	{prompt="X-pixel number for start of trail", mode="q"}
int	y1	{prompt="Y-pixel number for start of trail", mode="q"}
int	x2	{prompt="X-pixel number for end of trail", mode="q"}
int	x2	{prompt="Y-pixel number for end of trail", mode="q"}
int	width   {prompt="Width of trail in pixels", mode="q"}
real	value   {prompt="replace pixels with value", mode="q"}

begin

  string infile
  int	xpix, ypix, y_cen, x_cen
  real  x_1, y_1, x_2, y_2, x_bot, x_top, y_bot, y_top, x_dimen, y_dimen
  int   trail_width, trail_radius
  real  new_value
  
  # Print out the version banner.
  # print ("\nThis is PLANEBEGONE version 2006 May 26.")

  # Get the outstanding parameters.
  infile = input
  x_1 = real(x1)
  y_1 = real(y1)
  x_2 = real(x2)
  y_2 = real(y2)
  trail_width = width
  new_value = value 

  # Check validity of inputs
  if(!access(input)) {
    print ("\nERROR: No access to input file.")
    bye
  }
  if (x_1 < 1){
    print ("\nERROR: x1 out of range.")
    print ("Please retry with a value > 1.")
    bye
  }
  if (y_1 < 1){
    print ("\nERROR: y1 out of range.")
    print ("Please retry with a value > 1.")
    bye
  }
  if (x_2 < 1){
    print ("\nERROR: x2 out of range.")
    print ("Please retry with a value > 1.")
    bye
  }
  if (y_2 < 1){
    print ("\nERROR: y2 out of range.")
    print ("Please retry with a value > 1.")
    bye
  }
  if (trail_width < 1){
    print ("\nERROR: width out of range.")
    print ("Please retry with a value > 1.")
    bye
  }

  trail_radius = int(trail_width / 2)

  # Delete this temporary file
  if(access(infile//"_satbgon_pixels")) {
     delete(infile//"_satbgon_pixels",ver-)
  } 

  # Get the exact dimensions of the input file
  imgets(infile,"i_naxis1")
  x_dimen = int(imgets.value)
  imgets(infile,"i_naxis2")
  y_dimen = int(imgets.value)

  # Loop over pixels in range 
  if (abs(x_1 - x_2) >= abs(y_1 - y_2)) {
      if (x_1 < x_2) 
      for (xpix=x_1; xpix <= x_2; xpix+=1) {
          y_cen = y_1 + (xpix - x1) * (y_2 - y_1) / (x_2 - x_1) 
          y_top = y_cen + trail_radius 
          y_bot = y_cen - trail_radius 
          if (y_bot >= 1 && y_top <= y_dimen)
              print (xpix,xpix,y_bot,y_top, >>& infile//"_satbgon_pixels")
          if (y_bot < 1)
              print (xpix,xpix,1,y_top, >>& infile//"_satbgon_pixels")
          if (y_top > y_dimen)
              print (xpix,xpix,y_bot,y_dimen, >>& infile//"_satbgon_pixels")
      }
      if (x_1 > x_2) 
      for (xpix=x_2; xpix <= x_1; xpix+=1) {
          y_cen = y_2 + (xpix - x2) * (y_1 - y_2) / (x_1 - x_2) 
          y_top = y_cen + trail_radius 
          y_bot = y_cen - trail_radius 
          if (y_bot >= 1 && y_top <= y_dimen)
              print (xpix,xpix,y_bot,y_top, >>& infile//"_satbgon_pixels")
          if (y_bot < 1)
              print (xpix,xpix,1,y_top, >>& infile//"_satbgon_pixels")
          if (y_top > y_dimen)
              print (xpix,xpix,y_bot,y_dimen, >>& infile//"_satbgon_pixels")
      }
  } else {
      if (y_1 < y_2) 
          for (ypix=y_1; ypix <= y_2; ypix+=1) {
              x_cen = x_1 + (ypix - y1) * (x_2 - x_1) / (y_2 - y_1) 
              x_top = x_cen + trail_radius 
              x_bot = x_cen - trail_radius 
              if (x_bot >= 1 && x_top <= x_dimen)
                  print (x_bot,x_top,ypix,ypix, >>& infile//"_satbgon_pixels")
              if (x_bot < 1)
                  print (1,x_top,ypix,ypix, >>& infile//"_satbgon_pixels")
              if (x_top > x_dimen)
                  print (x_bot,x_dimen,ypix,ypix, >>& infile//"_satbgon_pixels")
          }
      if (y_1 > y_2) 
          for (ypix=y_2; ypix <= y_1; ypix+=1) {
              x_cen = x_2 + (ypix - y2) * (x_1 - x_2) / (y_1 - y_2) 
              x_top = x_cen + trail_radius 
              x_bot = x_cen - trail_radius 
              if (x_bot >= 1 && x_top <= x_dimen)
                  print (x_bot,x_top,ypix,ypix, >>& infile//"_satbgon_pixels")
              if (x_bot < 1)
                  print (1,x_top,ypix,ypix, >>& infile//"_satbgon_pixels")
              if (x_top > x_dimen)
                  print (x_bot,x_dimen,ypix,ypix, >>& infile//"_satbgon_pixels")
          }
  }

  # Convert the just-created text file to a maskfile
  if(access(infile//"_satbgon_pixels.pl")) {
     delete(infile//"_satbgon_pixels.pl",ver-)
  } 
  text2mask(infile//"_satbgon_pixels", infile//"_satbgon_pixels.pl", ncols=x_dimen, nlines=y_dimen)
  if(access(infile//"_satbgon_pixels.fits")) {
     delete(infile//"_satbgon_pixels.fits",ver-)
  } 
  imcopy(infile//"_satbgon_pixels.pl", infile//"_satbgon_pixels.fits")

  if (access("tmp.fits"))
      delete("tmp.fits",ver-)
  imexpr ("b==0?a : -10000", "tmp.fits", infile, infile//"_satbgon_pixels.fits")
  delete (infile, ver-)
  imrename ("tmp.fits", infile)

   # Clean up
   if(access(infile//"_satbgon_pixels")) {
      delete(infile//"_satbgon_pixels",ver-)
   } 
   if(access(infile//"_satbgon_pixels.pl")) {
      delete(infile//"_satbgon_pixels.pl",ver-)
   } 
   if(access(infile//"_satbgon_pixels.fits")) {
      delete(infile//"_satbgon_pixels.fits",ver-)
   } 

end