;+
;-------------------------------------------------------------------------------
; NAME		fcc.pro
;
; PURPOSE	Display a FITS image & perform cavity measurements.
;
; SYNTAX	fc, fits_image, cm='colormap.lut', wmin=1, wmax=240, /gif
;
; PARAMETERS
;		fits_image	filename of Mk4 FITS image.
;		cm		pathname of ASCII colormap file.
;				Each line has the syntax: index red green blue
;				where index = 0, 1, 2, ... 255,
;				and red/green/blue are in the range [0:255].
;		wmin		display minimum value.
;		wmax		display maximum value.
;		gif		write octant image & plots as GIF files.
;		nolabel		If set, do NOT display the position # label.
;               notext    	If set, do NOT write position information
;				to a text file.
;
; EXAMPLES	fc, '19981101.1234.mk4.rpb.fts', $
;		       cm='/home/stanger/color/bwy.lut'
; 		fc, '19981101.1234.mk4.rpb.fts', /gif
; 		fc, '19981101.1234.mk4.rpb.fts', /text
; 		fc, '19981101.1234.mk4.rpb.fts', wmin=0, wmax=2500 
; 		fc, '19981101.1234.mk4.rpb.fts', /nolabel
;
; EXTERNAL
; PROCEDURES
;		readfits		read FITS image
;		headfits		read FITS header
;		fxpar			read FITS keyword parameter
;		fitsdispw 		display  FITS image
;		fits_annotatew.pro	annotate FITS image
;		cavity_width.pro	position/label procedure
;
; DEFINITION	PA			position angle: CCW from solar north.
;
; HISTORY	Andrew L. Stanger   HAO/NCAR   14 September 2001
;		26 Sep 2001: [ALS] Add cursor position.
;		17 Oct 2001: [ALS] If "CRRADIUS" not present, use "SOLAR_R".
;		14 Jan 2003: [ALS] Reverse roll direction for mouse_pos_lab.
;		19 Sep 2003: [ALS] disp_label option added.
;		 7 Oct 2004: [ALS] Replace mouse_pos_lab with pos_lab_time.pro.
;		14 Oct 2004: [ALS] Modify for cavity measurements.
;-------------------------------------------------------------------------------
;-

PRO fcc, fits_name, cm=cm, wmin=wmin, wmax=wmax,	$
	 gif=gif, notext=notext, nolabel=nolabel

   ;--- Define variables.

   gif  = 1
   xlab = ''
   ylab = ''
   type_obs = 'model'
   xb = 160 			; x-axis border for annotation.
   yb =   0			; y-axis border for annotation.
   xdim_prev = 0		; previous image x-axis dimension.
   ydim_prev = 0		; previous image y-axis dimension.


   ;--- Color table: annotation colors.

   ov_white  = 255
   ov_red    = 254
   ov_green  = 253
   ov_blue   = 252
   ov_grey   = 251
   ov_yellow = 250
   ov_black  =   0

   ;--- Set text file flag.

   IF (KEYWORD_SET (notext)) THEN		$
         text = 0				$
   ELSE	 text = 1

   ;--- Set display label flag.

   IF (KEYWORD_SET (nolabel)) THEN		$
         disp_label = 0				$
   ELSE	 disp_label = 1

   ;--- Set character size.

   SET_PLOT, 'X'
   WINDOW, xs=10, ys=10
   WDELETE, !d.window
   DEVICE, SET_CHARACTER_SIZE = [8, 12]

   SET_PLOT, 'X'

   ;----------------------------------------------------------------------------
   ; Load color table.
   ;----------------------------------------------------------------------------

   IF (KEYWORD_SET (cm))  THEN 		$
   BEGIN
      PRINT, 'cm: ', cm
      dirend = -1

      FOR i = 0, strlen (cm) - 1 DO	$
      BEGIN
         dirloc = STRPOS (cm, '/', i)
         IF (dirloc GE 0)  THEN dirend = dirloc
      END

      IF (dirend NE -1)  THEN		$
      BEGIN
         PRINT, 'dirend: ', dirend
         coldir = STRMID (cm, 0, dirend)
         PRINT, 'coldir: ', coldir
         ccm = STRMID (cm, dirend+1, strlen (cm) - dirend - 1)
         PRINT, 'ccm: ', ccm
      END

      IF (dirend EQ -1)  THEN				$
         lct, '/home/cordyn/color/' + cm + '.lut'	$
      ELSE						$
         lct, cm				; Load specified colormap.
   END					$
   ELSE  lct, '/home/cordyn/color/quallab.lut'	; Load default colormap.
 
   red   = bytarr (256)
   green = bytarr (256)
   blue  = bytarr (256)
   TVLCT, red, green, blue, /GET	; Fetch RGB color look-up tables.

   ;--- Extract base name from specified FITS file name.

   ftspos = STRPOS (fits_name, '.fts')
   basename = STRMID (fits_name, 0, ftspos)
   print, 'basename: ', basename

   IF (text) THEN	$
   BEGIN
      pfile = basename + '.pos'
      CLOSE,  21
      OPENW,  21, pfile
      PRINTF, 21, fits_name, '   Position Measurement[s]'
      CLOSE,  21
   END

   ;----------------------------------------------------------------------------
   ; Read FITS image header & pixel data.
   ;----------------------------------------------------------------------------

   img = readfits (fits_name, hdu)

   ;--- Extract information from header.

   xdim     = fxpar (hdu, 'NAXIS1')
   ydim     = fxpar (hdu, 'NAXIS2')

   ;--- Get information from FITS header.

   orbit_id = fxpar (hdu, 'ORBIT-ID')
   image_id = fxpar (hdu, 'IMAGE-ID')
   telescop = fxpar (hdu, 'TELESCOP')
   instrume = fxpar (hdu, 'INSTRUME')
   dateobs  = fxpar (hdu, 'DATE-OBS')
   timeobs  = fxpar (hdu, 'TIME-OBS')
   typeobs  = fxpar (hdu, 'TYPE-OBS', count=count_type_obs)
   dataform = fxpar (hdu, 'DATAFORM')
   rsun     = fxpar (hdu, 'RSUN')
   bunit    = fxpar (hdu, 'BUNIT')
   datamin  = fxpar (hdu, 'DATAMIN')
   datamax  = fxpar (hdu, 'DATAMAX')
   dispmin  = fxpar (hdu, 'DISPMIN')
   dispmax  = fxpar (hdu, 'DISPMAX')
   crradius = fxpar (hdu, 'CRRADIUS', count=count_crradius)
   solar_r  = fxpar (hdu, 'SOLAR_R',  count=count_solar_r)

   srsun    = STRING (rsun, FORMAT='(F7.2)')
   print, 'count_crradius: ', count_crradius, '  crradius: ', crradius

   IF (count_type_obs EQ 1) THEN	$
      type_obs = typeobs		$
   ELSE					$
      type_obs = 'model'

   IF (count_crradius EQ 1)  THEN	$
      pixrs = crradius			$
   ELSE					$
   IF (count_solar_r EQ 1) THEN		$
      pixrs = solar_r

