;+
;NAME:		xmass
;
;PURPOSE:	Compute excess mass in a coronal image.
;
;CALLING SEQUENCE:	xmass,	log_file, scope_name,
;				ref_name, cme_name, imgref, imgcme,
;				xdim, ydim, nbypix,
;				xcen, ycen, pixrs, xloc, yloc, nloc, tscope
;
;INPUT:	imgref	- 2D reference (pre-event) image array.
;	imgcme	- 2D event image array.
;	xdim	- x-axis dimension
;	ydim	- y-axis dimension
;	nbypix	- # bytes per pixel
;	xcen	- x-axis center
;	ycen	- y-axis center
;	pixrs	- # pixels per solar radius
;	rsun    - arcseconds per solar radius
;	roll	- effective roll (Solar North w.r.t. +X axis)
;	xloc	- x-axis location for pixels within ROI (region of interest).
;	yloc	- y-axis location for pixels within ROI (region of interest).
;	nloc	- number of pixels within ROI
;	tscope	- =1, SMM C/P (B), =2, MK3 (pB).
;
;RETURNS:	excess mass (grams) for the ROI pixels provided.
;
;PROCEDURES USED:
;	rcoord
;	emdpb
;	emdsun
;	fcorona
;
;HISTORY:
;	Created 3 June 1998, Andrew L. Stanger, HAO/NCAR
;	 9 Sep 1998: Remove F-corona ONLY from SMM C/P images.
;	19 Oct 1998: Add scope_name parameter.
;	 8 Feb 1999: replace numth with num for brightness integration.
;	 2 Apr 1999: add rsun parameter.
;	             Replace stray & vignet with cpstray & cpvignet.
;	21 Apr 1999: calibration correction (cc) no longer applied
;	             to stray light TWICE !
;	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.
;-

FUNCTION xmass08, log_file, scope_name, ref_fits_file, cme_fits_file,	$
		imgref, imgcme,	bunit_cme, bunit_ref,			$
		xdim, ydim, nbypix, xcen, ycen, pixrs, rsun, roll,	$
		xloc, yloc, nloc, tscope

; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

   ;--- Determine whether we're doing absolute or excess mass.

   IF (N_ELEMENTS (imgref) EQ 1)  THEN  mx = 0  ELSE  mx = 1

   difunr   = 1.0e-10		; threshold for unreliable differences.
   emass    = 2.0e-24		; mass per light-scattering electron.
   RSUNCM   = 6.96e+10		; Average solar radius (cm).
   lstr     = 0
   lvig     = 0
   lcc      = 0
   ljc      = 0
   cbth     = 0.0
   cbthresh = cbth
   rspix    = 1.0 / pixrs
   pscale   = 1.0
   ixs      = 0
   iys      = 0
   n        = 0

;   print, 'pixrs: ', pixrs

   OPENU,  LUNLOG, log_file, /GET_LUN, /APPEND

   PRINTF, LUNLOG, scope_name + ' : Mass Measurement   [xmass v0.8]'
   PRINTF, LUNLOG, 'CME image : ', cme_fits_file
   IF (mx) THEN	$
      PRINTF, LUNLOG, 'REF image : ', ref_fits_file

   ;--- Stray Light & Vignetting removal required ? 

;   if (lstr EQ 2 AND lvig EQ 2) THEN mv = 1  ELSE mv = 0

   print, 'Please wait ... mass computation in progress.'
   print, '--------------------------------------------------------------------------------'
   PRINTF, LUNLOG, '--------------------------------------------------------------------------------'

   ;--- Compute the cross-sectional area (square cm) for one pixel. 

   cm2pix = (RSUNCM*rspix) * (RSUNCM*rspix)

   ;--- Initialize variables. 

   cc      = pscale
   fc      = 0.0
   str     = 0.0
   vig     = 1.0
   xocp    = 0
   yocp    = 0
   Nelec   = 0.0
   cmemass = 0.0
   kbrite  = 0.0
   kmass   = 0.0
   negmass = 0.0
   refmass = 0.0
   unmass  = 0.0
   xcm     = 0.0
   ycm     = 0.0

   num     = 0
   numneg  = 0
   numun   = 0
   numth   = 0

;   print, 'ixs: ', ixs, '  iys: ', iys
;   print, 'nloc: ', nloc, '  x/y loc [beg]:', xloc [0], yloc [0],	$
;	  '  x/y loc [end]:', xloc [nloc-1], yloc [nloc-1]

   ;--- Clear data buffer. 

   xbbuf = fltarr (nloc)
   xmbuf = fltarr (nloc)
   xbbuf [*] = 0.0
   xmbuf [*] = 0.0

   ;--- Read REF and CME images into memory arrays. 

