;+
; NAME:
;	RemoteView_CMEDensity
; PURPOSE:
;	Adds CME density to remoteview
; CATEGORY:
; CALLING SEQUENCE:
	FUNCTION RemoteView_CMEDensity, TT0, TT, nDim, rEdge, CMEmin, CMEmax,	$
		rot_cme=rot_cme, degrees=degrees, cme_smooth=cme_smooth, cme_average=cme_average
; INPUTS:
;	TT0			array[1]; type: time structure
;					time assigned to the CME density matrix read from file
;					Usually the same as the 3D solar wind density matrix with
;					which the CME density is merged.
;	TT			array[1]; type: time structure
;					time at which CME is displayed; the CME matrix expands
;					at a constant speed of V=400 km/s from time TT0 to TT
;					(usually TT will be a couple of days later than TT0)
;
;	The heliospheric matrix is defined by:
;
;	nDim		array[3]; type: integer
;					# longitudinal, latitudinal and radial grid points
;	rEdge		array[2]; type: float
;					inner and outer radius of radial grid (AU)
; OPTIONAL INPUTS:
;	rot_cme=rot_cme scalar or array[2]
;			extra rotations applied to the CME
; OUTPUTS:
;	CMEmin		scalar; type: float
;	CMEmax		scalar; type: float
;					min and max normalized density in CME matrix
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	EulerRotate, gridgen, TimeOp, TimeGet, TimeUnit, flt_read, TimeGet, ToRadians
; PROCEDURE:
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

; Rotation 1653   2443226.338378906       1977    MAR     23    82.838380
; Rotation 1654   2443253.616699219       1977    APR     20   110.116700

rpm = ToRadians(degrees=degrees)

dr  = 0.05								; Bin size of CME matrix (AU)
Sun = [9.5,9.5,0]						; Position of sun in matrix

status = flt_read(filepath(root=getenv('SSW_SMEI_UCSD'),subdir='image','catmay24.txt'), CME)

CME = reform(CME, 20,20,10) ;< 500

r  = gridgen([20,20,10], range=[-9.5,9.5, -9.5,9.5, 0,9]*dr,/multid)

CME = CME*total(r*r,1)				; (r^2)*n = constant in expansion

CMEmin = min(CME,max=CMEmax)			; Return arguments only

dCME = max(smooth(CME < 500,9),r)

dCME = ArrayLocation(r,size=size(CME))	; Location of maximum density

dCME = cv_coord(from_rect=dCME-Sun, /to_sphere)

message, /info, 'Maximum density at distance of'+strcompress(dCME[2])

JD = TimeOp(/subtract, TT[0], TT0, TimeUnit(/days)); Time difference in days

V = 4.									; CME velocity (100 km/s)
V = V/!sun.au*0.086400					; CME velocity (AU/day)

; The CME density matrix is defined in the heliographic coordinate system at time TT0.
; Needed is the CME at time TT in the heliographic coordinate system at time TT.
; A CME location at [lng,lat,rad] at time TT was at location
; (lng+PSun*JD,lat,rad-V*JD) on time TT0

; Heliographic longitude, latitude, radius grid at time TT. These are the grid
; locations of the remote view (NOT the grid in which the CME is defined).
; gridgen returns an array with dimension [3,nDim[0]*nDim[1]*nDim[2]]

r = gridgen(nDim, range=[2*!pi*[0,1], !pi/2*[-1,1], rEdge])

message, /info, '    CME matrix at '+TimeGet(TT0,/ymd)
message, /info, 'Display matrix at '+TimeGet(TT ,/ymd)

; Transform longitude to heliographic longitude at TT0
; Correct radial distance for CME expansion. This converts the
; remote view locations to locations in CME volume.

r[0,*] =  r[0,*]+JD/!sun.siderealp*2*!pi
r[2,*] = (r[2,*]-JD*V) > 0.

; The following Euler angles control the locations of the
; CME in its original grid. If no 'rot_cme' is specified then
; the CME will be put along the x-axis.

CASE n_elements(rot_cme) OF
1   : A = [rot_cme*rpm,0.0,0.0]
2   : A = [rot_cme*rpm,    0.0]
ELSE: A = [0.0        ,0.0,0.0]
ENDCASE

A = A+[0.,dCME[1],-dCME[0]]		; Rotate CME to x-axis

r = EulerRotate(A, r)

r = cv_coord(from_sphere=r, /to_rect)	; Convert to x,y,z coordinates

r = r/dr+Sun#replicate(1.,n_elements(r)); Convert to CME bin width of 0.05 AU

; Interpolate on original density matrix at TT0
; At TT0 the CME is located between r=.098 and r=.398 AU

; CME is a density in electrons/cm^-3, not yet corrected for expansion.
; The expanded density is CME*(r(2,*)/rr)^2.
; In addition the returned density should be normalized to 1 AU.
; The expanded normalized density is CME*r(2,*)^2

IF n_elements(cme_smooth ) NE 0 THEN BEGIN
	message, /info, 'smooth CME density by'+strcompress(cme_smooth)
	CME = smooth(CME, cme_smooth)
ENDIF

IF n_elements(cme_average) NE 0 THEN BEGIN
	cme_ratio = cme_average/mean( CME[ where(CME GT 0) ] )
	message, /info, 'multiply CME density by factor '+strcompress(cme_ratio)
	CME = cme_ratio*CME	
ENDIF

CME = interpolate(CME,r[0,*],r[1,*],r[2,*], missing=0.)
CME = reform(CME, nDim, /overwrite)

RETURN, CME  &  END