;   print, 'count_crradius: ', crradius, '  count_solar_r: ', count_solar_r
;   print, 'crradius: ', crradius
;   print, 'solar_r : ', solar_r
;   print, 'pixrs   : ', pixrs

   xcen    = fxpar (hdu, 'CRPIX1')
   ycen    = fxpar (hdu, 'CRPIX2')
   roll    = fxpar (hdu, 'CROTA1')
   revroll = -roll
   xcenf   = xcen + xb
   ycenf   = ycen + yb

   ;--- Resize window [if image size has changed].

   if (xdim NE xdim_prev OR ydim NE ydim_prev)  THEN	$
   BEGIN
       WINDOW, 1, xsize=xdim+xb, ys=ydim+yb, retain=2
       WSET,   1
   END

   print, 'xdim + xb: ', xdim + xb
   print, 'ydim + yb: ', ydim + yb

   xdim_prev = xdim
   ydim_prev = ydim

   ;--- Annotate image.

;   fits_annotatew, hdu, xdim, ydim, xb, yb

   ;----------------------------------------------------------------------------
   ; Display original image.
   ;----------------------------------------------------------------------------

   dmin = dispmin
   dmax = dispmax
;   IF (KEYWORD_SET (wmin) OR wmin EQ 0.0)  THEN dmin = wmin
   IF (KEYWORD_SET (wmin))  THEN dmin = wmin
   IF (KEYWORD_SET (wmax))  THEN dmax = wmax