;   IF (mx > 0) THEN imgref =  rdimg (namref)
;   imgcme = rdimg (namcme)

   ;--- ******* Pixel Loop ******* 

   i        = -1L
   num      = 0L
   numth    = 0L
   cbsum    = 0.0
   brefsum  = 0.0
   brefsum0 = 0.0
   bcmesum  = 0.0
   bcmesum0 = 0.0

   refbrtsum  = 0.0
   refmassum  = 0.0
   refbrtsum0 = 0.0
   refmassum0 = 0.0

   cmebrtsum  = 0.0
   cmemassum  = 0.0
   cmebrtsum0 = 0.0
   cmemassum0 = 0.0

   excbrtsum  = 0.0
   excmassum  = 0.0

   crrbrt   = 0.0

   FOR ipix = 0L, nloc - 1 DO	$
      BEGIN	;{

      ixg = xloc [ipix]		;--- Window x-coordinate. 
      iyg = yloc [ipix]		;--- Window y-coordinate. 

      iix = ixg - xocp			;--- Image x-coordinate. 
      iiy = yocp - iyg			;--- Image y-coordinate. 
      iiy = iyg - yocp			;--- Image y-coordinate.

      fx = float (iix)			;--- Float image x-coordinate. 
      fy = float (iiy)			;--- Float image y-coordinate. 

      ierr = rcoord (r, theta, fx, fy, -1, roll, xcen, ycen, pixrs)
      IF (ierr NE 0) THEN goto, oops

;      print, 'fx:', fx, ' fy:', fy, ' r:', r, ' theta:', theta

      ;--- Draw dot to indicate pixel has been integrated. 

;      gl_wvec (ixg, iyg, ixg, iyg, ov_blue)

;>>> print, 'fx:', fx, ' fy:', fy, ' . . . passed r,th test.'

      i = i + 1
      num = num + 1
;      xbbuf [i] = 0.0
;      xmbuf [i] = 0.0

      IF (tscope EQ 6) THEN 	$
         gcm2 = emdpb (r)	$		;--- grams/cm2/Bsun(pB).
      ELSE			$
      IF (tscope EQ 7) THEN 	$
	 BEGIN			
	 IF (STRPOS (bunit_cme, 'MASS') EQ -1) THEN	$
	    BEGIN
	    gcm2 = emdsun (r)			;--- grams/cm2/Bsun(B).
	    ljc  = 1
	    lstr = 1
	    lvig = 1
            IF (ljc)  THEN cc  = 2.0 * pscale	;--- Jupiter calib. correction. 
            IF (lstr) THEN str = cpstray  (r * rsun)	;--- Stray light.
            IF (lvig) THEN vig = cpvignet (r * rsun)	;--- Vignetting. 
            fc = fcorona (r, theta)		;--- F-corona.
	    END			$
         ELSE			$		;--- Mass image.
	    BEGIN
	    cc     = 1.0
	    gcm2   = 1.0
	    cm2pix = 1.0
	    END
	 END			$
      ELSE			$
	 gcm2 = emdsun (r)			;--- grams/cm2/Bsun.

      ;--- Absolute mass: Remove stray light from CME image. 

      IF (mx  EQ  0) THEN		$
         BEGIN	  ;{

	 ;--- Pixel absolute brightness & mass.

         cmebrt = ((float (imgcme [iix,iiy]) * cc - str) * vig) - fc
         xbbuf [i] = cmebrt			;--- Save pixel brightness.

	 cmemas = cmebrt * gcm2 * cm2pix
	 xmbuf [i] = cmemas			;--- Save pixel mass.

         ;--- Integrated absolute brightness & mass.

         cmebrtsum = cmebrtsum + cmebrt
	 cmemassum = cmemassum + cmemas

;>>>         print, '   imgcme (',iix,',',iiy,'):', imgcme [iix,iiy], '  str:',$
;>>>         str, 'vig:', vig, '  cmebrt: ', cmebrt


         ;--- Integrated brightness & mass above threshold.

         IF (cmebrt > cbthresh) THEN	$
            BEGIN	;{
            numth    = numth + 1
            cmebrtsum0 = cmebrtsum0 + cmebrt
	    cmemassum0 = cmemassum0 + cmemas
            END		;}

;>>>         print, 'gcm2: ', gcm2, '  cm2pix: ', cm2pix, $
;>>>                '  bcme: ', bcme, '  data [', i, ']]: ', data [i]
 
         END	$ ;}

      ;--- Excess mass: compute pixel coordinate for REF image. 

      ELSE	$
         BEGIN	;{
         iiyb = iiy - iys
         iixb = iix - ixs

         ;--- Pixel brightness.

         cmebrt  = ((float (imgcme [iix,iiy])   * cc - str) * vig) - fc
         refbrt  = ((float (imgref [iixb,iiyb]) * cc - str) * vig) - fc
         excbrt  =   float (imgcme [iix,iiy] - imgref [iixb,iiyb]) * cc * vig

	 cmemas = cmebrt * gcm2 * cm2pix
	 refmas = refbrt * gcm2 * cm2pix
	 excmas = excbrt * gcm2 * cm2pix

         ;--- Integrated brightness.

         cmebrtsum = cmebrtsum + cmebrt
         refbrtsum = refbrtsum + refbrt
         excbrtsum = excbrtsum + excbrt

         cmemassum = cmemassum + cmemas
         refmassum = refmassum + refmas
         excmassum = excmassum + excmas

         IF (excbrt GT cbthresh) THEN	$
 
            BEGIN	;{

            ;--- Brightness above threshold.

            numth    = numth + 1
            refbrtsum0 = refbrtsum0 + refbrt
            cmebrtsum0 = cmebrtsum0 + cmebrt

            ;--- Mass above threshold.

            refmassum0 = refmassum0 + refmas
            cmemassum0 = cmemassum0 + cmemas
            END		;}

         xbbuf [i] = excbrt			;--- Save pixel brightness.
         xmbuf [i] = excmas			;--- Save pixel mass.
         END	;}

;>>>   print, 'gcm2:', gcm2, ' cm2pix:', cm2pix
;>>>   print, 'imgcme:', imgcme [iix,iiy], ' imgref:', imgref [iixb,iiyb]
;>>>   print, 'xmbuf [', iix, iiy, ']: ', xmbuf [iix,iiy]

      ;--- K-corona brightness & mass.

      kbrite = kbrite + xbbuf [i]
      kmass  = kmass  + xmbuf [i]

      IF (mx GT 0) THEN		$
         BEGIN	;{
         IF (xbbuf [i] LT difunr AND xbbuf [i] GT (-difunr)) THEN	$
            BEGIN	;{			;--- Mass below threshold.
            unmass = unmass + xmbuf (i)
            numun = numun + 1
            END		;}

         IF (xmbuf [i] LT 0.0) THEN	$
            BEGIN	;{			;--- Negative mass.
            negmass = negmass + xmbuf [i]
            numneg = numneg + 1
            END		;}
         END	;}

      ;--- Center of mass: integrated values.

      xcm = xcm + xmbuf [i] * fx
      ycm = ycm + xmbuf [i] * fy

pixend:    n = n + 1		;--- print, 'pixend...iix;', iix 

recend: n = n + 1		;--- print, 'recend...iiy:', iiy 

      END	;}--- End pixel loop. 

   print, ' '
   PRINTF, LUNLOG, ''

