;+
; fm.pro   IDL procedure
;
; PURPOSE
;	Compute absolute/excess coronal mass.
;	Read FITS images and define a region of interest.
;
; SYNTAX
;	fm, 'cme_fits_file'
;	fm, 'cme_fits_file', 'ref_fits_file'
;
; PROCEDURES USED
;	fxpar		Get FITS parameter associated with a keyword
;	defroic		(tailored version of "defroi")
;	get_cal_name
;	readfits	Read FITS file
;	write_gif	Write GIF file
;	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).
;-

PRO fm06, cme_fits_file, ref_fits_file

   ;--- Initialize variables.

   i      = 0L
   nloc   = 0L
   regnum = 0

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

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

   ;--- 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
   xcen_cme     = fxpar (hdu_cme, "CRPIX1")
   ycen_cme     = fxpar (hdu_cme, "CRPIX2")
   cdelt1_cme   = fxpar (hdu_cme, "CDELT1")
   crradius_cme = fxpar (hdu_cme, "CRRADIUS")
   pixrs 	= crradius_cme
   rsun         = fxpar (hdu_cme, "RSUN")
   crota1_cme   = fxpar (hdu_cme, "CROTA1")
   scroll_cme   = fxpar (hdu_cme, "SCROLL")
   roll_cme = -crota1_cme + scroll_cme
   print, 'fm> crradius_cme: ', crradius_cme, '  pixrs: ', pixrs
   print, 'fm> rsun        : ', rsun

   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)
   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 = 7
      dspcme = imgcme ;
      END

   IF (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/'
;	    dte = STR2UTC (cmedate)
;	    yymmdd = UTC2YYMMDD (dte)
	    vig_fn = get_cal_name (sd + 'C2_cl*vig*.dat', yymmdd)
	    IF (vig_fn NE '') THEN	$
	    BEGIN		;{
	       vig = FLTARR (1024,1024)
	       OPENR, lucal, vig_fn, /GET_LUN
	       READU, lucal, vig
;	       IF (szlog (0) NE 0)   THEN PRINTF, lulog, 'Used ' + vig_fn
	       PRINT, 'Used ' + vig_fn
	       CLOSE, lucal
	       FREE_LUN, lucal
	       dspcme = imgcme / vig 
	    ENDIF	$	;}
	    ELSE BEGIN		;{
	       print,'ERROR: c2_calibrate - No ' + sd + 'C2_cl*vig*.dat file'
;	       IF (szlog(0) NE 0)   THEN	$
;	          PRINTF, lulog, 'ERROR: c2_calibrate - No ' + sd +	$
;				 'C2_cl*vig*.dat file'
	          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/'
;	    dte = STR2UTC (cmedate)
;	    yymmdd = UTC2YYMMDD (dte)
;	    print, 'fm> dte: ', dte, '  yymmdd: ', yymmdd
	    vig_fn = get_cal_name (sd + 'C3_cl*vig*.dat', yymmdd)
	    IF (vig_fn NE '') THEN		$
	    BEGIN	;{
	       vig = FLTARR (1024,1024)
	       OPENR, lucal, vig_fn,/GET_LUN
	       READU, lucal, vig
;	       IF (szlog (0) NE 0)   THEN PRINTF, lulog, 'Used ' + vig_fn
	       PRINT, 'Used ' + vig_fn
	       CLOSE, lucal
	       FREE_LUN, lucal
	       dspcme = imgcme / vig 
	    ENDIF	$;}
	    ELSE BEGIN	;{
	       print, 'ERROR: c3_calibrate - No ' + sd + 'C3_cl*vig*.dat file'
;	       IF (szlog (0) NE 0)   THEN	$
;	          PRINTF, lulog, 'ERROR: c3_calibrate - No ' + sd + 	$
;				 'C3_cl*vig*.dat file'
	          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		; 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
      xcen_ref     = fxpar (hdu_ref, "CRPIX1")
      ycen_ref     = fxpar (hdu_ref, "CRPIX2")
      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")

      refday    = STRMID (refdate, 0, 2)
      refmonth  = STRMID (refdate, 3, 2)
      refyear   = STRMID (refdate, 6, STRLEN (refdate) - 6)
      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
      IF (err NE 0) THEN RETURN

      if (STRPOS (telescop_ref, "SMM C/P") NE -1)  THEN	$	; SMM C/P image.
         BEGIN
	 tscope = 7
         dspref = imgref
         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, '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/'
;	       dte = STR2UTC (cmedate)
;	       yymmdd = UTC2YYMMDD (dte)
	       vig_fn = get_cal_name (sd + 'C2_cl*vig*.dat', yymmdd)
	       IF (vig_fn NE '') THEN	$
	       BEGIN		;{
	          vig = FLTARR (1024,1024)
	          OPENR, lucal, vig_fn, /GET_LUN
	          READU, lucal, vig
;	          IF (szlog (0) NE 0)   THEN PRINTF, lulog, 'Used ' + vig_fn
	          PRINT, 'Used ' + vig_fn
	          CLOSE, lucal
	          FREE_LUN, lucal
	          dspref = imgref / vig 
	       ENDIF	$	;}
	       ELSE BEGIN		;{
	          print,'ERROR: c2_calibrate - No ' + sd + 'C2_cl*vig*.dat file'
;	          IF (szlog(0) NE 0)   THEN	$
;	             PRINTF, lulog, 'ERROR: c2_calibrate - No ' + sd +	$
;				    'C2_cl*vig*.dat file'
	             PRINT, '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/'
;	       dte = STR2UTC (cmedate)
;	       yymmdd = UTC2YYMMDD (dte)
;	       print, 'fm> dte: ', dte, '  yymmdd: ', yymmdd
	       vig_fn = get_cal_name (sd + 'C3_cl*vig*.dat', yymmdd)
	       IF (vig_fn NE '') THEN		$
	       BEGIN	;{
	          vig = FLTARR (1024,1024)
	          OPENR, lucal, vig_fn,/GET_LUN
	          READU, lucal, vig
;	          IF (szlog (0) NE 0)   THEN PRINTF, lulog, 'Used ' + vig_fn
	          PRINT, 'Used ' + vig_fn
	          CLOSE, lucal
	          FREE_LUN, lucal
	          dspref = imgref / vig 
	       ENDIF	$;}
	       ELSE BEGIN	;{
	          print, 'ERROR: c3_calibrate - No ' + sd + 'C3_cl*vig*.dat file'
;	          IF (szlog (0) NE 0)   THEN	$
;	             PRINTF, lulog, 'ERROR: c3_calibrate - No ' + sd + 	$
;				    'C3_cl*vig*.dat file'
	             PRINT, '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, 'dif min/max: ', imin, imax
      imin = imin * 1.5

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

      IF (ABS (imin) LT abs (imax))  THEN imax = - imin
      IF (ABS (imin) GT abs (imax))  THEN imin = - imax
      print, 'dif min/max: ', imin, imax

      dspdif = BYTSCL (dspdif, min=imin, max=imax, top=250)
      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.
      imin = MIN (dspcme, max=imax)
;      print, '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, '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)

      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 = 'smm_cp'
      ELSE: scope_name = 'Unknown'
   ENDCASE

   ;--- Action Loop

loop: print, 'Request Action'

   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

   ;--- Request Action Command.

   print, 'Refresh Display                       [d]'
   print, 'Keyboard input to define curved wedge [kc]'
   print, 'Mouse cursor   to define curved wedge [mc]'
   print, 'Keyboard input to define wedge        [kw]'
   print, 'Mouse cursor   to define wedge        [mw]'
   print, 'Mouse cursor   to define polygon      [p]'
   print, 'Read disk file to define polygon      [f]'
   print, 'Quit program                          [q]'

   action = ''
   read, 'action: ', action

   CASE action OF
   'd': BEGIN
;	IF (mx EQ 0) THEN TV, dspcmeb	$	; Display CME image.
;	ELSE TV, dspdif				; Display difference image.
	regnum = regnum - 1
	GOTO, loop
	END

   'kc': BEGIN
	kb_cyl5,      xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $
		      regnum, log_file, pv_file, xv, yv, nver, polypix
	END

   'mc': BEGIN
	mouse_cyl5,   xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $
		      regnum, log_file, pv_file, xv, yv, nver, polypix
	END

   'kw': BEGIN
	kb_wedge5,    xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $
		      regnum, log_file, pv_file, xv, yv, nver, polypix
	END

   'mw': BEGIN
	mouse_wedge5, xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $
		      regnum, log_file, pv_file, xv, yv, nver, polypix
        END

   'p': BEGIN
	mouse_poly5,  xdim_cme, ydim_cme,	$
		      regnum, log_file, pv_file, xv, yv, nver, polypix
        END

   'f': BEGIN
	file_poly5,   xdim_cme, ydim_cme,	$
		      regnum, log_file, pv_file, xv, yv, nver, polypix
	END
   'q': RETURN
   ELSE: RETURN
   ENDCASE

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

   size_roi = SIZE (polypix)

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

   print, 'size (polypix): ', size_roi
   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 = xmass06 (log_file, scope_name, ref_fits_file, cme_fits_file,	$
	       imgref, imgcme,						$
	       xdim_cme, ydim_cme, nbypix, xcen_cme, ycen_cme,		$
	       pixrs, rsun, xloc, yloc, nloc, tscope)

   print, 'xm: ', xm

   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

   GOTO, loop

   RETURN
END