;+
;NAME:		xmass
;
;PURPOSE:	Compute excess mass in a coronal image.
;
;CALLING SEQUENCE:	xmass, 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
;	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.
;
;HISTORY:	Created 3 June 1998, Andrew L. Stanger, HAO/NCAR
;-

FUNCTION xmass, imgref, imgcme, xdim, ydim, nbypix,	$
                xcen, ycen, pixrs, xloc, yloc, nloc, tscope


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

;cpmassb (namref, namcme, mx, cbth, pgmout, winpos, np)

;char	namref []	;--- REF image file name. 
;char	namcme []	;--- CME image file name. 
;int	mx		;--- Mass: =1: CME-REF otherwise, Absolute: CME. 
;float	cbth		;--- Coronal brightness threshold (B/Bsun).
;FILE *pgmout		;--- Output destination. 
;short winpos [][2]	;--- Window x/y coordinate positions of pixels in AOI. 
;int np			;--- # pixels in AOI. 

;{
;#define	MAXPIX	896
;#define RSUN 6.96e+10	;--- Solar radius (cm). 
;extern	char *malloc()
;extern	char *rdimg ()
;extern  int    cphskp(), hkline() ;--- Housekeeping decommutation & display. 
;extern	int    min (), max ()	;--- Integer min & max functions. 
;extern	float  xmin (), xmax ()	;--- Specialized float min & max functions. 
;extern	float  stray ()		;--- Declare stray light function. 
;extern	float  vignet ()	;--- Declare vignetting function. 
;extern	float  fcorona ()	;--- F-corona function. 
;extern	float  calcor ()	;--- Munro calibration correction function. 
;extern	double medsun ()	;--- Mass density (grams/cm2) per 1.0 Bsun. 
;
;static char	answer [2]	;--- Storage for y/n response. 
;char	*imgref, *imgcme
;
;
;int	cmefd, reffd	;--- File descriptors for CME & REF images. 
;int	ctype		;--- Data type for CME image. 
;int	i		;--- General loop index. 
;int	icp, irp	;--- Array indices for imgcme and imgref. 
;int	iix, iiy	;--- Index variables for pixel and line. 
;int	iixb		;--- X-axis index for REF image. 
;int	iiyb		;--- Y-axis index for REF image. 
;int	ipix
;int	ipc,irc,lpc,lrc	;--- Limits for CME image. 
;int	ipb,irb,lpb,lrb	;--- Limits for REF image. 
;int	irad		;--- Index for elec array. 
;int	ixs, iys	;--- Sun center shift of REF image relative to CME. 
;int	ixcm, iycm	;--- Center of mass coordinates. 
;int	ixcmg, iycmg	;--- Center of mass coordinates for display. 
;int	ixg, iyg	;--- Display coordinates of integrated pixel. 
;int	ix11,ix12,ix21,ix22	;--- Pixel limits for wedge. 
;int	iy11,iy12,iy21,iy22	;--- Line limits for wedge. 
;int	mv		;--- Vignetting flag: =1 vignetting has been removed. 
;int	n		;--- temporary debug variable. 
;int	nbyrec		;--- Number of bytes per record. 
;int	num, numun	;--- Number of values in kmass & unmass. 
;int	numneg		;--- Number of negative values in integration. 
;int	numth		;--- Number of values above 'cpth' brightness threshold.
;int	rtype		;--- Data type for REF image. 
;int	xocp = 63	;--- Window x origin. 
;int	yocp = 510	;--- Window y origin. 
;
;float	*xbbuf		;--- Data buffer: excess brightness. 
;float	*xmbuf		;--- Data buffer: excess mass. 
;float	bcme		;--- Brightness of CME image pixel. 
;float	bref		;--- Brightness of REF image pixel. 
;float	cbrite		;--- Coronal mass. 
;float	cbsum		;--- Integrated coronal brightness. 
;float	cbthresh	;--- Threshold level (B/Bsun). 
;float	cc		;--- Munro calibration correction factor. 
;float	cm2pix		;--- Area (cm**2) per pixel. 
;float	bcmesum		;--- Integrated coronal mass. 
;float	bcmesum0		;--- Integrated coronal mass above a threshold level. 
;float	cmemass		;--- CME image mass. 
;float	crrbrt		;--- Integrated brightness: (CME - REF) / REF. 
;float	cxcen, cycen	;--- Sun center for CME image. 
;float	fx, fy		;--- Floating variables for iix & iiy. 
;float	kbrite		;--- Integrated absolute or excess brightness value. 
;float	kmass		;--- Integrated absolute or excess mass value. 
;float	masspct		;--- % mass change. 
;float	negmass		;--- Negative mass. 
;float	negmassp	;--- Negative mass percentage. 
;float	negpct		;--- % Mass which is negative. 
;float	p11,p12,p21,p22	;--- Pixel limits for wedge. 
;float	r11,r12,r21,r22	;--- Line limits for wedge. 
;float	r, theta	;--- Radius & theta. 
;float	rcm, thcm	;--- (r,th) coordinates of center-of-mass.  
;float	brefsum		;--- Integrated coronal mass. 
;float	brefsum0	;--- Integrated coronal mass above a threshold level. 
;float	refmass		;--- REF image mass. 
;float	rxcen, rycen	;--- Sun center for REF image. 
;float	str, vig, fc	;--- Stray light, vignetting & F-corona  values. 
;float	xcm, ycm	;--- (x,y) coordinates of center-of-mass. 
;float	unmass		;--- Mass from unreliably low differences. 
;float	unmassp		;--- % Mass from unreliably low differences. 
;float	unpct		;--- Percentage from unreliably low differences. 
;float	wtheta		;--- Wedge theta coordinate. 
;
;double	gcm2		;--- grams per cm**2 along line of sight. 
;double	Nelec		;--- # electrons. 
;
;static float difunr = 5.0e-10	;--- Threshold for unreliable differences. 
;static double emass = 2.0e-24	;--- Mass per light-scattering electron. 
;
;--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   difunr   = 1.0e-10		; threshold for unreliable differences.
   emass    = 2.0e-24		; mass per light-scattering electron.
   RSUN     = 6.96e+10
   lstr     = 0
   lvig     = 0
   lcc      = 0
   ljc      = 0
   cbth     = 0.0
   cbthresh = cbth
   rspix    = 1.0 / pixrs
   pscale   = 1.0
   scroll   = 0.0
   ixs      = 0
   iys      = 0
   n        = 0

   ;--- 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.'

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

   cm2pix = (RSUN*rspix) * (RSUN*rspix)

   ;--- Initialize variables. 

   cc      = pscale
   fc      = 0.0
   str     = 0.0
   vig     = 1.0
   mx      = 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

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

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

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

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

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

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

   FOR ipix = 0, nloc - 1 DO	$
      BEGIN	;{
;>>>      fprintf (stderr, "xmass...pixel # %3d\n", ipix) 

      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, scroll, $
		     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 2) THEN 	$
         gcm2 = emdpb (r)	$			;--- grams/cm2/Bsun(pB).
      ELSE			$
      IF (tscope EQ 1) THEN 	$
	 BEGIN			
	 gcm2 = emdsun (r)				;--- grams/cm2/Bsun(B).
	 ljc = 1
	 lstr = 1
	 lvig = 1
	 END			$
      ELSE			$
	 gcm2 = emdsun (r)				;--- grams/cm2/Bsun.

      ;--- Use external stray/vignet functions 

;      IF (lcc) THEN  cc  = calcor (r) * pscale		;--- Munro cal. corr. 

      IF (ljc)  THEN cc  = 2.0 * pscale		;--- Jupiter calib. correction. 
      IF (lstr) THEN str = stray (r)		;--- stray light.
      IF (lvig) THEN vig = vignet (r)		;--- vignetting. 
      fc = fcorona (r, theta)			; F-corona.

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

      IF (mx  EQ  0) THEN		$
         BEGIN	  ;{

	 ;--- Pixel absolute brightness.

         cbrite  = ((float (imgcme [iix,iiy]) - str) * cc * vig) - fc
         bcme    = ((float (imgcme [iix,iiy]) - str) * cc * vig) - fc

         ;--- Integrated absolute brightness.

         bcmesum = bcmesum + bcme
         cbsum   = cbsum   + cbrite

;>>>         print, '   imgcme(',iix,',',iiy,'):', imgcme [iix,iiy], '  str:', $
;>>>         str, 'vig:', vig, '  cbrite: ', cbrite
 
         ;--- Integrated absolute brightess & mass.

         IF (cbrite > cbthresh) THEN	$
            BEGIN	;{
            numth    = numth + 1
            bcmesum0 = bcmesum0 + cbrite 
            cmemass  = cmemass + cbrite * gcm2 * cm2pix * cbrite
            END		;}

         xbbuf (i) = cbrite				;--- Save brightness.
         xmbuf (i) = float (cbrite * gcm2 * cm2pix)	;--- Save mass.

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

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

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

         ;--- Pixel brightness.

         cbrite  =  float (imgcme [iix,iiy] - imgref [iixb,iiyb]) * cc * vig
         bref    = ((float (imgref [iixb,iiyb]) - str) * cc * vig) - fc
         bcme    = ((float (imgcme [iix,iiy])   - str) * cc * vig) - fc

         ;--- Integrated brightness.

         cbsum   = cbsum   + cbrite
         brefsum = brefsum + bref
         bcmesum = bcmesum + bcme


         IF (cbrite GT cbthresh) THEN	$
 
            BEGIN	;{

            ;--- Brightness above threshold.

            numth    = numth + 1
            brefsum0 = brefsum0 + bref
            bcmesum0 = bcmesum0 + bcme

            ;--- Mass above threshold.

            refmass  = refmass + bref * gcm2 * cm2pix
            cmemass  = cmemass + bcme * gcm2 * cm2pix
            END		;}

         xbbuf (i) = cbrite				;--- Save brightness.
         xmbuf (i) = float (cbrite * gcm2 * cm2pix)	;--- Save mass.
         END	;}

;>>>   print, 'gcm2:', gcm2, ' cm2pix:', cm2pix, ' v:', v
;>>>   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, '# pixels : ',  num

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

   IF (lcc) THEN	$
      print, '>>> NOTE: Munro calibration correction has been applied. <<<'
   IF (ljc) THEN	$
      print, '>>> NOTE: Jupiter calibration correction has been applied. <<<'

   IF (brefsum GT 0.0) THEN			$
      crrbrt = (bcmesum - brefsum) / brefsum	$
   ELSE						$
      crrbrt = 0.0

   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, scroll, $
		  xcen, ycen, pixrs)
   IF (ierr NE 0) THEN  goto, oops
