;+
; fms.pro   IDL procedure
;
; PURPOSE
;	Compute absolute/excess coronal mass.
;	Read FITS images and define a region of interest.
;
; SYNTAX
;	fms, 'fits_list'
;
; PROCEDURES USED
;	readfits	Read FITS file
;	fxpar		Get FITS parameter associated with a keyword
;	write_gif	Write GIF file
;	defroic		(tailored version of "defroi")
;	dte		LASCO use only
;	get_cal_name	LASCO use only
;       fms_cyl		curved wedge.
;	kb_cyl		Keyboard input to define curved wedge
;	mouse_cyl	Mouse cursor   to define curved wedge
;	kb_wedge	Keyboard input to define straight wedge
;	mouse_wedge	Mouse cursor   to define straight wedge
;	mouse_poly	Mouse cursor   to define polygon
;	file_poly	File           to define polygon
;	xmass		compute excess/absolute mass
;
; HISTORY
;	Author: Andrew L. Stanger   2 June 1998
;	 1 Oct 1998: Add logic for C1, C2, C3, EIT telescopes.
;	19 Oct 1998: Add scope_name argument to xmass function.
;	10 Nov 1998: Change 'defroi' name to 'defroic' to distinguish
;		     it from the RSI version.
;	28 Jan 1999: version 3: Add ability to store/fetch polygon vertex file.
;	29 Jan 1999: version 4: Add ability to define a wedge area.
;	 8 Feb 1999: version 5: Add ability to define a cylindrical wedge area.
;	11 Feb 1999: Use FINDFILE to avoid overwriting *.pv, *.gif, *.fm files.
;	12 Feb 1999: Make cmeyear 4 digits, not 2 (i.e., Y2K compliance).
;	14 May 1999: Accomodate "mass" SMM C/P images.
;	20 May 1999: version 0.7: use crota1 & scroll to establish solar north. 
;	23 Nov 1999: version 0.8: do NOT use threshold in computing mass.
;	 5 Feb 2001: version 0.9: Add Mk4.
;	 5 Dec 2006: version 1.0: Absolute mass only.  Process list of images.
;-

PRO fms, fits_list

cme_fits_file = ''
ref_fits_file = ''
index = 0
floc  = LONG (0)
fpos  = LONARR (1200)

;--- Image list loop.

