;+
; NAME		fc [fcav.pro]
;
; PURPOSE	Display a FITS image & perform cavity measurements.
;
; SYNTAX	fcav, fits_image, /gif, cm='colormap.lut', wmin=1, wmax=240
;
; PARAMETERS
;		fits_image	filename of Mk4 FITS image.
;		gif		write displayed image as a GIF file.
;		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.
;		nolabel		If set, do NOT display the position # label.
;               text    	write position information to a text file.
;
; EXAMPLES	fcav, '19981101.1234.mk4.rpb.fts', $
;		       cm='/home/stanger/color/bwy.lut'
; 		fcav, '19981101.1234.mk4.rpb.fts', /gif
; 		fcav, '19981101.1234.mk4.rpb.fts', /text
; 		fcav, '19981101.1234.mk4.rpb.fts', wmin=0, wmax=2500 
; 		fcav, '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
;		mouse_pos_cav.pro	position/label procedure
;
; 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 fc, fits_name, gif=gif, cm=cm, wmin=wmin, wmax=wmax, text=text, $
	  nolabel=nolabel

   ;--- Set character size.

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

   disp_label = 1			; Set display label option variable.
   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.

   xb = 160 
   yb =   0
   xdim_prev = 0
   ydim_prev = 0

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

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

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

   img = readfits (fits_name, hdu)

   ;--- Extract information from header.

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

   ;--- 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

   ;--- 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')
   type_obs = fxpar (hdu, '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_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') + xb
   ycen     = fxpar (hdu, 'CRPIX2') + yb
   roll     = fxpar (hdu, 'CROTA1')
   revroll  = -roll

;   telescop = STRTRIM (telescop, 2)
;      IF (telescop NE 'SPARTAN 201' AND instrume NE 'WLC')  THEN         $
;      BEGIN
;       PRINT, 'TELESCOP: ', telescop, ' should be "SPARTAN 201".'
;       PRINT, 'INSTRUME: ', instrume, ' should be "WLC".'
;       RETURN
;      END

   ;--- Display 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, xcen, ycen, 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 = middle
   IF (angle GE   0.0 AND angle LE  22.5)  THEN oct = 1	;   0
   IF (angle GE  22.5 AND angle LT  77.5)  THEN oct = 2	;  45
   IF (angle GE  77.5 AND angle LT 112.5)  THEN oct = 3	;  90
   IF (angle GE 112.5 AND angle LT 157.5)  THEN oct = 4	; 135
   IF (angle GE 157.5 AND angle LT 202.5)  THEN oct = 5	; 180
   IF (angle GE 202.5 AND angle LT 247.5)  THEN oct = 6	; 225
   IF (angle GE 247.5 AND angle LT 292.5)  THEN oct = 7	; 270
   IF (angle GE 292.5 AND angle LT 337.5)  THEN oct = 8	; 315

   print, 'oct: ', oct
   ixcen = xdim / 2
   iycen = ydim / 2
   oxdim = xdim / 2
   oydim = ydim / 2

   ascanmin =  0.0
   ascanmax = 40.0

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

   CASE oct OF
   1:    BEGIN	;   0 degrees octant
	     xmin = ixcen - xdim / 4
	     ymin = iycen
	     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 = iycen
	     oxcen = oxdim + 0.5
	     oycen = - 0.5
	     ascanmin =  25.0
	     ascanmax =  65.0
         END
   3:    BEGIN	;  90 degrees octant
	     xmin = 0
	     ymin = iycen - 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 = ixcen - 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 = ixcen
	     ymin = 0
	     oxcen = - 0.5
	     oycen = oydim + 0.5
	     ascanmin = 205.0
	     ascanmax = 245.0
         END
   7:    BEGIN	; 270 degrees octant
	     xmin = ixcen
	     ymin = iycen - ydim / 4
	     oxcen = - 0.5
	     oycen = oydim / 2 - 0.5
	     ascanmin = 250.0
	     ascanmax = 290.0
         END
   8:    BEGIN	; 335 degrees octant
	     xmin = ixcen
	     ymin = iycen
	     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

   ;--- 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, rxcur, rycur, /DEVICE, WAIT=3
   ierror = rcoord (radius, angle, rxcur, rycur, -1, roll, sxcen, sycen, spixrs)

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

   radtop = radius
   angtop  = angle
   radmin  = 1.1
   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.

   IF (KEYWORD_SET (nolabel)) THEN					$
      cavity_width, sxdim, sydim, sxcen, sycen, spixrs, revroll,	$
		    rbuf, abuf, ascanmin, ascanmax,			$
		    radtop, angtop, r1, r2,				$
		    numpos, pos=1, pfile=pfile				$
   ELSE									$
      cavity_width, sxdim, sydim, sxcen, sycen, spixrs, revroll,	$
		    rbuf, abuf, ascanmin, ascanmax,			$
		    radtop, angtop, r1, r2,				$
		    numpos, pos=1, pfile=pfile,				$
		    disp_label=disp_label

   ;--- Do position measurements.

