pro trace_tmap, index, data, T_guess=T_guess, tmap=tmap, em=em, chisq=cs_new, $
    debug=debug, maxit=maxit, dT=dT, ALPHA=ALPHA, BETA=BETA, GAMMA=GAMMA

;Volume expansion factors for 3 fundamental simplex operations
if n_elements(ALPHA) eq 0 then ALPHA = 1.0  ;reflection
if n_elements(BETA) eq 0 then BETA = 0.5    ;contraction
if n_elements(GAMMA) eq 0 then GAMMA = 2.0  ;expansion

if n_elements(maxit) eq 0 then maxit=10

datasize = size(data)
if datasize(0) ne 3 then message, $
    'First parameter must be a data cube (array of 2 or more 2D images).'
nx = datasize(1)
ny = datasize(2)
n = datasize(3) ;number of images



;Construct inital guess using input value of T_guess, if defined.
if n_elements(T_guess) eq 0 then T_guess=1.1e6  ;default temperature
if n_elements(T_guess) eq 1 then T_guess = fltarr(nx,ny) > T_guess
if n_elements(T_guess) ne nx*ny then message, $
    'T_guess must be either a scalar or an array with the same size as data(*,*,0).'
T_old = T_guess

delta_p = sqrt(float(data)/12.0)>1.0  ;photon counting error in DN for each pixel

;Create a trivial 1D simplex, i.e. two points where chisquared is evaluated.
if not keyword_set(dT) then dT = 5e4 ;initial size of 1D simplex
tmap = T_old+dT
tt_badness, data, index, datasize, delta_p, T_old, em, cs_old
tt_badness, data, index, datasize, delta_p, tmap, em, cs_new

;==== Here begins the main loop, where optimization is done. ====
for i = 0L, maxit do begin
    ;Keep 'em sorted so that the new one has lowest chisquared.
    ss = where(cs_new ge cs_old, count)
    if count ne 0 then begin
        matt_20_16, cs_new(ss), cs_old(ss)
        matt_20_16, tmap(ss), T_old(ss)
    endif
    
    ;Now start looking for a replacement for the worst point (T_old)
    
    ;(1) reflect from T_old through tmap
    T_guess = tmap + ALPHA * (tmap - T_old)
    tt_badness, data, index, datasize, delta_p, T_guess, em, cs_guess
    
    ;Assess the results
    ss = where(cs_guess le cs_new, count)   ;best case
    ss1 = where(cs_guess ge cs_old, count1)  ;worst case
    
    ;(2) if cs_guess le cs_new, amplify that reflection
    if count ne 0 then T_guess(ss) = tmap(ss) + GAMMA * (tmap(ss) - T_old(ss))
    ;*****gotta do something with ss here for data, index, delta_p!!
    ;tt_badness, data, index, datasize, delta_p, T_guess(ss), em, cs_guess(ss)
    
    ;(3) else shrink the reflection if it made cs_guess worse than cs_old
    if count1 ne 0 then T_guess(ss1) = tmap(ss1) + BETA * (tmap(ss1) - T_old(ss1))
    ;*****gotta do something with ss here for data, index, delta_p!!
    ;tt_badness, data, index, datasize, delta_p, T_guess(ss1), em, cs_guess(ss1)
    
    ;Since it was hard to vectorize in (2) and (3),
    tt_badness, data, index, datasize, delta_p, T_guess, em, cs_guess
    
    ;Move "new" into "old" and "guess" into "new"
    T_old = tmap
    cs_old = cs_new
    tmap = T_guess
    cs_new = cs_guess
endfor



end


;****** Auxiliary procedures *******

;Switch the values of a pair of variables