;   TV, BYTSCL (img, min=dispmin, max=dispmax, top=249), xb, yb

   fitsdispw, fits_name, xdim_prev, ydim_prev, gif=gif, wmin=dmin, wmax=dmax, $
	      xb, yb

   ;--- Get approximate cavity position angle.

   PRINT,  'Click left mouse button at approximate position of the cavity.'
   CURSOR, rxcur, rycur, /DEVICE, WAIT=3
   ierror = rcoord (radius, angle, rxcur, rycur, -1, roll, xcenf, ycenf, pixrs)

   IF (angle LT   0.0) THEN angle = angle + 360.0
   IF (angle GT 360.0) THEN angle = angle - 360.0

   ;--- Determine which octant to display.

   IF (angle GE 337.5 AND angle LE 360.0)  THEN oct = 1	 ;   0 deg=octant middle
   IF (angle GE   0.0 AND angle LE  22.5)  THEN oct = 1	 ;   0 deg
   IF (angle GE  22.5 AND angle LT  77.5)  THEN oct = 2	 ;  45 deg
   IF (angle GE  77.5 AND angle LT 112.5)  THEN oct = 3	 ;  90 deg
   IF (angle GE 112.5 AND angle LT 157.5)  THEN oct = 4	 ; 135 deg
   IF (angle GE 157.5 AND angle LT 202.5)  THEN oct = 5	 ; 180 deg
   IF (angle GE 202.5 AND angle LT 247.5)  THEN oct = 6	 ; 225 deg
   IF (angle GE 247.5 AND angle LT 292.5)  THEN oct = 7	 ; 270 deg
   IF (angle GE 292.5 AND angle LT 337.5)  THEN oct = 8	 ; 315 deg

   print, 'oct: ', oct

   ;--- Octant image array center & dimension.

   ixcen = xdim / 2
   iycen = ydim / 2
   oxdim = xdim / 2
   oydim = ydim / 2

   ;--- Compute offset of sun center from array center.

   xcenoff = xcen - ((xdim / 2) - 0.5)
   ycenoff = ycen - ((ydim / 2) - 0.5)

   ascanmin =  0.0			; start angle for azimuthal scans.
   ascanmax = 40.0			; stop  angle

   ;--- Specify the lower left corner (x,y) of the subset image (octrant).

   CASE oct OF
   1:    BEGIN	;   0 degrees octant
	     xmin = xdim / 4
	     ymin = ydim / 2
	     oxcen = oxdim / 2 - 0.5
	     oycen = - 0.5
	     ascanmin = 340.0
	     ascanmax =  20.0 + 360.0
         END
   2:    BEGIN	;  45 degrees octant
	     xmin = 0
	     ymin = ydim / 2
	     oxcen = oxdim + 0.5
	     oycen = - 0.5
	     ascanmin =  25.0
	     ascanmax =  65.0
         END
   3:    BEGIN	;  90 degrees octant
	     xmin = 0
	     ymin = ydim / 4
	     oxcen = oxdim + 0.5
	     oycen = oydim / 2 - 0.5
	     ascanmin =  70.0
	     ascanmax = 110.0
         END
   4:    BEGIN	; 135 degrees octant
	     xmin = 0
	     ymin = 0
	     oxcen = oxdim + 0.5
	     oycen = oydim + 0.5
	     ascanmin = 115.0
	     ascanmax = 155.0
         END
   5:    BEGIN	; 180 degrees octant
	     xmin = xdim / 4
	     ymin = 0
	     oxcen = oxdim / 2 - 0.5
	     oycen = oydim + 0.5
	     ascanmin = 160.0
	     ascanmax = 200.0
         END
   6:    BEGIN	; 225 degrees octant
	     xmin = xdim / 2
	     ymin = 0
	     oxcen = - 0.5
	     oycen = oydim + 0.5
	     ascanmin = 205.0
	     ascanmax = 245.0
         END
   7:    BEGIN	; 270 degrees octant
	     xmin = xdim / 2
	     ymin = ydim / 4
	     oxcen = - 0.5
	     oycen = oydim / 2 - 0.5
	     ascanmin = 250.0
	     ascanmax = 290.0
         END
   8:    BEGIN	; 335 degrees octant
	     xmin = xdim / 2
	     ymin = ydim / 2
	     oxcen = - 0.5
	     oycen = - 0.5
	     ascanmin = 295.0
	     ascanmax = 335.0
         END
   ELSE: BEGIN
	     PRINT, 'Invalid octant: ', oct
	     RETURN
	 END
   ENDCASE

   ;--- Specify the upper right corner (x,y) of the subset image (octrant).

   ;--- Note: subset image (octant) to be extracted is 1/2 the size
   ;---       of the original image.

   xmax = xmin + xdim / 2 - 1
   ymax = ymin + ydim / 2 - 1

   print, 'xmin, ymin: ', xmin, ymin
   print, 'xmax, ymax: ', xmax, ymax
   print, 'oxcen, oycen: ', oxcen, oycen

   ;--- Adjust subset image array center so that it is now sun center.

   oxcen = oxcen + xcenoff
   oycen = oycen + ycenoff

   ;----------------------------------------------------------------------------
   ; Extract octant image.
   ;----------------------------------------------------------------------------

   octimg = img [xmin:xmax, ymin:ymax]

   ;--- Double the size of the octant image.

   sxdim = xdim
   sydim = ydim
   sxcen = oxcen * 2.0
   sycen = oycen * 2.0
   simg = REBIN (octimg, sxdim, sydim)
   spixrs = pixrs * 2.0

   print, 'dmin/dmax: ', dmin, dmax
   simg = BYTSCL (simg, min=dmin, max=dmax, top=249)

   imin = MIN (octimg, max=imax)
   print, 'octimg min/max: ', imin, imax

   ;----------------------------------------------------------------------------
   ; Display octant image.
   ;----------------------------------------------------------------------------

   WDELETE, 1
   WINDOW,  2, xsize=sxdim, ysize=sydim, retain=2
   WSET,    2