;   print, 'rcm: ', rcm, ' thcm: ', thcm, ' xcm: ', xcm, ' ycm: ', ycm, $
;	  ' scroll: ', scroll

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

   IF (mx EQ 0) THEN	$		;--- Absolute brightness/mass. 
      BEGIN	;{
;      print, '[CME brightness] (', numth, ' pixels) : ', bcmesum, ' B/Bsun'

      IF (cbthresh GT 0.0) THEN	$
         print, '[Brightness > ', cbthresh, 'B/Bsun] (', numth, ') :  ', $
	        bcmesum0, ' B/Bsun'

      print, '[Integrated absolute brightness] : ', cbsum, ' B/Bsun'

      print, '[Integrated absolute mass] :       ', kmass, ' grams <->', $
	     Nelec, 'electrons.'
      END	$ ;}
   ELSE				$
   IF (mx GT 0) THEN		$	;--- Excess brightness/mass. 
      BEGIN	;{
      print, '[CME brightness] (', numth, ' pixels) :', bcmesum, ' B/Bsun'
      print, '[REF brightness] (', numth, ' pixels) :', brefsum, ' B/Bsun'
      print, '[CME - REF brightness] :            ', cbsum, $
	     ' B/Bsun (excess brightness)'
      print, '[(CME - REF) / REF brightness] :', crrbrt, $
	     '        (delta B / B)'

      masspct = ((cmemass - refmass) / refmass) * 100.0
      print, 'CME image mass :          ', cmemass, ' grams'
      print, 'REF image mass :          ', refmass, ' grams'
      print, '[Integrated excess mass] :', kmass, ' grams <-> ', Nelec,	$
	     ' electrons (', 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, 'Unreliably low points [', unpct, '% (', numun, '/', num, ') with'
      print, '   -', difunr, ' < B < ', difunr, ']'
      print, 'contributed ', unmass, ' grams [', unmassp, '%%] to the total.'

      print, 'Negative mass points [', negpct, '% (', numneg, '/', num,	$
	     ') with B < 0.0]'
      print, 'contributed ', negmass, ' grams [', negmassp, '%] to the total.'

      END	;}

   print, 'The center of mass within this area is at ', rcm, ' Rsun   ', $
	  thcm, ' degrees.'
   print, '--------------------------------------------------------------------------------'

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