procedure satzap( object, x1, y1, x2, y2 )

# procedure to mask out satellite trails in images. Note that 
# this procedure runs on an object frame, but modifies the 
# BPM file named by the BPM header keyword! 
#
# Takes about 90 seconds / projected MOSAIC frame on a Ultra-60
#
# Written by James Rhoads
# Modified by Arjun Dey to allow partial trails to be flagged
# 08/04/99

string object  { prompt="Input image"}
real x1 { prompt="Satellite trail: location 1 - X "}
real y1 { prompt="Satellite trail: location 1 - Y "}
real x2 { prompt="Satellite trail: location 2 - X "}
real y2 { prompt="Satellite trail: location 2 - Y "}
real hwidth=15. {prompt="half-width of flagged track?"}

begin
   real l_hwide, dx12, dy12, d12, theta, ctheta, stheta
   real x_min, y_min, x_max, y_max
   string pdexp, ypr, minmaxt
   l_hwide = hwidth

   dx12 = x1-x2
   dy12 = y1-y2
   d12 = sqrt( (x1-x2)**2 + (y1-y2)**2 )

   theta = atan2( (y2-y1),(x2-x1) )
   ctheta = cos(theta)
   stheta = sin(theta)

   # find minimum and max x & y values: 
   x_min = x1-hwidth-1.
   x_max = x2+hwidth+1.
   if (x2 < x1) { 
	x_min = x2-hwidth-1.
	x_max = x1+hwidth+1.
   }
   y_min = y1-hwidth-1.
   y_max = y2+hwidth+1.
   if (y2 < y1) { 
	y_min = y2-hwidth-1.
	y_max = y1+hwidth+1.
   }
	
   ypr = "(y-"//str(y1)//")*"//str(ctheta)//" - (x-"//str(x1)//")*"//str(stheta)
   minmaxt="((x>="//str(x_min)//")&&(x<="//str(x_max)//")&&(y>="//str(y_min)//")&&(y<="//str(y_max)//"))"
   pdexp = "if ((abs("//ypr//") < "//str(hwidth)//")&&"//minmaxt//") then im1+1 else im1"

   print( pdexp )
   imgets( object, "BPM" )
   imcalc( imgets.value, imgets.value, pdexp, verbose+ )

end