;>>>   print, 'F-corona has NOT been removed.' 

   IF (lcc) THEN	$
      BEGIN
      print, '>>> NOTE: Munro calibration correction has been applied. <<<'
      PRINTF, LUNLOG, '>>> NOTE: Munro calibration correction has been applied. <<<'
      END

   IF (ljc) THEN	$
      BEGIN
      print, '>>> NOTE: Jupiter calibration correction has been applied. <<<'
      PRINTF, LUNLOG, '>>> NOTE: Jupiter calibration correction has been applied. <<<'
      END

   IF (refbrtsum NE 0.0) THEN			$
      crrbrt = (cmebrtsum - refbrtsum) / refbrtsum	$
   ELSE						$
      crrbrt = 0.0
   print, '>>>>>>>>>>>>>>>> crrbrt: ', crrbrt

   Nelec = double (kmass) / emass

   ;--- Center of mass: normalize using integrated mass. 

   IF (kmass  NE  0.0) THEN	$
      BEGIN	;{
      ixcm = FIX (xcm / kmass)
      iycm = FIX (ycm / kmass)
      END	;}

   ;--- Cartesian coordinates of center of mass.

   xcm = float (ixcm)
   ycm = float (iycm)

   ;--- Convert to polar coordinates.

   ierr = rcoord (rcm, thcm, xcm, ycm, -1, roll, xcen, ycen, pixrs)
   IF (ierr NE 0) THEN  goto, oops
