;+
; fr.pro   IDL procedure
;
; PURPOSE
;	Read FITS images and define a region of interest.
;
; SYNTAX
;	froi, 'fits_file'
;
; PROCEDURES USED
;	readfits
;	writefits
;	fxaddpar
;
; HISTORY
;	Author: Andrew L. Stanger   2 June 1998
;-

PRO fr, cme_fits_file, ref_fits_file

   ;--- Initialize variables.

   mx       = 0			; =0 cme image only, =1 cme & ref images

   xdim     = 0			; x-axis dimension
   ydim     = 0			; y-axis dimension
   nbypix   = 0			; # bytes/pixel
   xcen     = 0			; x-axis center
   ycen     = 0			; y-axis center
   crradius = 0.0		; # pixels/Rsun

   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

   ;--- Read FITS image.

   imgcme   = readfits (cme_fits_file, hdu_cme)

   xdim     = fxpar (hdu_cme, "NAXIS1")
   ydim     = fxpar (hdu_cme, "NAXIS2")
   nbypix   = fxpar (hdu_cme, "BITPIX") * 8
   xcen     = fxpar (hdu_cme, "CRPIX1")
   ycen     = fxpar (hdu_cme, "CRPIX2")
   crradius = fxpar (hdu_cme, "CRRADIUS")
   date     = fxpar (hdu_cme, "DATE-OBS")
   time     = fxpar (hdu_cme, "TIME-OBS")
   date_time = date + ' ' + time
   solar_info = pb0r (date_time,/arcsec)	; Get solar P, B0, Rsun.
   Pangle  = solar_info [0]			; P-angle
   B0      = solar_info [1]			; B0 angle
   Rsunarc = solar_info [2]			; Solar radius (arcsec).

   IF (N_PARAMS() EQ 2) THEN	$
      BEGIN		;{
      imgref = readfits (ref_fits_file, hdu_ref)
      xdim_ref     = fxpar (hdu_cme, "NAXIS1")
      ydim_ref     = fxpar (hdu_cme, "NAXIS2")
      nbypix_ref   = fxpar (hdu_cme, "BITPIX") * 8
      xcen_ref     = fxpar (hdu_cme, "CRPIX1")
      ycen_ref     = fxpar (hdu_cme, "CRPIX2")
      crradius_ref = fxpar (hdu_cme, "CRRADIUS")

      IF (xdim NE xdim_ref) THEN RETURN, -1
      IF (ydim NE ydim_ref) THEN RETURN, -2
      IF (nbypix NE nbypix_ref) THEN RETURN, -3
      IF (xcen NE xcen_ref) THEN RETURN, -4
      IF (ycen NE ycen_ref) THEN RETURN, -5
      mx = 1
      lct, '/home/stanger/color/dif.lut'
      tvscl, imgcme - imgref
      END	$	;}
   ELSE		$
      BEGIN	;{
      mx = 0
      imgref = 0
      lct, '/home/stanger/color/bwy.lut'
      tvscl, imgcme
      END	;}

   pixrs = crradius

   print, 'fr> crradius: ', crradius, '  pixrs: ', pixrs

   ;--- Display image.

   tvscl, imgcme

   ;--- Define Region-of-interest (ROI).

   polypix = DEFROI (xdim, ydim, xv, yv)	; array of pixel coordinates.

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

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

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

   size_xv = SIZE (xv)
   nv = size_xv [1]
   print, 'size [xv]: ', size [xv]
   print, 'size [yv]: ', size [yv]

   FOR i = 0, nv - 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 = xmass (imgref, imgcme, xdim, ydim, nbypix,	$
          xcen, ycen, pixrs, xloc, yloc, nloc)

   print, 'xm: ', xm

END