;   lct, '/home/cordyn/color/quallab.lut'
   TV, simg

   ;--- Get the polar coordinates of the top of the cavity.

   PRINT,  'Click left mouse button at the TOP EDGE of the cavity.'
   CURSOR, xtop, ytop, /DEVICE, WAIT=3

   ;--- Convert position to polar coordinates [radius, angle].

   ierror = rcoord (radius, angle, xtop, ytop, -1, roll, sxcen, sycen, spixrs)
   print, '[x,y]:',  xtop, ytop, '  [r,th]: ', radius, angle

   ;--- Mark the position of the cavity top on the octant image.

   pos = 1
   draw_pos, pos, xtop, ytop, 3, ov_red, cs
   pos = pos + 1

   ;--- Ensure that angle is between 0 and 360 degrees.

   IF (angle LT   0.0) THEN angle = angle + 360.0
   IF (angle GT 360.0) THEN angle = angle - 360.0

   ;--- Write the cavity height information to the text file.

   CLOSE,  21
   OPENU,  21, pfile, /APPEND
   PRINTF, 21, 'Cavity Height: ', radius, ' Rsun  Position Angle: ', angle, $
	       ' degrees'
   PRINTF, 21, 'pixel cartesian coordinates: ', xtop, ytop
   PRINTF, 21, 'pixel magnitude:             ', simg [xtop, ytop]
   CLOSE,  21

   ;--- Compute 1/3 & 2/3 cavity heights.

   radtop = radius
   angtop  = angle
   radmin  = 1.1	; Assume inner FOV limit is the bottom of the cavity.
   rerr =  (radtop - radmin) * 0.15			; radial tolerance
   r1   = ((radtop - radmin) * 0.33333333) + radmin     ; 1/3 height
   r2   = ((radtop - radmin) * 0.66666667) + radmin     ; 2/3 height

   ;----------------------------------------------------------------------------
   ; Get cavity width [use mouse cursor].
   ;----------------------------------------------------------------------------

   cw, sxdim, sydim, sxcen, sycen, spixrs, revroll,		$
	         rpos, apos, xpos, ypos, ascanmin, ascanmax,		$
		 radtop, angtop, r1, r2,				$
		 numpos, pos=pos, pfile=pfile,				$
		 disp_label=disp_label

   ;--- Determine which position measurements belong to each azimuthal arc.

   print, 'numpos: ', numpos
   print, 'rpos:   ', rpos
   print, 'apos:   ', apos

;   cavmin =    360.0
;   cavmax =      0.0
;   cavtop =      0.0
;
;   FOR i = 0, numpos - 1  DO	$
;   BEGIN
;      print, 'i: ', i, ' rpos: ', rpos [i], ' apos: ', apos [i]
;      IF (apos [i] LT cavmin) THEN $
;	 BEGIN
;	    cavmin = apos [i]
;	    minpos = i
;         END
;      IF (apos [i] GT cavmax) THEN $
;	 BEGIN
;	    cavmax = apos [i]
;	    maxpos = i
;         END
;      IF (rpos [i] GT cavtop) THEN $
;	 BEGIN
;	    cavtop = rpos [i]
;	    toppos = i
;         END
;   END

   amax1 = -1000.0
   amax2 = -1000.0
   amin1 =  1000.0
   amin2 =  1000.0
   ss = FLTARR (numpos)			; Indicates which scan height was used.

   ;--- Identify which scan height belongs to each cavity position measurement.

   FOR i = 0, numpos - 1  DO	$
   BEGIN
      IF (abs (rpos [i] - r1) LT rerr) THEN ss [i] = 1 
      IF (abs (rpos [i] - r2) LT rerr) THEN ss [i] = 2
   END

   FOR i = 0, numpos - 1 DO	$
   BEGIN
      ;--- Identify the max PA cavity boundary angle for the lower height.

      IF (ss [i] EQ 1 AND apos [i] GT amax1) THEN	$
      BEGIN
         amax1 = apos [i]
	 rmax1 = rpos [i]
	 sxmax1 = xpos [i]
	 symax1 = ypos [i]
      END

      ;--- Identify the min PA cavity boundary angle for the lower height.

      IF (ss [i] EQ 1 AND apos [i] LT amin1) THEN	$
      BEGIN
         amin1 = apos [i]
	 rmin1 = rpos [i]
	 sxmin1 = xpos [i]
	 symin1 = ypos [i]
      END

      ;--- Now identify the max PA cavity boundary angle for the upper height.

      IF (ss [i] EQ 2 AND apos [i] GT amax2) THEN	$
      BEGIN
         amax2 = apos [i]
	 rmax2 = rpos [i]
	 sxmax2 = xpos [i]
	 symax2 = ypos [i]
      END

      ;--- Now identify the min PA cavity boundary angle for the upper height.

      IF (ss [i] EQ 2 AND apos [i] LT amin2) THEN	$
      BEGIN
         amin2 = apos [i]
	 rmin2 = rpos [i]
	 sxmin2 = xpos [i]
	 symin2 = ypos [i]
      END

   END

   print, 'rmin1/amin1: ', rmin1, amin1
   print, 'rmax1/amax1: ', rmax1, amax1
   print, 'rmin2/amin2: ', rmin2, amin2
   print, 'rmax2/amax2: ', rmax2, amax2

   ;--- Convert cavity edge positions from polar to cartesian coordinates.
   
   ierror = rcoord (r1, amin1, xmin1, ymin1,  1, roll, xcen, ycen, pixrs)
   ierror = rcoord (r1, amax1, xmax1, ymax1,  1, roll, xcen, ycen, pixrs)
   ierror = rcoord (r2, amin2, xmin2, ymin2,  1, roll, xcen, ycen, pixrs)
   ierror = rcoord (r2, amax2, xmax2, ymax2,  1, roll, xcen, ycen, pixrs)

   ;--- Write information to text file.

   CLOSE,  21
   OPENU,  21, pfile, /APPEND
   PRINTF, 21, 'Lower Height:     ', rmin1, ' Rsun'
   PRINTF, 21, 'Min/Max angle:    ', amin1, amax1, ' degrees'
   PRINTF, 21, 'Pixel magnitudes: ', simg [sxmin1,symin1], simg [sxmax1,symax1]
   PRINTF, 21, 'Pixel magnitudes: ', img  [ xmin1, ymin1],  img [ xmax1, ymax1]
   PRINTF, 21, 'Upper Height:     ', rmin2, ' Rsun'
   PRINTF, 21, 'Min/Max angle:    ', amin2, amax2, ' degrees'
   PRINTF, 21, 'Pixel magnitudes: ', simg [sxmin2,symin2], simg [sxmax2,symax2]
   PRINTF, 21, 'Pixel magnitudes: ', img  [ xmin2, ymin2], img  [ xmax2, ymax2]
   CLOSE,  21

   ;--- Compute center of cavity [angle].
   ;    Average the centers found from the two scan heights.

   angcen1 = (amax1 - amin1) * 0.5 + amin1
   angcen2 = (amax2 - amin2) * 0.5 + amin2
   angcen = (angcen1 + angcen2) * 0.5
   sangcen = STRTRIM (STRING (angcen, FORMAT='(F7.2)'), 2)

   ;--- Set scan parameters.

   anginc = 0.5
   radbeg = 1.1
   radend = 2.0
   radinc = 0.01
   radmin = 1.1
   radsun = 1.0

   srad1   = STRTRIM (STRING (r1, FORMAT='(F5.2)'), 2)
   srad2   = STRTRIM (STRING (r2, FORMAT='(F5.2)'), 2)

