pro trace_dem, index, data, basis=basis, emc=emc, A=A, Tavg=Tavg, K_inv=K_inv, $
    emtot=emtot, volume=voume, tau=tau, debug=debug

datasize = size(data)
if datasize(0) ne 3 or datasize(3) ne 3 then message, $
    'First parameter must be a data cube containing three 2D images.'
nx = datasize(1)
ny = datasize(2)
npix = nx*ny

;First, calculate K_inv and A for the chosen (or default) basis.
trace_dem_setup, index,[0.0,0.0,0.0], basis=basis, A=A, K_inv=K_inv, volume=voume, tau=tau

;Convert pixel values to DN/s
p = fltarr(nx,ny,3)
for i=0, 2 do p(*,*,i) = (data(*,*,i)/index(i).sht_mdur)

;Calculate emission measure per basis element and EM averaged temperature
emc = fltarr(nx,ny,3)
Tavg = fltarr(nx,ny)
for x = 0L, nx-1 do begin
for y = 0L, ny-1 do begin
    emc(x,y,*) = K_inv # reform(p(x,y,*),3)
    Tavg(x,y) = total(emc(x,y,*)*tau)/total(emc(x,y,*))
endfor
endfor

emtot = total(emc, 3)

end