;   print, 'rcm: ', rcm, ' thcm: ', thcm, ' xcm: ', xcm, ' ycm: ', ycm, $
;	  ' roll: ', roll

   ;--- Draw "+" on display to indicate center of mass location. 
   ;--- Draw lines connecting 4 wedge corners. 

;   IF (igr EQ 1) THEN	$
;      BEGIN	;{
;      ixcmg = ixcm/nres + xocp
;      iycmg = yocp - iycm/nres
;      gl_wvec (ixcmg-5, iycmg, ixcmg+5, iycmg, ov_red)
;      gl_wvec (ixcmg, iycmg-5, ixcmg, iycmg+5, ov_red)
;
;      ix11 = FIX (p11) / nres + xocp
;      ix12 = FIX (p12) / nres + xocp
;      ix21 = FIX (p21) / nres + xocp
;      ix22 = FIX (p22) / nres + xocp
;      iy11 = yocp - FIX (r11) / nres
;      iy12 = yocp - FIX (r12) / nres
;      iy21 = yocp - FIX (r21) / nres
;      iy22 = yocp - FIX (r22) / nres
;      END	;}

   print, ' '
   PRINTF, LUNLOG, ''

   IF (mx EQ 0) THEN	$		;--- Absolute brightness/mass. 
      BEGIN	;{

      IF (cbthresh GT 0.0) THEN	$
	 BEGIN
         print, '[Brightness > ', cbthresh, 'B/Bsun] (', num, ') :  ', $
	        bcmesum0, ' B/Bsun'
         PRINTF, LUNLOG, '[Brightness > ', cbthresh, 'B/Bsun] (', num, ') :  ', $
	        cmebrtsum0, ' B/Bsun'
         END

      print, FORMAT = '("[Integrated absolute brightness] : ", e10.3, " B/Bsun")', cmebrtsum
      PRINTF, LUNLOG, FORMAT = '("[Integrated absolute brightness] : ", e10.3, " B/Bsun")', cmebrtsum

      print, FORMAT = '("[Integrated absolute mass] :       ", e10.3, " grams <-> ", e10.3, " electrons")', kmass, Nelec
      PRINTF, LUNLOG, FORMAT = '("[Integrated absolute mass] :       ", e10.3, " grams <-> ", e10.3, " electrons")', kmass, Nelec

      print, ' '
      PRINTF, LUNLOG, ''
      END	$ ;}
   ELSE				$
   IF (mx GT 0) THEN		$	;--- Excess brightness/mass. 
      BEGIN	;{
      print, FORMAT = '("------- Brightness Integration -------    # pixels : ", I7)', num
      PRINTF, LUNLOG, FORMAT = '("------- Brightness Integration -------    # pixels : ", I7)', num
      print, FORMAT = '(" CME              : ", e10.3, " B/Bsun")', cmebrtsum
      PRINTF, LUNLOG, FORMAT = '(" CME              : ", e10.3, " B/Bsun")', cmebrtsum

      print, FORMAT = '(" REF              : ", e10.3, " B/Bsun")', refbrtsum
      PRINTF, LUNLOG, FORMAT = '(" REF              : ", e10.3, " B/Bsun")', refbrtsum

      print, FORMAT = '(" CME - REF        : ", e10.3, " B/Bsun")', excbrtsum 
      PRINTF, LUNLOG, FORMAT = '(" CME - REF        : ", e10.3, " B/Bsun")', excbrtsum 

      print, FORMAT = '("(CME - REF) / REF : ", f7.4,  " %", 13x, "(delta B / B)")', crrbrt * 100.0
      PRINTF, LUNLOG, FORMAT = '("(CME - REF) / REF : ", f8.4,  " %", 13x, "(delta B / B)")', crrbrt * 100.0

      masspct = ((cmemassum - refmassum) / refmassum) * 100.0
      print, ' '
      PRINTF, LUNLOG, ''

      print, FORMAT = '("---------- Mass Integration ----------    # pixels : ", I7)', num
      PRINTF, LUNLOG, FORMAT = '("---------- Mass Integration ----------    # pixels : ", I7)', num
      print, FORMAT = '(" CME              : ", e10.3, " grams")', cmemassum
      PRINTF, LUNLOG, FORMAT = '(" CME              : ", e10.3, " grams")', cmemassum

      print, FORMAT = '(" REF              : ", e10.3, " grams")', refmassum
      PRINTF, LUNLOG, FORMAT = '(" REF              : ", e10.3, " grams")', refmassum

      print, FORMAT = '(" CME - REF        : ", e10.3, " grams <-> ", e10.3, " electrons ")', kmass, Nelec
      PRINTF, LUNLOG, FORMAT = '(" CME - REF        : ", e10.3, " grams <-> ", e10.3, " electrons ")', kmass, Nelec

      print, FORMAT = '("(CME - REF) / REF : ", f8.4,  " %", 13x, "(delta M / M)")', masspct
      PRINTF, LUNLOG, FORMAT = '("(CME - REF) / REF : ", f8.4,  " %", 13x, "(delta M / M)")', masspct

      unpct    = ( float (numun)  / float (num) ) * 100.0
      negpct   = ( float (numneg) / float (num) ) * 100.0
      unmassp  = (unmass  / kmass) * 100.0
      negmassp = (negmass / kmass) * 100.0
      print, ' '
      PRINTF, LUNLOG, ''

      print, FORMAT = '("Unreliably low points [", F7.2, " %", " (", I6, "/", I6, ")]")', unpct, numun, num
      print, FORMAT = '("with -", e10.3, " < B < ", e10.3, "]")', difunr, difunr
      print, FORMAT = '("contributed ", e10.3, " grams [", f8.3, "%] to the total")', unmass, unmassp
      PRINTF, LUNLOG, FORMAT = '("Unreliably low points [", F7.2, " %", " (", I6, "/", I6, ")]")', unpct, numun, num
      PRINTF, LUNLOG, FORMAT = '("with -", e10.3, " < B < ", e10.3, "]")', difunr, difunr
      PRINTF, LUNLOG, FORMAT = '("contributed ", e10.3, " grams [", f8.3, "%] to the total")', unmass, unmassp

      print, ' '
      PRINTF, LUNLOG, ''

      print, FORMAT = '("Negative mass  points [", f7.2, " %", " (", I6, "/", I6, ")]")', negpct, numneg, num
      print, FORMAT = '("contributed ", e10.3, " grams [", f8.3, "%] to the total")', negmass, negmassp
      PRINTF, LUNLOG, FORMAT = '("Negative mass  points [", f7.2, " %", " (", I6, "/", I6, ")]")', negpct, numneg, num
      PRINTF, LUNLOG, FORMAT = '("contributed ", e10.3, " grams [", f8.3, "%] to the total")', negmass, negmassp


      print, ' '
      PRINTF, LUNLOG, ''
      END	;}

   print, FORMAT = '("Center of mass : ", F7.2, " Rsun  ", F6.1, " degrees")', $
                   rcm, thcm
   PRINTF, LUNLOG, FORMAT = '("Center of mass : ", F7.2, " Rsun  ", F6.1, " degrees")', $
                   rcm, thcm
   print, '--------------------------------------------------------------------------------'
   PRINTF, LUNLOG, '--------------------------------------------------------------------------------'

   fx = 0.0
   fy = 447.0
   ierr = rcoord (r, theta, fx, fy, -1, roll, xcen, ycen, pixrs)
