function scc_cme_massimg2total,fn,roi,MAXSCL=maxscl,MINSCL=minscl,SAVE=save, $
         SECTOR=SECTOR,RADII=RADII,ANGLES=ANGLES, batch=batch, HT_INFO=HT_INFO, $
	 VERSION=version, PAN=pan
;+
; $Id: scc_cme_massimg2total.pro,v 1.12 2012/04/12 20:05:52 nathan Exp $
;
; NAME:
;	scc_CME_MASSIMG2TOTAL
;
; PURPOSE:
;	This function calculates the mass of a CME from mass images.
;
; CATEGORY:
;	CME
;
; CALLING SEQUENCE:
;	Result = scc_CME_MASSIMG2TOTAL(Fn)
;
; INPUTS:
;	Fn:	String containing the filename of the CME mass image
;
; KEYWORD PARAMETERS:
;	MINSCL:	Set this keyword with the value to use for the minimum value
;		in scaling the image. The default value is to use -1.e-11 msb
;
;	MAXSCL:	Set this keyword with the value to use for the maximum value
;		in scaling the image. The default value is to use +1.e-11 msb
;
;	SAVE:	Set this keyword with the filename to save the mass information to
;	SECTOR:	Set this keyword if the ROI is a sector, centered on the sun
;		The default is to draw the boundary of the CME using the cursor
;	RADII:	Set this keyword with a 2-element array of the inner and
;		outer radii (in solar radii) for the sector ROI
;	ANGLES:	Set this keyword with a 2-element array of the left and right
;		hand boundaries (viewed from sun-center) for the sector ROI
;   	    	(PA = deg CCW from solar north)
;   	HT_INFO=[PA, width, height] Used in output 
;   	/BATCH	Do not prompt for user input.
;   	VERSION=    Return version from Id
;   	PAN=	Reduce/expand image size by this fraction before displaying
;
; OUTPUTS:
;	This function returns the mass calculated for the ROI selected
;   	roi 	ROI used is returned in this argument
;
; RESTRICTIONS:
;	Only tested on C3, but should work for any telescope
;
; EXTERNAL CALLS:
;       mREADFITS, scc_SUN_CENTER, scc_ROI_SECTOR, AWIN, one2two
;
; PROCEDURE:
;	The files for the CME mass image is read in.  The image is displayed and XLOADCT
;	is called to permit the contrast to be adjusted.  Then RDPIX is called to permit
;	individual pixel values to be displayed.
;	If neither SECTOR, RADII, nor ANGLES are set, then the ROI is selected using the
;	cursor to draw the boundary by using DEFROI.
;	If only SECTOR is set then the ROI is an annular sector whose vertex is the sun 
;	center and the radii and angles of the sector are determined interactively.  If
;	RADII is set then the radial values are set to the input and similarly if ANGLES
;	is set.
;	
; EXAMPLE:
;	To find the mass of a CME, where the CME mass image is in
;	'320005.fts', and saving the mass information in 'mass.lst':
;
;		Mass = scc_CME_MASSIMG2TOTAL ('320005.fts',save='mass.lst')
;
;	To use an annular sector, with boundaries defined interactively:
;
;		Mass = scc_CME_MASSIMG2TOTAL (''320005.fts',save='mass.lst',/SECTOR)
;
;	To use an annular sector, with angular boundaries pre-set :
;
;		Mass = scc_CME_MASSIMG2TOTAL ('320005.fts',save='mass.lst',/SECTOR,ANGLES=[250,280])
;
;	or, since the SECTOR keyword is not necessary:
;
;		Mass = scc_CME_MASSIMG2TOTAL ('320005.fts',save='mass.lst',ANGLES=[250,280])
;
; MODIFICATION HISTORY:
; $Log: scc_cme_massimg2total.pro,v $
; Revision 1.12  2012/04/12 20:05:52  nathan
; change default format of printed output to match HT_INFO= format
;
; Revision 1.11  2012/03/05 16:00:53  nathan
; fix spacing format of T1 T2 output
;
; Revision 1.10  2012/03/02 23:31:01  nathan
; Use WCS routines to get coordinates for HI images; fix implementation of
; adjusting display size using PAN= keyword; use different scaling for
; non-mass input
;
; Revision 1.9  2012/03/02 18:54:30  nathan
; convert input position angles to angle from solar north based on input header
;
; Revision 1.8  2012/03/01 16:35:16  nathan
; fix adj for roi_sector call; call get_sun_center with /silent; conditional detstroy wbase
;
; Revision 1.7  2011/05/27 19:22:03  avourlid
; Fixed widget blocking by registering slide_image
;
; Revision 1.6  2011/02/07 18:36:30  nathan
; format dte compatible with previous output of CME mass files
;
; Revision 1.5  2010/12/20 23:06:49  nathan
; return VERSION=Id; allow undefined ha.rsun
;
; Revision 1.4  2010/12/15 20:02:42  nathan
; Generalize for SECCHI or LASCO; make consistent with ed cme_massimg2total_new.pro;
; compute and print also potential, cm_r, cm_a
;
; 	Written by:	RA Howard, NRL, 5/19/97
;	Modified:	RAH, NRL, 3/14/98, comments added
;                       RCC, NRL, 10/31/08 adapted for SECCHI
;
;	@(#)scc_cme_massimg2total.pro	1.0 10/31/08 LASCO IDL LIBRARY
;-
;loadct,0
;tvlct,rrr,ggg,bbb,/get
;ggg[255]=0
;bbb[255]=0
;tvlct,rrr,ggg,bbb

COMMON cme_massimg,base,hbase,calbase,bname, debug

version= '$Id: scc_cme_massimg2total.pro,v 1.12 2012/04/12 20:05:52 nathan Exp $' ; SECCHI IDL LIBRARY

PRINT,'PROCESSING File: ',fn
MREADFITS,fn,ha,a,/dash2underscore
;
;  Check for MASS image
;
tags = TAG_NAMES(ha)
w = WHERE (tags EQ 'BUNIT',nw)
date_obs=ha.date_obs
IF tag_exist(ha,'time_obs') THEN date_obs=date_obs+' '+ha.time_obs
utdt=anytim2utc(date_obs)
dte=utc2str(utdt,/ecs,/trunc)

PRINT,'TIME             ',dte
nx = ha.naxis1
ny = ha.naxis2
;
;  Display the byte scale difference image
;

mn=-0.5e11
mx=1e11
amed=median(a)
IF amed GT mx THEN mx=mx*15

IF (nw EQ 0)  THEN BEGIN
   PRINT,'The units of the input file are not defined.'
   PRINT,'Continuing anyway'
ENDIF ELSE IF (STRMID(ha.bunit,0,5) NE 'GRAMS')  THEN BEGIN
   PRINT,'The units of the input file is not mass:  ',ha.bunit
   PRINT,'Continuing anyway'
   mn=0
   mx=0.9*max(a)
ENDIF
IF (KEYWORD_SET(minscl))  THEN mn=minscl  
IF (KEYWORD_SET(maxscl))  THEN mx=maxscl 
help,mn,mx

IF keyword_set(PAN) THEN adj=1./pan ELSE adj=1.

IF (nx/adj gt 1024) or (ny/adj gt 1024) THEN BEGIN
  SLIDE_IMAGE,rebin(bytscl(a,min=mn,max=mx),nx/adj,ny/adj),xvisible=1024,yvisible=1024,show_full=0,slide_window=sw_a,/register, top_id=wbase
  wset,sw_a
  junk = ''
  read,'Center CME in window. When done press return:',junk
ENDIF ELSE BEGIN
    wnd,30,BYTSCL(a,min=mn,max=mx),1/adj,/tv  
    wait,1
ENDELSE

; mask area so can read time
tv,make_array(177,18,value=128)
XYOUTS,5,5,dte+'UT' ,/device,charsi=1.5
;
;  Adjust scaling
;
IF NOT keyword_set(batch) THEN XLOADCT, /block

;
;  Select a region of interest, compute the size and
;  location of the ROI.
;
IF tag_exist(ha,'RSUN') THEN rsun=ha.rsun ELSE rsun=0.
IF rsun LE 0 THEN rsun=get_solar_radius(ha)
pixrad = rsun/ha.cdelt1
IF ha.cunit1 EQ 'deg' THEN pixrad = pixrad/3600.

IF ha.ctype1 EQ 'HPLN-AZP' THEN BEGIN
; need to do this for all images that are not sun-pointed, such as HI.
    ;--Compute coordinates of cursor location
    wcs=fitshead2wcs(ha)

    cxy=wcs_get_pixel(wcs,[0,0])
    cx=cxy[0]
    cy=cxy[1]
    cunit=wcs.cunit[0]
    IF cunit EQ 'deg' THEN unitf=1.
    arcsc=ABS(wcs.cdelt[0])*3600./unitf

ENDIF ELSE BEGIN

    sun = get_SUN_CENTER(ha, ROLL=sroll, /deg,/silent)
    cx=sun.xcen
    cy=sun.ycen
ENDELSE
;
;  Is a sector desired?
;
help,cx,cy,sroll
ans = 'n'
REPEAT begin
; until ans eq 'n'

IF (KEYWORD_SET(SECTOR) OR KEYWORD_SET(RADII) OR KEYWORD_SET(ANGLES)) THEN BEGIN
;
;  An annular sector is desired.
;  Get the inner and outer radii of the annulus
;  Repeat to make sure they are not equal
;
   roitype = 'SECTOR'
   REPEAT BEGIN
      r1 = 0
      r2 = 0
      IF (KEYWORD_SET(RADII))  THEN BEGIN
         sz = SIZE(radii)
         IF (sz(1) EQ 2)  THEN BEGIN
            r1 = radii(0)
            r2 = radii(1)
         ENDIF
      ENDIF
      IF (r1 EQ r2)  THEN BEGIN
    	PRINT,'SELECT INNER AND OUTER RADII'
    	CURSOR,x,y,3,/device
	x=x*adj
	y=y*adj
    	IF ha.ctype1 EQ 'HPLN-AZP' THEN BEGIN
    	    xy=wcs_get_coord(wcs,[x,y])
    	    IF wcs.projection EQ 'ARC' THEN $
    	    rth=[xy[0],xy[1]+90] ELSE $
    	    rth=hpc2hpr(xy, DEGREES=(unitf LT 3600),/silent)
    ; output always degrees
    	    r1=rth[1]*3600./rsun
	ENDIF ELSE $
            r1 = SQRT((x-cx)^2+(y-cy)^2)/pixrad
    	CURSOR,x,y,3,/device
	x=x*adj
	y=y*adj
    	IF ha.ctype1 EQ 'HPLN-AZP' THEN BEGIN
    	    xy=wcs_get_coord(wcs,[x,y])
    	    IF wcs.projection EQ 'ARC' THEN $
    	    rth=[xy[0],xy[1]+90] ELSE $
    	    rth=hpc2hpr(xy, DEGREES=(unitf LT 3600),/silent)
    ; output always degrees
    	    r2=rth[1]*3600./rsun
	ENDIF ELSE $
    	    r2 = SQRT((x-cx)^2+(y-cy)^2)/pixrad
      ENDIF
      IF (r1 GT r2) THEN BEGIN
         x = r1
         r1 = r2
         r2 = x
      ENDIF
   ENDREP UNTIL (r1 NE r2)
   ar = 0.5*(r1+r2)		; the center distance
   ht = r2-r1			; the height of the annulus
;
;  Get the left and right boundaries of the sector
;  Repeat to make sure they are not equal
;
   REPEAT BEGIN
      t1 = 0
      t2 = 0
      IF (KEYWORD_SET(ANGLES)) THEN BEGIN
         sz = SIZE(angles)
         IF (sz(1) EQ 2)  THEN BEGIN
            t1 = angles(0)
            t2 = angles(1)
         ENDIF
      ENDIF
      IF (t1 EQ t2) THEN BEGIN
        PRINT,'SELECT LEFT AND RIGHT HAND BOUNDARIES (ANGLES)'
        CURSOR,x,y,3,/device
	x=x*adj
	y=y*adj
    	IF ha.ctype1 EQ 'HPLN-AZP' THEN BEGIN
    	    xy=wcs_get_coord(wcs,[x,y])
    	    IF wcs.projection EQ 'ARC' THEN $
    	    rth=[xy[0],xy[1]+90] ELSE $
    	    rth=hpc2hpr(xy, DEGREES=(unitf LT 3600),/silent)
    ; output always degrees
    	    t1=rth[0]
	ENDIF ELSE BEGIN
            t1 = ATAN(x-cx,y-cy)*!radeg
            IF (t1 LT 0)   THEN t1=t1+360.
            IF (t1 GT 360)   THEN t1=t1-360.
            t1a = 360-t1		; convert to position angles
	    t1=t1a+sroll
	ENDELSE
	;
    	CURSOR,x,y,3,/device
    	x=x*adj
	y=y*adj
    	IF ha.ctype1 EQ 'HPLN-AZP' THEN BEGIN
    	    xy=wcs_get_coord(wcs,[x,y])
    	    IF wcs.projection EQ 'ARC' THEN $
    	    rth=[xy[0],xy[1]+90] ELSE $
    	    rth=hpc2hpr(xy, DEGREES=(unitf LT 3600),/silent)
    ; output always degrees
    	    t2=rth[0]
	ENDIF ELSE BEGIN
            t2 = ATAN(x-cx,y-cy)*!radeg
            IF (t2 LT 0)   THEN t2=t2+360.
            IF (t2 GT 360)   THEN t2=t2-360.
            t2a = 360-t2		; convert to position angles
	    t2=t2a+sroll
	ENDELSE
      ENDIF
   ENDREP UNTIL (t1 NE t2)
   ac = 0.5*(t1+t2)		; the central angle
   wt = t1-t2			; the width of the sector
   IF (wt LT 0)  THEN wt=wt+360
;
;  Convert the definition of the annular sector into pixel numbers
;
help,r1,r2,t1,t2,t1a,t2a

    roi =scc_ROI_SECTOR(ha,r1,r2,t1,t2,/draw,adjust=[adj,adj])
    
    help,roi    
                                ; Fudge to make format consistent with fluxrope IDL procs
    ONE2TWO, roi, [nx, ny], cols, rows
    rr = SQRT( (cols - cx)^2 + (rows - cy)^2 )
    ar = ac
    ac = wt
    ht = r1
    wt = r2
                                ;End Fudge
    help,rr,ar,ac,ht,wt
    
ENDIF ELSE BEGIN
;
;  Since a sector is not desired, define the ROI using the cursor
;

   roitype = 'DEFROI'
print,'SELECT ROI'
   roi = DEFROI (nx,ny,/NOFILL)
    one2two, roi, [nx, ny], cols, rows
    rr = SQRT( (cols - cx)^2 + (rows - cy)^2 )
;
;  Calculate the central row and column numbers in radii from sun center
;  Calculate the maximum height and width in radii
;
    ht = max(rr)/pixrad			; the height of the annulus
    wt = min(rr)/pixrad			; the width of the sector
    ac = average(rr)/pixrad		; the central angle
    ar = pixrad
  

ENDELSE
np = n_elements(roi)
;
;  Finally, compute the mass (grams) and print out the mass,
;  median and average values of grams/cm^2
;
GMR_sun = 1.9167e15             ;G=6.67e-8 dyn*cm2/gr2, M=2e33gr,Rsun=6.96e10cm

mass = a[roi]

potential = GMR_sun*TOTAL(mass*(1. - pixrad/rr)) ;ergs

mass = TOTAL(mass)
cm_x = TOTAL(a[roi]*(cols - cx))/mass
cm_y = TOTAL(a[roi]*(rows - cy))/mass
cm_r = SQRT(cm_x^2 + cm_y^2)/pixrad
cm_a = -ATAN(cm_x, cm_y)*!radeg ;CCW from Solar North   -AV@99/06/01
;dte = ha.date_obs+' '+STRMID(ha.time_obs,0,8)

     PRINT,'Date-Obs            Cen-PA  Width Height ROI-Typ     Pts PixRad       Mass    Mass/Pt  Potential    Cm_R    Cm_A'
     PRINT, dte, ac, wt, ht, roitype, np, pixrad, mass, mass/np, potential, cm_r, cm_a, $
     FORMAT='(a19,3f7.2,a8,i8,f7.2,3e11.3,2f8.2)'

    wait,1

IF keyword_set(batch) THEN ans = 'n' $
ELSE read,  'Repeat measurement (y,n)? ', ans

ENDREP UNTIL ( ans EQ 'n')

IF (KEYWORD_SET(SAVE)) THEN BEGIN
  sz = SIZE(save)
  ff = ha.filename
  len = strlen(ff)
  SAVE, file=STRMID(ff, 0, len-4) + '.roi', roi
  IF (sz(1) EQ 7)  THEN file=save ELSE $
    file = PICKFILE(/write,FILTER='*.mas')
  
  ff= FIND_FILE(file)

  OPENW,lu,file,/get_lun,/append

  IF (KEYWORD_SET(HT_INFO)) THEN BEGIN
     IF((SIZE(ff))(0) EQ 0) THEN $ 
     PRINTF, lu,'Date-Obs            Cen-PA  Width Height ROI-Type (R1     R2     T1     T2       Pts) PixRad     Mass    Potential    Cm_R   Cm_A'
     PRINTF, lu, dte, ht_info,roitype, r1, r2, t2, t1, np, pixrad, mass, potential, cm_r, cm_a, $
     FORMAT='(a19,3f7.2,a7,2f7.2,2f8.2,i8,f7.2,2e11.3,2f8.2)'
  ENDIF ELSE BEGIN
     IF((SIZE(ff))(0) EQ 0) THEN $ 
     PRINTF, lu,'Date-Obs            Cen-PA  Width Height ROI-Typ     Pts PixRad       Mass    Mass/Pt  Potential    Cm_R    Cm_A'
     PRINTF, lu, dte, ac, wt, ht, roitype, np, pixrad, mass, mass/np, potential, cm_r, cm_a, $
     FORMAT='(a19,3f7.2,a8,i8,f7.2,3e11.3,2f8.2)'

  ENDELSE
  CLOSE,lu
  FREE_LUN,lu
ENDIF
;
IF datatype(wbase) NE 'UND' THEN widget_control, /destroy, wbase
RETURN,mass
END