;   IF (KEYWORD_SET (nolabel)) THEN					$
;   BEGIN
;     mouse_pos_cav, sxdim, sydim, sxcen, sycen, spixrs, revroll,	$
;		    rbuf, abuf, numpos,			         	$
;		    pos=1, pfile=pfile
;   ENDIF								$
;   ELSE								$
;   BEGIN
;     print, 'mouse_pos_cav'
;     mouse_pos_cav, sxdim, sydim, sxcen, sycen, spixrs, revroll,	$
;		    rbuf, abuf, numpos,			         	$
;		    pos=1, pfile=pfile, /disp_label
;   ENDELSE

   ;--- Write displayed 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

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

   print, 'numpos: ', numpos
   print, 'rbuf:   ', rbuf
   print, 'abuf:   ', abuf

;   cavmin =    360.0
;   cavmax =      0.0
;   cavtop =      0.0
;
;   FOR i = 0, numpos - 1  DO	$
;   BEGIN
;      print, 'i: ', i, ' rbuf: ', rbuf [i], ' abuf: ', abuf [i]
;      IF (abuf [i] LT cavmin) THEN $
;	 BEGIN
;	    cavmin = abuf [i]
;	    minpos = i
;         END
;      IF (abuf [i] GT cavmax) THEN $
;	 BEGIN
;	    cavmax = abuf [i]
;	    maxpos = i
;         END
;      IF (rbuf [i] GT cavtop) THEN $
;	 BEGIN
;	    cavtop = rbuf [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.

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

   FOR i = 0, numpos - 1 DO	$
   BEGIN
      IF (ss [i] EQ 1 AND abuf [i] GT amax1) THEN	$
      BEGIN
         amax1 = abuf [i]
	 rmax1 = rbuf [i]
      END

      IF (ss [i] EQ 1 AND abuf [i] LT amin1) THEN	$
      BEGIN
         amin1 = abuf [i]
	 rmin1 = rbuf [i]
      END

      IF (ss [i] EQ 2 AND abuf [i] GT amax2) THEN	$
      BEGIN
         amax2 = abuf [i]
	 rmax2 = rbuf [i]
      END

      IF (ss [i] EQ 2 AND abuf [i] LT amin2) THEN	$
      BEGIN
         amin2 = abuf [i]
	 rmin2 = rbuf [i]
      END

   END

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

   angcen = (amax2 - amax1) * 0.5 + amax1

   ;--- 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)

;   aa1 = cavmin 
;   aa2 = cavmax
;   IF (aa2 LT aa1) THEN aa2 = aa2 + 360.0

;   print, 'aa1/aa2: ', aa1, aa2

;   angcen  = (aa2 - aa1) / 2.0 + aa1
   sangcen = STRTRIM (STRING (angcen, FORMAT='(F7.2)'), 2)

;   print, 'angcen: ', angcen

   ;--- Do azimuthal scan at first height.

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

   ;--- Smooth scan data.

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

   ;--- Do azimuthal scan at second 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)

   ;--- Do radial scan at center of cavity.

   rscan, fits_name, img, pixrs, revroll, xcen, ycen,	$
	  radbeg, radend, radinc, angcen,	$
	  as,  asdx,  n3

   ;--- Smooth scan data.

   as_ave = SMOOTH (as, 3)
   as_min = MIN (as_ave)
   as_max = MAX (as_ave)

   WDELETE, 2

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

   SET_PLOT, 'Z'
   LOADCT, 0
   fcol =   0
   bcol = 255
   !P.BACKGROUND=bcol
   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

   ;--- Save plot as a GIF file.

   gif_file = basename + '.ts1.gif'
   img_gif = TVRD ()
   WRITE_GIF, gif_file, img_gif, red, green, blue

   ;--- Plot azimuthal scan at second 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

   ;--- Save plot as a GIF file.

   gif_file = basename + '.ts2.gif'
   img_gif = TVRD ()
   WRITE_GIF, gif_file, img_gif, red, green, blue

   ;--- Plot radial scan.

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

   ;--- Save plot as a GIF file.

   gif_file = basename + '.rs3.gif'
   img_gif = TVRD ()
   WRITE_GIF, gif_file, img_gif, red, green, blue

   !P.BACKGROUND=0
   SET_PLOT, 'X'

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

END