;   print, 'angcen: ', angcen

   ;----------------------------------------------------------------------------
   ; Radial scan at middle of cavity.
   ;----------------------------------------------------------------------------

   rscan, fits_name, img, pixrs, revroll, xcen, ycen,	$
	  radbeg, radend, radinc, angcen,	$
	  rs,  rsdx,  n3

   ;--- Smooth scan data.

   rs_ave = SMOOTH (rs, 3)
   rs_min = MIN (rs_ave)
   rs_max = MAX (rs_ave)

   ;--- Draw "X" symbols to mark the radial scan path on the octant image.

   WSET, 2
;   print, 'radial scan:'
   angle = angcen
   FOR i = 0, n3-1 DO			$
   BEGIN
     radius = radbeg + i * radinc
     ierror = rcoord (radius, angle, xpos, ypos,  1, roll, sxcen, sycen, spixrs)
;     print, '[x,y]:',  xpos, ypos, '  [r,th]: ', radius, angle
     PLOTS, [xpos-1, xpos+1], [ypos-1, ypos+1], /DEVICE, color=ov_green
     PLOTS, [xpos-1, xpos+1], [ypos+1, ypos-1], /DEVICE, color=ov_green
   END

   ;--- Draw azimuthal arc at top of cavity.

   WSET, 2
   print, 'theta scan at top of cavity:'
   angbeg = angcen - 20.0
   anginc = 1.0
   FOR i = 0, 40 DO			$
   BEGIN
     angle = angbeg + i * anginc
     ierror = rcoord (radtop, angle, xpos, ypos,  1, roll, sxcen, sycen, spixrs)