index = 0
OPENR, 11, fits_list
WHILE (NOT EOF (11)) DO $
BEGIN  ;{
   POINT_LUN, -11, floc
   fpos [index] = floc

   ;--- Read FITS name from file list.

   READF, 11, cme_fits_file
   PRINT, 'FITS image: ', cme_fits_file

   ;--- Initialize variables.

   i      = 0L
   nloc   = 0L
   regnum = 0

   red    = bytarr (256)
   green  = bytarr (256)
   blue   = bytarr (256)

   imgcme = 0
   imgref = 0

   mx     = 0			; =0 cme image only, =1 cme & ref images
   tscope = 0			; Telescope designation
				; 1: C1   2: C2   3: C3   4: EIT   5: SXT
				; 6: MK3  7: Mk4  8: SMM C/P

   xdim_cme     = 0		; X-axis dimension
   ydim_cme     = 0		; Y-axis dimension
   nbypix_cme   = 0		; # bytes/pixel
   xcen_cme     = 0		; X-axis center
   ycen_cme     = 0		; Y-axis center
   crradius_cme = 0.0		; # pixels/Rsun
   hdu_cme      = ''		; CME header
   telescop_cme = ''		; Telescope
   instrume_cme = ''		; Instrument
   detector_cme = ''		; Detector
   roll_cme     = 0.0		; Image "roll".
   scroll_cme   = 0.0		; S/C roll.

   xdim_ref     = 0		; X-axis dimension
   ydim_ref     = 0		; Y-axis dimension
   nbypix_ref   = 0		; # bytes/pixel
   xcen_ref     = 0		; X-axis center
   ycen_ref     = 0		; Y-axis center
   crradius_ref = 0.0		; # pixels/Rsun
   hdu_ref      = ''		; REF header
   telescop_ref = ''		; Telescope
   instrume_ref = ''		; Instrument
   detector_ref = ''		; Detector
   vig_fn       = ''

   ;--- Read CME FITS image.

;   imgcme = rfits (cme_fits_file, header=hdu_cme, /scale)

   imgcme = readfits (cme_fits_file, hdu_cme)
   dspcme = imgcme

   xdim_cme     = fxpar (hdu_cme, "NAXIS1")
   ydim_cme     = fxpar (hdu_cme, "NAXIS2")
   nbypix_cme   = fxpar (hdu_cme, "BITPIX") * 8
   crpix1_cme   = fxpar (hdu_cme, "CRPIX1",   count=p_crpix1)
   crpix2_cme   = fxpar (hdu_cme, "CRPIX2",   count=p_crpix2)
   xcen_cme     = fxpar (hdu_cme, "XCEN",     count=p_xcen)
   ycen_cme     = fxpar (hdu_cme, "YCEN",     count=p_ycen)
   cdelt1_cme   = fxpar (hdu_cme, "CDELT1",   count=p_cdelt1)
   cdelt2_cme   = fxpar (hdu_cme, "CDELT1",   count=p_cdelt1)

   print, 'fm> xcen_cme   : ', xcen_cme
   print, 'fm> ycen_cme   : ', ycen_cme
   print, 'fm> cdelt1_cme : ', cdelt1_cme
   print, 'fm> cdelt2_cme : ', cdelt2_cme

   IF (xcen_cme GT xdim_cme) THEN	$
   BEGIN
       xcen_cme = xcen_cme / cdelt1_cme
       ycen_cme = ycen_cme / cdelt2_cme
   END					$
   ELSE					$
   BEGIN
       xcen_cme = crpix1_cme
       ycen_cme = crpix2_cme
   END
   print, 'fm> xcen_cme: ', xcen_cme
   print, 'fm> ycen_cme: ', ycen_cme

   crradius_cme = fxpar (hdu_cme, "CRRADIUS", count=p_crradius)
   pixrs 	= crradius_cme
   rsun         = fxpar (hdu_cme, "RSUN", count=p_rsun)
   crota1_cme   = fxpar (hdu_cme, "CROTA1")
   scroll_cme   = fxpar (hdu_cme, "SCROLL")
   roll_cme     = -crota1_cme
   print, 'fm> crradius_cme: ', crradius_cme
   print, 'fm> cdelt1/2    : ', cdelt1_cme, cdelt2_cme
   print, 'fm> rsun        : ', rsun
   IF (p_crradius LT 1) THEN		$
   BEGIN
       pixrs = rsun / cdelt1_cme
   END
   print, 'fm> pixrs       :', pixrs

   bunit_cme    = fxpar (hdu_cme, "BUNIT")
   dispmin_cme  = fxpar (hdu_cme, "DISPMIN")
   dispmax_cme  = fxpar (hdu_cme, "DISPMAX")
;   print, 'CME dispmin/max: ', dispmin_cme, dispmax_cme

   telescop_cme = fxpar (hdu_cme, "TELESCOP")
   instrume_cme = fxpar (hdu_cme, "INSTRUME")
   detector_cme = fxpar (hdu_cme, "DETECTOR")
   cmedate   = fxpar (hdu_cme, "DATE-OBS")
   cmetime   = fxpar (hdu_cme, "TIME-OBS")
   print, 'cmedate: ', cmedate
   print, 'cmetime: ', cmetime

;--- Extract day, month & year from 'DATE-OBS'.

;   cmeday    = STRMID (cmedate, 0, 2)
;   cmemonth  = STRMID (cmedate, 3, 2)
;   cmeyear   = STRMID (cmedate, 6, STRLEN (cmedate) - 6)

   cmeyear   = STRMID (cmedate, 0, 4)
   cmemonth  = STRMID (cmedate, 5, 2)
   cmeday    = STRMID (cmedate, 8, 2)
   print, 'STRLEN (cmedate): ', STRLEN (cmedate)
   print, 'cmeyear:...', cmeyear, '...'

   IF (STRLEN (cmeyear) EQ 2) THEN			$
   BEGIN
      IF (cmeyear GE 50) THEN cmeyear = cmeyear + 1900	$
      ELSE                    cmeyear = cmeyear + 2000
   END
   print, 'STRLEN (cmeyear): ', STRLEN (cmeyear)
   IF (STRLEN (cmeyear) GT 4)  THEN $
      cmeyear = STRCOMPRESS (STRTRIM (STRMID (cmeyear, STRLEN (cmeyear) - 4, 4)))
   print, 'cmeyear:...', cmeyear, '...'
   yymmdd = cmeyear + cmemonth + cmeday
   print, 'fm> yymmdd: ', yymmdd

;--- LASCO/C1/C2/C3:

   dte = STR2UTC (cmedate)
   yymmdd = UTC2YYMMDD (dte)
;   cmeyear   = STRMID (yymmdd, 0, 2)
;   cmemonth  = STRMID (yymmdd, 2, 2)
;   cmeday    = STRMID (yymmdd, 4, 2)

   cmehour   = STRMID (cmetime, 0, 2)
   cmeminute = STRMID (cmetime, 3, 2)
   cmesecond = STRMID (cmetime, 6, 2)
   PRINT, 'cme year, month, day: ', cmeyear, cmemonth, cmeday
   PRINT, 'cme hour, minute, second: ', cmehour, cmeminute, cmesecond

   IF (STRPOS (telescop_cme, "SMM C/P") NE -1)  THEN	$	; SMM C/P image.
      BEGIN
      tscope = 8
      dspcme = imgcme ;
      END

   IF (STRPOS (telescop_cme, "MK4") NE -1   OR		$
       STRPOS (telescop_cme, "Mk4") NE -1)  THEN	$	; Mk4 image.
   BEGIN	;{
      tscope = 7
      imgvig = FLTARR (960, 960)
      OPENR, 1, '/home/stanger/vignet/mk4vig.template'
      READU, 1, imgvig
      CLOSE, 1
;      dspcme = imgcme * imgvig			; Apply "vignetting"
      PRINT, 'Apply MK4 vignetting to dspcme.'
   END		;}

   IF (STRPOS (telescop_cme, "MK3") NE -1   OR		$
       STRPOS (telescop_cme, "Mk3") NE -1)  THEN	$	; Mk3 image.
   BEGIN	;{
      tscope = 6
      imgvig = FLTARR (512, 512)
      OPENR, 1, '/home/stanger/vignet/mk3vig.template'
      READU, 1, imgvig
      CLOSE, 1
      dspcme = imgcme * imgvig				; Apply "vignetting"
      PRINT, 'Apply MK3 vignetting to dspcme.'
   END		;}

   IF (STRPOS (telescop_cme, "SOHO") NE -1)  THEN		$	; SOHO
   BEGIN	;{
      IF (STRPOS (instrume_cme, "LASCO") NE -1) THEN		$	; LASCO
      BEGIN		;{

	 IF (STRPOS (detector_cme, "C1") NE -1)  THEN		$	; C1
	 BEGIN	;{
	    tscope = 1
	 END	$	;}

         ELSE IF (STRPOS (detector_cme, "C2") NE -1) THEN	$	; C2
	 BEGIN	;{
	    tscope = 2
;	    sd = GETENV ('NRL_LIB') + '/lasco/data/calib/'
	    sd = GETENV ('SSW') + '/soho/lasco/lasco/data/calib/'
	    sd = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/'
	    sv = sd + 'C2_cl*vig*.dat'
	    print, 'sv: ', sd 
	    dte = STR2UTC (cmedate)
	    yymmdd = UTC2YYMMDD (dte)
	    print, 'yymmdd: ', yymmdd
            dispmin_cme  = fxpar (hdu_cme, "DATAMIN")
	    dispmin_cme  = 0.0
            dispmax_cme  = fxpar (hdu_cme, "DATAMAX") / 10.0
	    vig_fn = ''
;	    vig_fn = get_cal_name (sd + 'C2_cl*vig*.dat', yymmdd)
            vig_fn = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/C2_cl_vig_960611_v01.dat'
	    IF (vig_fn NE '') THEN	$
	    BEGIN		;{
	       dte_vig = FILE_DATE_MOD (vig_fn, /DATE_ONLY)
	       print, 'vig_fn: ', vig_fn
	       vig = READFITS (vig_fn, hvig)
	       zv = WHERE (vig EQ 0.0)
	       vig [zv] = 1.0
	       dspcme = imgcme / vig 
	    ENDIF	$	;}
	    ELSE BEGIN		;{
	       dspcme = imgcme
	       print,'ERROR: c2_calibrate - No ' + sd + 'C2_cl*vig*.dat file'
;	       RETURN
	    END		;}
         END	$	;}

         ELSE IF (STRPOS (detector_cme, "C3") NE -1) THEN	$	; C3
         BEGIN	;{
            tscope = 3
	    sd = GETENV ('NRL_LIB') + '/lasco/data/calib/'
	    sd = GETENV ('SSW') + '/soho/lasco/lasco/data/calib/'
	    sd = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/'
	    dte = STR2UTC (cmedate)
	    yymmdd = UTC2YYMMDD (dte)
	    print, 'fm> dte: ', dte, '  yymmdd: ', yymmdd
	    vig_fn = ''
;	    vig_fn = get_cal_name (sd + 'C3_cl*vig*.dat', yymmdd)
            vig_fn = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/C3_cl_vig_951201_v03.dat'
	    IF (vig_fn NE '') THEN		$
	    BEGIN	;{
	       dte_vig = FILE_DATE_MOD (vig_fn, /DATE_ONLY)
	       print, 'vig_fn: ', vig_fn
	       vig = READFITS (vig_fn, hvig)
	       dspcme = imgcme / vig 
	    ENDIF	$;}
	    ELSE BEGIN	;{
	       dspcme = imgcme
	       print, 'ERROR: c3_calibrate - No ' + sd + 'C3_cl*vig*.dat file'
;	       RETURN
            END			;}
         END		;}

      END	$	;}

      ELSE IF (STRPOS (instrume_cme, "EIT") NE -1) THEN		$	; EIT
      BEGIN		;{
	 tscope = 4
      END		;}

   END	;}

   IF (STRPOS (telescop_cme, "SXT") NE -1)  THEN		$	; SXT
      BEGIN
      tscope = 5
      END

;--- Use NRL_LIB version of pb0r:

   date_time = cmedate + ' ' + cmetime
   solar_info = pb0r (date_time)		; Get solar P, B0, Rsun.
   Pangle  = solar_info [0]			; P-angle
   B0      = solar_info [1]			; B0 angle
   Rsunarc = solar_info [2] * 60.0		; Solar radius (arcsec).
;   pixrs = Rsunarc / cdelt1_cme			; Pixels / Rsun.
;   print, 'fm> pixrs: ', pixrs

   WINDOW, xsize=xdim_cme, ysize=ydim_cme, retain=2	; Create image window.

   ;--- Reference Image ?

   IF (N_PARAMS() EQ 2) THEN	$	; Reference (pre-event) image supplied.
      BEGIN		;{
;      imgref = rfits (ref_fits_file, header=hdu_ref, /scale)
      imgref = readfits (ref_fits_file, hdu_ref)
      dspref = imgref

      xdim_ref     = fxpar (hdu_ref, "NAXIS1")
      ydim_ref     = fxpar (hdu_ref, "NAXIS2")
      nbypix_ref   = fxpar (hdu_ref, "BITPIX") * 8
      crpix1_ref   = fxpar (hdu_ref, "CRPIX1",   count=p_crpix1)
      crpix2_ref   = fxpar (hdu_ref, "CRPIX2",   count=p_crpix2)
      xcen_ref     = fxpar (hdu_ref, "XCEN")
      ycen_ref     = fxpar (hdu_ref, "YCEN")
      cdelt1_ref   = fxpar (hdu_ref, "CDELT1",   count=p_cdelt1)
      cdelt2_ref   = fxpar (hdu_ref, "CDELT1",   count=p_cdelt1)

      print, 'fm> xdim_ref   : ', xdim_ref
      print, 'fm> ydim_ref   : ', ydim_ref
      print, 'fm> xcen_ref   : ', xcen_ref
      print, 'fm> ycen_ref   : ', ycen_ref
      print, 'fm> cdelt1_ref : ', cdelt1_ref
      print, 'fm> cdelt2_ref : ', cdelt2_ref

      IF (xcen_ref GT xdim_ref) THEN	$
      BEGIN
         xcen_ref = xcen_ref / cdelt1_ref
         ycen_ref = ycen_ref / cdelt2_ref
      END					$
      ELSE					$
      BEGIN
         xcen_ref = crpix1_ref
         ycen_ref = crpix2_ref
      END
      print, 'fm> xcen_ref: ', xcen_ref
      print, 'fm> ycen_ref: ', ycen_ref

      crradius_ref = fxpar (hdu_ref, "CRRADIUS")
      telescop_ref = fxpar (hdu_ref, "TELESCOP")
      instrume_ref = fxpar (hdu_ref, "INSTRUME")
      detector_ref = fxpar (hdu_ref, "DETECTOR")
      refdate      = fxpar (hdu_ref, "DATE-OBS")
      reftime      = fxpar (hdu_ref, "TIME-OBS")

      bunit_ref    = fxpar (hdu_ref, "BUNIT")
;      refday    = STRMID (refdate, 0, 2)
;      refmonth  = STRMID (refdate, 3, 2)
;      refyear   = STRMID (refdate, 6, STRLEN (refdate) - 6)

      refyear   = STRMID (refdate, 0, 4)
      refmonth  = STRMID (refdate, 5, 2)
      refday    = STRMID (refdate, 8, 2)
      IF (STRLEN (refyear) EQ 2) THEN			$
      BEGIN
         IF (refyear GE 50) THEN refyear = refyear + 1900	$
         ELSE                    refyear = refyear + 2000
      END
      yymmdd = refyear + refmonth + refday

;--- LASCO/C1/C2/C3:

;      dte = STR2UTC (refdate)
;      yymmdd = UTC2YYMMDD (dte)
;      refyear   = STRMID (yymmdd, 0, 2)
;      refmonth  = STRMID (yymmdd, 2, 2)
;      refday    = STRMID (yymmdd, 4, 2)

      refhour   = STRMID (reftime, 0, 2)
      refminute = STRMID (reftime, 3, 2)
      refsecond = STRMID (reftime, 6, 2)

;      base_name = STRTRIM (STRING (cmeyear)			$
;                + STRING (cmemonth) + STRING (cmeday)		$
;		+ '.'						$
;		+ STRING (cmehour) + STRING (cmeminute) + '-'	$
;                + STRING (refhour) + STRING (refminute)
      base_name = cmeyear + cmemonth  + cmeday + '.' +	$
	          cmehour + cmeminute + '-'    +	$
	          refhour + refminute
      err = 0
      IF (xdim_cme NE xdim_ref)     THEN err = -1
      IF (ydim_cme NE ydim_ref)     THEN err = -2
      IF (nbypix_cme NE nbypix_ref) THEN err = -3
      IF (xcen_cme NE xcen_ref)     THEN err = -4
      IF (ycen_cme NE ycen_ref)     THEN err = -5
      print, 'fm> xcen_cme/ycen_cme: ', xcen_cme, ycen_cme
      print, 'fm> xcen_ref/ycen_ref: ', xcen_ref, ycen_ref
      print, 'fm> err: ', err
      IF (err NE 0) THEN RETURN

      if (STRPOS (telescop_ref, "SMM C/P") NE -1)  THEN	$	; SMM C/P image.
         BEGIN
	 tscope = 8
         dspref = imgref
         END

      if (STRPOS (telescop_ref, "MK4") NE -1)  THEN	$	; Mk4 image.
         BEGIN
	 tscope = 7
         imgvig = fltarr (960, 960)
         openr, 1, '/home/stanger/vignet/mk4vig.template'
         readu, 1, imgvig
         close, 1
;         dspref = imgref * imgvig			; Apply "vignetting"
         print, 'fm> Apply MK4 vignetting to dspref.'
         END

      if (STRPOS (telescop_ref, "MK3") NE -1)  THEN	$	; MK3 image.
         BEGIN
	 tscope = 6
         imgvig = fltarr (512, 512)
         openr, 1, '/home/stanger/vignet/mk3vig.template'
         readu, 1, imgvig
         close, 1
         dspref = imgref * imgvig			; Apply "vignetting"
         print, 'fm> Apply MK3 vignetting to dspref.'
         END

      IF (STRPOS (telescop_ref, "SOHO") NE -1)  THEN		$	; SOHO
      BEGIN	;{
         IF (STRPOS (instrume_ref, "LASCO") NE -1) THEN		$	; LASCO
         BEGIN		;{

	    IF (STRPOS (detector_ref, "C1") NE -1)  THEN	$	; C1
	    BEGIN	;{
	       tscope = 1
	    END	$	;}

            ELSE IF (STRPOS (detector_ref, "C2") NE -1) THEN	$	; C2
	    BEGIN	;{
	       tscope = 2
	       sd = GETENV ('NRL_LIB') + '/lasco/data/calib/'
	       sd = GETENV ('SSW') + '/soho/lasco/lasco/data/calib/'
	       sd = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/'
	       dte = STR2UTC (cmedate)
	       yymmdd = UTC2YYMMDD (dte)
	       vig_fn = ''
;	       vig_fn = get_cal_name (sd + 'C2_cl*vig*.dat', yymmdd)
               vig_fn = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/C2_cl_vig_960611_v01.dat'
	       IF (vig_fn NE '') THEN	$
	       BEGIN		;{
	          dte_vig = FILE_DATE_MOD (vig_fn, /DATE_ONLY)
	          print, 'fm> vig_fn: ', vig_fn
	          vig = READFITS (vig_fn, hvig)
	          zv = WHERE (vig EQ 0.0)
	          vig [zv] = 1.0
	          dspref = imgref / vig 
	       ENDIF	$	;}
	       ELSE BEGIN		;{
		  dspref = imgref
	          print,'fm> ERROR: c2_calibrate - No ' + sd + 'C2_cl*vig*.dat file'
;	          RETURN
	       END		;}
            END	$	;}

            ELSE IF (STRPOS (detector_ref, "C3") NE -1) THEN	$	; C3
            BEGIN	;{
               tscope = 3
	       sd = GETENV ('NRL_LIB') + '/lasco/data/calib/'
	       sd = GETENV ('SSW') + '/soho/lasco/lasco/data/calib/'
	       sd = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/'
	       dte = STR2UTC (cmedate)
	       yymmdd = UTC2YYMMDD (dte)
	       print, 'fm> dte: ', dte, '  yymmdd: ', yymmdd
	       vig_fn = ''
;	       vig_fn = get_cal_name (sd + 'C3_cl*vig*.dat', yymmdd)
               vig_fn = '/project/cordyn/sw/ssw/soho/lasco/idl/data/calib/C3_cl_vig_951201_v03.dat'
	       IF (vig_fn NE '') THEN		$
	       BEGIN	;{
	          dte_vig = FILE_DATE_MOD (vig_fn, /DATE_ONLY)
	          print, 'fm> vig_fn: ', vig_fn
	          vig = READFITS (vig_fn, hvig)
	          zv = WHERE (vig EQ 0.0)
	          vig [zv] = 1.0
	          dspref = imgref / vig 
	       ENDIF	$;}
	       ELSE BEGIN	;{
		  dspref = imgref
	          print, 'fm> ERROR: c3_calibrate - No ' + sd + 'C3_cl*vig*.dat file'
;	          RETURN
               END			;}
            END		;}

         END	$	;}

         ELSE IF (STRPOS (instrume_ref, "EIT") NE -1) THEN		$	; EIT
         BEGIN		;{
	    tscope = 4
         END		;}
      END	;}

      IF (STRPOS (telescop_ref, "SXT") NE -1)  THEN		$	; SXT
         BEGIN
         tscope = 5
         END

      mx = 1
      lct, '/home/stanger/color/dif.lut'	; Load difference color map.
      dspdif = dspcme - dspref
      imin = MIN (dspdif, max=imax)
      print, 'fm> dif imin/imax: ', imin, imax

      wmin = imin
      wmax = imax

      IF (tscope EQ 2) THEN	$
      BEGIN
	 wmin = imin * 0.2
         wmax = imax * 0.2
      END

;      IF (tscope EQ 8) THEN wmax = imax * 0.5
      IF (tscope EQ 7) THEN wmax = imax * 0.2
      IF (tscope EQ 6) THEN wmax = imax * 0.2

      IF (ABS (wmin) LT abs (wmax))  THEN wmax = - wmin
      IF (ABS (wmin) GT abs (wmax))  THEN wmin = - wmax
      print, 'fm> dif wmin/wmax: ', wmin, wmax

      dspdif = BYTSCL (dspdif, min=wmin, max=wmax, top=249)
      TV, dspdif				; Display difference image.
      END	$	;}

   ;--- No Reference image.

   ELSE		$		; No reference image supplied.
      BEGIN	;{
      mx = 0
      imgref = 0
      lct, '/home/stanger/color/bwy.lut'   ; Load blue-white-yellow color map.
      lct, '/home/stanger/color/bwcp.lut'  ; Load black/white       color map.
      print, 'fm> size (dspcme): ', size (dspcme)
      imin = MIN (dspcme, max=imax)
      print, 'fm> CME min/max range: ', imin, imax
      IF (imin LT dispmin_cme) THEN imin = dispmin_cme
      IF (imax GT dispmax_cme) THEN imax = dispmax_cme
      print, 'CME min/max range: ', imin, imax

      print, 'fm> CME disp min/max : ', dispmin_cme, dispmax_cme
      IF (tscope EQ 6 OR tscope EQ 7)  THEN imin = imin * 0.5 ;
      IF (tscope EQ 6 OR tscope EQ 7)  THEN imax = imax * 0.75 ;
      print, 'fm> CME disp min/max : ', dispmin_cme, dispmax_cme

      IF (dispmin_cme NE dispmax_cme)  THEN				$
	 BEGIN
         dspcmeb = bytscl (dspcme, min=imin, max=imax, top=249)
	 END								$
      ELSE								$
	 dspcmeb = bytscl (dspcme)

      IF (STRPOS (detector_cme, "C2") NE -1) THEN 			$
      BEGIN
         wmin = imax / 200.0
         wmax = imax /  40.0
         print, 'wmin/wmax: ', wmin, wmax
         dspcmeb = bytscl (dspcme, min=wmin, max=wmax, top=249)
      END

      TV, dspcmeb			   ; Display image.
;      TVSCL, dspcme

;      base_name  = STRTRIM (STRING (cmeyear))			$
;                 + STRING (cmemonth) + STRING (cmeday) + '.'	$
;                 + STRING (cmehour) + STRING (cmeminute)

      base_name  = cmeyear + cmemonth  + cmeday + '.' +	cmehour + cmeminute 
      END	;}

   CASE tscope OF
      1: scope_name = 'soho_c1'
      2: scope_name = 'soho_c2'
      3: scope_name = 'soho_c3'
      4: scope_name = 'soho_eit'
      5: scope_name = 'yohkoh_sxt'
      6: scope_name = 'mlso_mk3'
      7: scope_name = 'mlso_mk4'
      8: scope_name = 'smm_cp'
      ELSE: scope_name = 'Unknown'
   ENDCASE

   IF (mx EQ 0) THEN TV, dspcmeb	$	; Display CME image.
   ELSE TV, dspdif				; Display difference image.

   regnum = regnum + 1			; Increment region number.

   ;--- Determine valid region number & define output file names.

outfiles: IF (regnum LE 9) THEN	$
   BEGIN
   pv_file  = scope_name + '.' + base_name + '.r0' + $
              STRTRIM (STRING (regnum, FORMAT='(I1)')) + '.pv'
   aoi_file = scope_name + '.' + base_name + '.r0' + $
              STRTRIM (STRING (regnum, FORMAT='(I1)')) + '.gif'
   log_file = scope_name + '.' + base_name + '.r0' + $
              STRTRIM (STRING (regnum, FORMAT='(I1)')) + '.fm'
   END				$
   ELSE				$
   BEGIN
   pv_file  = scope_name + '.' + base_name + '.r' + $
              STRTRIM (STRING (regnum, FORMAT='(I2)')) + '.pv'
   aoi_file = scope_name + '.' + base_name + '.r' + $
              STRTRIM (STRING (regnum, FORMAT='(I2)')) + '.gif'
   log_file = scope_name + '.' + base_name + '.r' + $
              STRTRIM (STRING (regnum, FORMAT='(I2)')) + '.fm'
   END

   ;--- Check for pre-existing mass log file.
   ;    If it exists, increment the region number and try again.

   fname = FINDFILE (log_file, COUNT=nfiles)
   if (nfiles > 0)  THEN	$
   BEGIN
   regnum = regnum + 1
   GOTO, outfiles
   END

   print, 'fm>  pv_file: ',  pv_file
   print, 'fm> aoi_file: ', aoi_file
   print, 'fm> log_file: ', log_file

   ;--- Choose an annulus area 1.2 - 2.2 Rsun, 0 - 359.5 degrees.

   rmin   =   1.2
   rmax   =   2.2
   angmin =   0.0
   angmax = 360.0

   fms_cyl, rmin, rmax, angmin, angmax,					$
	    xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme,	$
	    regnum, log_file, pv_file, xv, yv, nver, polypix

   ;--- Determine coordinates of all pixel within polygon area.

   size_roi = SIZE (polypix)

   nloc = size_roi [1]
   xloc = intarr (nloc)
   yloc = intarr (nloc)

   print, 'nloc: ', nloc

   ;--- Convert 1-D coordinates to 2-D coordinates.

   FOR i = 0L, nloc-1 DO		$
      BEGIN
	 yloc [i] = polypix (i) / xdim_cme
	 xloc [i] = polypix (i) - yloc (i) * xdim_cme
      END

   size_xv = SIZE (xv)
   nver = size_xv [1]

;   print, 'size [xv]: ', size (xv)
;   print, 'size [yv]: ', size (yv)

;   FOR i = 0L, nver - 1 DO		$
;      print, 'i, xv, yv: ', i, xv [i], yv [i]

;   print, 'size (polypix): ', size (polypix)
;   FOR i = 0, nloc - 1 DO	$
;      print, 'polypix, xloc, yloc: ', polypix (i), xloc (i), yloc (i)

   ;--- Compute integrated absolute/excess mass inside polygon.

   xm = xmass10 (log_file, scope_name, ref_fits_file, cme_fits_file,	$
	       imgref, imgcme, bunit_cme, bunit_ref,			$
	       xdim_cme, ydim_cme, nbypix, xcen_cme, ycen_cme,		$
	       pixrs, rsun, roll_cme,					$
	       xloc, yloc, nloc, tscope)

   print, 'xm: ', xm

   ;--- Annotate image.

   xtop = xdim_cme - 1
   ytop = ydim_cme - 1
   XYOUTS,   2, xtop-10, cme_fits_file, /device, color=254, charsize=0.6
   IF (mx EQ 0) THEN	$
      XYOUTS,   2, xtop-25, telescop_cme,  /device, color=254, charsize=0.6 $
   ELSE			$
   BEGIN
      XYOUTS,   2, xtop-25, ref_fits_file, /device, color=254, charsize=0.6
      XYOUTS,   2, xtop-40, telescop_cme,  /device, color=254, charsize=0.6	
   END

   XYOUTS, 420, xtop-10, cmedate,  /device, color=254, charsize=0.6
   XYOUTS, 420, xtop-25, cmetime + ' UT',  /device, color=254, charsize=0.6
   sregnum = STRING (regnum, FORMAT='(I2)')
   sregnum = STRTRIM (sregnum, 2)
   print, 'regnum: ', regnum, ' sregnum: ', sregnum
   XYOUTS,   2,  5, 'Region: ' + sregnum, /device, color=254, charsize=0.6

   aoi_image = TVRD ()
   dotpos = STRPOS (cme_fits_file, '.fts')

;   gifname = cme_fits_file
;   STRPUT, gifname, '.gif', dotpos
;   print, 'Output Image : ', gifname

   print, "AOI file: ", aoi_file

   TVLCT, red, green, blue, /GET
   WRITE_GIF, aoi_file, aoi_image, red, green, blue

ENDWHILE  ;}
CLOSE, 11

   RETURN
END