;   print, '[0,447] -> [', r, ',', theta, ']'

   CLOSE, LUNLOG
   FREE_LUN, LUNLOG

;   IF (mx GT 0) THEN CLOSE, reffd
;   CLOSE, cmefd
   RETURN, kmass 

   ;--- Coordinate out-of-range error. 

oops:	$
   print, 'Invalid range: radius or theta'
   RETURN, -10
END ;}

;--- xmin (x1, x2, x3, x4)
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;* Find minimum of 4 floating point variables.
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;* Andrew L. Stanger
;* 1 March 1984
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

FUNCTION xmin, x1, x2, x3, x4

;float x1, x2, x3, x4

;{
   x = 0.0

   x = x1
   IF (x GT x2) THEN  x = x2
   IF (x GT x3) THEN  x = x3
   IF (x GT x4) THEN  x = x4

   RETURN, x
END	;}

;--- xmax (x1, x2, x3, x4)
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;* Find maximum of 4 floating point variables.
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;* Andrew L. Stanger
;* 1 March 1984
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
FUNCTION xmax, x1, x2, x3, x4

;float x1, x2, x3, x4

;{
   x = 0.0

   x = x1
   IF (x GT x2) THEN x = x2
   IF (x GT x3) THEN x = x3
   IF (x GT x4) THEN x = x4

   RETURN, x

END	;}