;     print, '[x,y]:',  xpos, ypos, '  [r,th]: ', radtop, angle
     PLOTS, [xpos-1, xpos+1], [ypos-1, ypos+1], /DEVICE, color=ov_red
     PLOTS, [xpos-1, xpos+1], [ypos+1, ypos-1], /DEVICE, color=ov_red
   END

   ;--- Write displayed octant image as a GIF file (if "gif" keyword is set).

   IF (KEYWORD_SET (gif))  THEN		$
   BEGIN
      gif_file = basename + '.cav.gif'
      img_gif = TVRD ()
      WRITE_GIF, gif_file, img_gif, red, green, blue
   END

   ;----------------------------------------------------------------------------
   ; Azimuthal scan at lower height.
   ;----------------------------------------------------------------------------

   tscan, fits_name, img, pixrs, revroll, xcen, ycen,	$
	  ascanmin, ascanmax, anginc, r1,		$
	  ts1, ts1dx, n1

   ;--- Scan only inside the cavity.

   tscan, fits_name, img, pixrs, revroll, xcen, ycen,	$
	  amin1, amax1, anginc, r1,		$
	  ts11, ts11dx, n11

   ts11_ave = SMOOTH (ts11, 3)				; Smooth theta scan.

   ;--- Find minimum intensity along theta scan inside cavity.

   pixlow1 = 1.0e+20
   FOR i = 0, n11 - 1 DO		$
   BEGIN
      IF (ts11 [i] LT pixlow1) THEN	$
      BEGIN
	 pixlow1 = ts11_ave [i]
	 anglow1 = ts11dx   [i]
      END
   END

   ;--- Compute cartesian coordinates for minimum intensity at lower height.
   ;--- Use original image.

   ierror = rcoord (r1, anglow1, xlow1, ylow1,  1, roll, xcen, ycen, pixrs)
   ixlow1 = FIX (xlow1 + 0.5)
   iylow1 = FIX (ylow1 + 0.5)

   ;--- Compute cartesian coordinates for minimum intensity at lower height.
   ;--- Use subset image.

   ierror = rcoord (r1, anglow1, sxlow1, sylow1,  1, roll, sxcen, sycen, spixrs)
   isxlow1 = FIX (sxlow1 + 0.5)
   isylow1 = FIX (sylow1 + 0.5)
   draw_pos, pos, isxlow1, isylow1, 3, ov_red, cs
   pos = pos + 1

   ;--- Smooth scan data.

   ts1_ave = SMOOTH (ts1, 3)
   ts1_min = MIN (ts1_ave)
   ts1_max = MAX (ts1_ave)

   ;----------------------------------------------------------------------------
   ; Azimuthal scan at upper height.
   ;----------------------------------------------------------------------------

   tscan, fits_name, img, pixrs, revroll, xcen, ycen,	$
	  ascanmin, ascanmax, anginc, r2,		$
	  ts2, ts2dx, n2

   ;--- Smooth scan data.

   ts2_ave = SMOOTH (ts2, 3)
   ts2_min = MIN (ts2_ave)
   ts2_max = MAX (ts2_ave)

   ;--- Scan only inside the cavity.

   tscan, fits_name, img, pixrs, revroll, xcen, ycen,	$
	  amin2, amax2, anginc, r2,		$
	  ts22, ts22dx, n22

   ts22_ave = SMOOTH (ts22, 3)				; Smooth theta scan.

   ;--- Find minimum intensity along theta scan inside cavity.

   pixlow2 = 1.0e+20
   FOR i = 0, n22 - 1 DO		$
   BEGIN
      IF (ts22 [i] LT pixlow2) THEN	$
      BEGIN
	 pixlow2 = ts22_ave [i]
	 anglow2 = ts22dx   [i]
      END
   END

   ;--- Compute cartesian coordinates for minimum intensity at upper height.
   ;--- Use original image.

   ierror = rcoord (r2, anglow2, xlow2, ylow2,  1, roll, xcen, ycen, pixrs)
   ixlow2 = FIX (xlow2 + 0.5)
   iylow2 = FIX (ylow2 + 0.5)

   ;--- Compute cartesian coordinates for minimum intensity at lower height.
   ;--- Use subset image.

   ierror = rcoord (r2, anglow2, sxlow2, sylow2,  1, roll, sxcen, sycen, spixrs)
   isxlow2 = FIX (sxlow2 + 0.5)
   isylow2 = FIX (sylow2 + 0.5)
   draw_pos, pos, isxlow2, isylow2, 3, ov_red, cs
   pos = pos + 1

   ;--- Write displayed octant image as a GIF file (if "gif" keyword is set).

   IF (KEYWORD_SET (gif))  THEN		$
   BEGIN
      gif_file = basename + '.cav.gif'
      img_gif = TVRD ()
      WRITE_GIF, gif_file, img_gif, red, green, blue
   END

   ;--- Write minimum intensity value & position to text file.

   CLOSE,  21
   OPENU,  21, pfile, /APPEND
   PRINTF, 21, 'Height ', r1
   PRINTF, 21, 'Minimum pixel intensity: ', img [ixlow1, iylow1]
   PRINTF, 21, 'Minimum pixel intensity: ', pixlow1
   PRINTF, 21, 'Minimum pixel angle:     ', anglow1
   PRINTF, 21, 'Height ', r2
   PRINTF, 21, 'Minimum pixel intensity: ', img [ixlow2, iylow2]
   PRINTF, 21, 'Minimum pixel intensity: ', pixlow2
   PRINTF, 21, 'Minimum pixel angle:     ', anglow2
   CLOSE,  21

;   WDELETE, 2

   ;----------------------------------------------------------------------------

   SET_PLOT, 'Z'
   SET_PLOT, 'X'
   WINDOW, 3, xs=640, ys=480, retain=2
;   LOADCT, 0
   fcol = ov_black
   bcol = ov_white
   !P.BACKGROUND=bcol
   print, 'type_obs: ', type_obs
   print, 'bunit: ', bunit
   ylab = type_obs + ' [' + bunit + ']'
   print, 'ylab: ', ylab

