;+
; NAME		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 fcav, fits_name, gif=gif, cm=cm, wmin=wmin, wmax=wmax, text=text, $
	  nolabel=nolabel

   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')
   srsun    = STRING (rsun, FORMAT='(F7.2)')
   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)
   print, 'count_crradius: ', count_crradius, '  crradius: ', crradius
   solar_r  = fxpar (hdu, 'SOLAR_R', count=count_solar_r)

   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 on 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 quadrant to display.

   IF (angle GE 337.5 AND angle LE 360.0)  THEN quad = 1
   IF (angle GE   0.0 AND angle LE  22.5)  THEN quad = 1
   IF (angle GE  22.5 AND angle LT  77.5)  THEN quad = 2
   IF (angle GE  77.5 AND angle LT 112.5)  THEN quad = 3
   IF (angle GE 112.5 AND angle LT 157.5)  THEN quad = 4
   IF (angle GE 157.5 AND angle LT 202.5)  THEN quad = 5
   IF (angle GE 202.5 AND angle LT 247.5)  THEN quad = 6
   IF (angle GE 247.5 AND angle LT 292.5)  THEN quad = 7
   IF (angle GE 292.5 AND angle LT 337.5)  THEN quad = 8

   ixcen = xdim / 2
   iycen = ydim / 2
   qxdim = xdim / 2
   qydim = ydim / 2

   ascanmin =  0.0
   ascanmax = 40.0

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

   CASE quad OF
   1:    BEGIN
	     xmin = ixcen - xdim / 4
	     ymin = iycen
	     qxcen = qxdim / 2 - 0.5
	     qycen = - 0.5
	     ascanmin = 340.0
	     ascanmax =  20.0 + 360.0
         END
   2:    BEGIN
	     xmin = 0
	     ymin = iycen
	     qxcen = qxdim + 0.5
	     qycen = - 0.5
	     ascanmin =  25.0
	     ascanmax =  65.0
         END
   3:    BEGIN
	     xmin = 0
	     ymin = iycen - ydim / 4
	     qxcen = qxdim + 0.5
	     qycen = qydim / 2 - 0.5
	     ascanmin =  70.0
	     ascanmax = 110.0
         END
   4:    BEGIN
	     xmin = 0
	     ymin = 0
	     qxcen = qxdim + 0.5
	     qycen = qydim + 0.5
	     ascanmin = 115.0
	     ascanmax = 155.0
         END
   5:    BEGIN
	     xmin = ixcen - xdim / 4
	     ymin = 0
	     qxcen = qxdim / 2 - 0.5
	     qycen = qydim + 0.5
	     ascanmin = 160.0
	     ascanmax = 200.0
         END
   6:    BEGIN
	     xmin = ixcen
	     ymin = 0
	     qxcen = - 0.5
	     qycen = qydim + 0.5
	     ascanmin = 205.0
	     ascanmax = 245.0
         END
   7:    BEGIN
	     xmin = ixcen
	     ymin = iycen - ydim / 4
	     qxcen = - 0.5
	     qycen = qydim / 2 - 0.5
	     ascanmin = 250.0
	     ascanmax = 290.0
         END
   8:    BEGIN
	     xmin = ixcen
	     ymin = iycen
	     qxcen = - 0.5
	     qycen = - 0.5
	     ascanmin = 295.0
	     ascanmax = 335.0
         END
   ELSE: BEGIN
	     PRINT, 'Invalid Quadrant: ', quad
	     RETURN
	 END
   ENDCASE

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

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

   print, 'xmin, ymin: ', xmin, ymin
   print, 'xmax, ymax: ', xmax, ymax
   print, 'qxcen, qycen: ', qxcen, qycen

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

   sxdim = xdim
   sydim = ydim
   sxcen = qxcen * 2.0
   sycen = qycen * 2.0
   simg = rebin (quadimg, sxdim, sydim)
   spixrs = pixrs * 2.0

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

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

   ;--- Display quadrant image.

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

   ;--- 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 array element is cavity min, cavity max, cavity top.

   cavmin =  1000.0
   cavmax = -1000.0
   cavtop =     1.0
   print, 'numpos: ', numpos
   print, 'rbuf:   ', rbuf
   print, 'abuf:   ', abuf

   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

   ;--- Set scan parameters.

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

   r1 = (rbuf [toppos] - radmin) * 0.33 + radmin 
   r2 = (rbuf [toppos] - radmin) * 0.67 + radmin
   r3 =  rbuf [toppos]

   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, [cavmin, cavmin], [ts1_min, ts1_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [cavmax, cavmax], [ts1_min, ts1_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [cavmin, cavmax], [ts1_min, ts1_min], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [cavmin, cavmax], [ts1_max, ts1_max], /DATA, LINESTYLE=2, color=fcol
   XYOUTS, (cavmax-cavmin)*0.3 + cavmin, 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, [cavmin, cavmin], [ts2_min, ts2_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [cavmax, cavmax], [ts2_min, ts2_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [cavmin, cavmax], [ts2_min, ts2_min], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [cavmin, cavmax], [ts2_max, ts2_max], /DATA, LINESTYLE=2, color=fcol
   XYOUTS, (cavmax-cavmin)*0.3 + cavmin, 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, [cavtop, cavtop], [as_min, as_max], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [radmin, cavtop], [as_min, as_min], /DATA, LINESTYLE=2, color=fcol
   PLOTS, [radmin, cavtop], [as_max, as_max], /DATA, LINESTYLE=2, color=fcol
   XYOUTS, (cavtop-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