function CMEDensity, TT, nLng, nLat, nR, rDat,	CMEmin, CMEmax
;+
; NAME:
;	CMEDensity
; PURPOSE:
; CATEGORY:
; CALLING SEQUENCE:
; INPUTS:
;	TT			array[1]; type: time structure
;
;	The heliospheric matrix is defined by:
;	nLng
;	nLat
;	nR			scalar; type: int
;					# radial grid points
;	rDat		array[2]; type: float
;					range of radial grid (AU)
; OUTPUTS:
;	CMEmin		scalar; type: float
;	CMEmax		scalar; type: float
;					min an max normalized density in CME matrix
; CALLS:
;	EulerRotate, flt3dgen, TimeGet, TimeUnit, flt_read, TimeOp
; PROCEDURE:
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-
on_error, 2							; On error, return to caller

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

TT0 = TimeSet(yr=1977, doy=82.83-1.); Time assigned to CME matrix

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('DAT'), subdir='misc','catmay24.txt'), CME )

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

r = flt3dgen([20,20,10], norm=[-9.5,9.5, -9.5,9.5, 0,9]*dr,/threed)
CME = total(r*r,1)*CME			; r^2*n	= constant in expansion

CMEmin = min(CME,max=CMEmax)

dCME = max(smooth(CME<500,9),r)
;dCME = [(r mod 400) mod 20, (r mod 400)/20, r/400]	; Location of maximum density
dCME = ArrayLocation(r,size=size(CME))				; Location of maximum density

dCME = (cv_coord(from_rect=dCME-Sun, /to_sphere))[0:1]

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 TT

; Heliographic longitude, latitude, radius grid at time TT
; flt3dgen returns an array with dimension [3,nLng*nLat*nR]

r	= flt3dgen([nLng,nLat,nR], norm=[0,2*!pi, -!pi/2,!pi/2, rDat])

; Transform longitude to heliographic longitude at TT0
; Correct radial distance for CME expansion

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

; Rotate CME to x-axis

r[0:1,*] = EulerRotate([0.,dCME[1],-dCME[0]], r[0:1,*])

r = cv_coord(from_sphere=r, /to_rect)	; Convert to x,y,z coordinates
										; Convert to CME bin width of 0.05 AU
r = r/dr+Sun#replicate(1.,n_elements(r))

; 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

CME = Interpolate(CME,r[0,*],r[1,*],r[2,*], missing=0.)

return, reform(CME,nLng,nLat,nR)  &  end