;   print, 'cavmin/cavmax:   ', cavmin,  cavmax
   print, 'ts1_min/ts1_max: ', ts1_min, ts1_max

   ;----------------------------------------------------------------------------
   ; Plot azimuthal scan at first height.
   ;----------------------------------------------------------------------------

   PLOT, ts1dx, ts1_ave,						$
	 title='Theta Scan ' + fits_name + ' ' + srad1 + ' Rsun',	$
	 xtitle = 'Position Angle',					$
	 ytitle = ylab,							$
	 color=fcol
   PLOTS, [amin1, amin1], [ts1_min, ts1_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [amax1, amax1], [ts1_min, ts1_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [amin1, amax1], [ts1_min, ts1_min], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [amin1, amax1], [ts1_max, ts1_max], /DATA, LINESTYLE=2, color=fcol
   XYOUTS, (amax1-amin1)*0.3 + amin1, ts1_min * 1.1, 'CAVITY', $
           color=254, charsize=cs, /DATA

   ;--- Use cursor to get intensity maxima of cavity legs.

   ;--- Left leg maximum.

   PRINT, 'Click left mouse button at the MAXIMUM of the left cavity leg.'
   CURSOR, angleg1l, pixleg1l, /DATA, WAIT=3
   WAIT, 0.5

   ;--- Find the scan index corresponding to the intensity maxima of the legs.

   print, 'angleg1l: ', angleg1l
   ii = -1
   FOR i = 0, n1-1 DO			$
   BEGIN
      IF (ts1dx [i] LE angleg1l) THEN ii = i
      print, 'ts1dx: ', ts1dx [i], ' i: ', i
   END
   print, 'ii: ', ii

   IF (ii GE 0 AND ii LE n2-1) THEN	$
   BEGIN
      aleg1l = ts1dx   [ii]				; angle
      pleg1l = ts1_ave [ii]				; scan intensity
      print, 'aleg1l: ', aleg1l
      print, 'pleg1l: ', pleg1l
   END					$
   ELSE					$
   BEGIN
      print, 'Invalid left leg position: ', angleg1l
      print, 'Expected range: ', ts1dx [0], ts2dx [n1-1]
   END

   ;--- Draw "+" to mark left leg maximum on azimuthal plot.

   PLOTS, [angleg1l, angleg1l], [0.9*pleg1l,  1.1*pleg1l],  /DATA, color=ov_red
   PLOTS, [angleg1l-3, angleg1l+3], [pleg1l,  pleg1l],  /DATA, color=ov_red

   ;--- Right leg maximum.

   PRINT, 'Click left mouse button at the MAXIMUM of the right cavity leg.'
   CURSOR, angleg1r, pixleg1r, /DATA, WAIT=3
   WAIT, 0.5

   ii = -1
   FOR i = 0, n1-1 DO			$
   BEGIN
      IF (ts1dx [i] LE angleg1r) THEN ii = i
   END

   IF (ii GE 0 AND ii LE n1-1) THEN	$
   BEGIN
      aleg1r = ts1dx   [ii]				; angle
      pleg1r = ts1_ave [ii]				; scan intensity
      print, 'aleg1r: ', aleg1r
      print, 'pleg1r: ', pleg1r
   END					$
   ELSE					$
   BEGIN
      print, 'Invalid right leg position: ', angleg1r
      print, 'Expected range: ', ts1dx [0], ts1dx [n1-1]
   END

   ;--- Draw "+" to mark right leg maximum on azimuthal plot.

   PLOTS, [angleg1r, angleg1r], [0.9*pleg1r,  1.1*pleg1r],  /DATA, color=ov_red
   PLOTS, [angleg1r-3, angleg1r+3], [pleg1r,  pleg1r],  /DATA, color=ov_red

   ;--- Draw "+" to mark minimum intensity position on azimuthal plot.

   PLOTS, [anglow1,  anglow1 ], [0.9*pixlow1, 1.1*pixlow1], /DATA, color=ov_red
   PLOTS, [anglow1 -3, anglow1 +3], [pixlow1, pixlow1], /DATA, color=ov_red

   ;--- Save azimuthal scan plot [lower height] as a GIF file.

   IF (KEYWORD_SET (gif))  THEN		$
   BEGIN
      gif_file = basename + '.ts1.gif'
      img_gif = TVRD ()
      WRITE_GIF, gif_file, img_gif, red, green, blue
   END

   ;--- Write leg maxima information to the text file.

   CLOSE,  21
   OPENU,  21, pfile, /APPEND
   PRINTF, 21, 'Height: ', r1
   PRINTF, 21, 'Left  leg   angle: ', angleg1l, ' magnitude: ', pixleg1l
   PRINTF, 21, 'Right leg   angle: ', angleg1r, ' magnitude: ', pixleg1r
   CLOSE,  21

   ;----------------------------------------------------------------------------
   ; Plot azimuthal scan at upper height.
   ;----------------------------------------------------------------------------

   PLOT, ts2dx, ts2_ave,						$
	 title='Theta Scan ' + fits_name + ' ' + srad2 + ' Rsun',	$
	 xtitle = 'Position Angle',					$
	 ytitle = ylab,							$
	 color=fcol
   PLOTS, [amin2, amin2], [ts2_min, ts2_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [amax2, amax2], [ts2_min, ts2_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [amin2, amax2], [ts2_min, ts2_min], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [amin2, amax2], [ts2_max, ts2_max], /DATA, LINESTYLE=2, color=fcol
   XYOUTS, (amax2-amin2)*0.3 + amin2, ts2_min * 1.1, 'CAVITY', $
           color=254, charsize=cs, /DATA

   ;--- Use cursor to get intensity maxima of cavity legs.

   ;--- Left leg maximum.

   PRINT, 'Click left mouse button at the MAXIMUM of the left cavity leg.'
   CURSOR, angleg2l, pixleg2l, /DATA, WAIT=3
   WAIT, 0.5

   ;--- Find the scan index corresponding to the intensity maximum
   ;    of the left leg.

   ii = -1
   FOR i = 0, n2-1 DO			$
   BEGIN
      IF (angleg2l GE ts2dx [i]) THEN ii = i
   END

   IF (ii GE 0 AND ii LE n2-1) THEN	$
   BEGIN
      aleg2l = ts2dx   [ii]
      pleg2l = ts2_ave [ii]
      print, 'aleg2l: ', aleg2l
      print, 'pleg2l: ', pleg2l
   END					$
   ELSE					$
   BEGIN
      print, 'Invalid left leg position: ', angleg2l
      print, 'Expected range: ', ts2dx [0], ts2dx [n2-1]
   END

   PLOTS, [angleg2l, angleg2l], [0.9*pleg2l,  1.1*pleg2l],  /DATA, color=ov_red
   PLOTS, [angleg2l-3, angleg2l+3], [pleg2l,  pleg2l],  /DATA, color=ov_red

   ;--- Right leg maximum.

   PRINT, 'Click left mouse button at the MAXIMUM of the right cavity leg.'
   CURSOR, angleg2r, pixleg2r, /DATA, WAIT=3
   WAIT, 0.5

   ;--- Find the scan index corresponding to the intensity maximum
   ;    of the right leg.

   ii = -1
   FOR i = 0, n2-1 DO			$
   BEGIN
      IF (angleg2r GE ts2dx [i]) THEN ii = i
   END

   IF (ii GE 0 AND ii LE n2-1) THEN	$
   BEGIN
      aleg2r = ts2dx   [ii]
      pleg2r = ts2_ave [ii]
      print, 'aleg2r: ', aleg2r
      print, 'pleg2r: ', pleg2r
   END					$
   ELSE					$
   BEGIN
      print, 'Invalid right leg position: ', angleg2r
      print, 'Expected range: ', ts2dx [0], ts2dx [n2-1]
   END

   PLOTS, [angleg2r, angleg2r], [0.9*pleg2r,  1.1*pleg2r],  /DATA, color=ov_red
   PLOTS, [angleg2r-3, angleg2r+3], [pleg2r,  pleg2r],  /DATA, color=ov_red

   PLOTS, [anglow2,  anglow2 ], [0.9*pixlow2, 1.1*pixlow2], /DATA, color=ov_red
   PLOTS, [anglow2 -3, anglow2 +3], [pixlow2, pixlow2], /DATA, color=ov_red

   circle = FINDGEN (16) * (!PI * 2 / 16.0)
   USERSYM, COS (circle), SIN (circle)
;   PLOTS, [angleg2l, angleg2l], [pixleg2l, pixleg2l], /DATA, psym=4, color=fcol
;   PLOTS, [angleg2r, angleg2r], [pixleg2r, pixleg2r], /DATA, psym=4, color=fcol
;   PLOTS, [anglow2,  anglow2 ], [pixlow2,  pixlow2 ], /DATA, psym=8, color=fcol

   ;--- Save azimuthal scan plot [upper height] as a GIF file.

   IF (KEYWORD_SET (gif))  THEN		$
   BEGIN
      gif_file = basename + '.ts2.gif'
      img_gif = TVRD ()
      WRITE_GIF, gif_file, img_gif, red, green, blue
   END

   ;--- Write leg maxima information to the text file.

   CLOSE,  21
   OPENU,  21, pfile, /APPEND
   PRINTF, 21, 'Height: ', r2
   PRINTF, 21, 'Left  leg   angle: ', angleg2l, ' magnitude: ', pixleg2l
   PRINTF, 21, 'Right leg   angle: ', angleg2r, ' magnitude: ', pixleg2r
   CLOSE,  21

   ;----------------------------------------------------------------------------
   ; Plot radial scan.
   ;----------------------------------------------------------------------------

   PLOT, rsdx, rs_ave,						$
	 title='Radial Scan ' + fits_name + ' ' + sangcen + ' deg.',	$
	 xtitle = 'Radius [Rsun]',					$
	 ytitle = ylab,							$
	 color=fcol
   PLOTS, [radmin, radmin], [rs_min, rs_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [radtop, radtop], [rs_min, rs_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [radmin, radtop], [rs_min, rs_min], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [radmin, radtop], [rs_max, rs_max], /DATA, LINESTYLE=2, color=fcol
   XYOUTS, (radtop-radmin)*0.3 + radmin, rs_min * 1.1, 'CAVITY', $
           color=254, charsize=cs, /DATA

   ;--- Save radial scan plot as a GIF file.

   IF (KEYWORD_SET (gif))  THEN		$
   BEGIN
      gif_file = basename + '.rs3.gif'
      img_gif = TVRD ()
      WRITE_GIF, gif_file, img_gif, red, green, blue
   END

;   PLOT, ts11dx, ts11, color=fcol

   !P.BACKGROUND=0
   SET_PLOT, 'X'

   ;----------------------------------------------------------------------------

END