PRO Hxt_mcimg, index, data, index_out, data_out, p, Outfile=outfile, $
prompt_outfile=prompt_outfile, append_outfile=append_outfile, $
channels=channels, Nolook=nolook, Nonorm=nonorm, Bsmooth=bsmooth, $
Img_trange=img_trange, Bck_trange=bck_trange, alt_tsel=alt_tsel, $
Dt_resolution=dt_resolution, Accum_cnts=accum_cnts, $
Ch_accum=ch_accum, Noback=noback, Sys_err=sys_err, $
New_mod_patterns=new_mod_patterns, Xy=xy, Expfac=expfac, $
N_pix=n_pix, Patt_shape=patt_shape, Ph_correct=ph_correct, $
Gamma0=gamma0, Chilim=chilim, Lambda_limit=lambda_limit, $
lambda_max=lambda_max, delta_max=delta_max, iter_max=iter_max, $
Qlook=qlook, Bin=bin, Smooth=smooth, $
Pixons=pixons, Snr=snr, Guess=guess, Quiet=quiet, $
Verbose=verbose, Continued=continued, Partan=partan, $
Resolution=resolution, Pixon_sizes=pixon_sizes, $
pixonmap=pixonmap, imagedata=imagedata, $
chisqr=chisqr, imagesig=imagesig, sensitivity=sensitivity, $
Notty=notty, Lo2hi_guess=lo2hi_guess, mc_trials=mc_trials, $
poisson=poisson, Serr_only=serr_only, Norm_error=norm_error, $
No_serr=no_serr, patt_file=patt_file, grause_file=grause_file, $
St_times=st_times, Int_times=int_times, $
sato_patterns=sato_patterns, calc_btot=calc_btot
;
COMMON Mod_patterns, sc_paras, ictrl, hxi_index0
grdnum = 64
ch_ch = [' LO ',' M1 ',' M2 ',' HI ']
; Do we calculate new mod. patterns??
IF(N_ELEMENTS(p) EQ 0) THEN yn = 0 ELSE yn = 1
IF(KEYWORD_SET(new_mod_patterns)) THEN yn = 0
IF(KEYWORD_SET(xy)) THEN BEGIN
yn = 0
xy_p = xy
ENDIF ELSE xy_p = 0
IF(KEYWORD_SET(expfac)) THEN BEGIN
yn = 0
expf = expfac
ENDIF ELSE expf = 0
IF(KEYWORD_SET(n_pix)) THEN BEGIN
yn = 0
pxlnum = n_pix
ENDIF ELSE pxlnum = 64
IF(KEYWORD_SET(patt_shape)) THEN BEGIN
yn = 0
pshape = patt_shape
ENDIF ELSE pshape = 0
IF(KEYWORD_SET(ph_correct)) THEN BEGIN
yn = 0
phc = ph_correct
ENDIF ELSE phc = 0
IF(KEYWORD_SET(patt_file)) THEN BEGIN
yn = 0
pattfile = patt_file
ENDIF ELSE pattfile = 0
IF(KEYWORD_SET(grause_file)) THEN BEGIN
yn = 0
grfile = grause_file
ENDIF ELSE grfile = 0
IF(yn EQ 0) THEN BEGIN
initpatt = 1
; Read in the default parameters for the SC pattern generation
patt_default, hxi_index0, ph_correct = phc, patt_shape = pshape
pattern_maker, hxi_index0, p, sc_paras, xy = xy_p, $
expfac = expf, n_pix = pxlnum, date = index(0), patt_file = pattfile, $
grause_file = grfile, sato_patterns = sato_patterns
ENDIF ELSE BEGIN
initpatt = 0
pxlnum = N_ELEMENTS(p(*, 0, 0))
patt_default, hxi_index0, ph_correct = phc, patt_shape = pshape
hxi_index0.shape_sav(*) = pxlnum ;gotta have it, jmm 3-jan-1997
ENDELSE
pix2 = pxlnum^2
; Read in the default parameters for the image synthesis
img_default, hxi_index0, sys_err = sys_err, gamma0 = gamma0, chilim = chilim
IF(KEYWORD_SET(nolook) OR KEYWORD_SET(notty)) THEN ictrl = 0 ELSE ictrl = 1
; You need an initial version of hxi_index to pass into the accumulation routine
hxi_index = hxi_index0
;ok, get the time intervals for background and image, and accumulate counts
Hxt_mcaccum, index, data, index1, hxi_index, counts, sigma, bcounts, $
bsigma, trials = mc_trials, img_trange = img_trange, bck_trange = bck_trange, $
accum_cnts = accum_cnts, ch_accum = ch_accum, serr_only = serr_only, $
norm_error = norm_error, no_serr=no_serr, st_times = st_times, $
int_times = int_times, noback = noback
; subtract background
IF(KEYWORD_SET(noback)) THEN zb = 1 ELSE zb = 0 ;probably unneccesary, jmm, 6/22/95
hxtsbt_mult, hxi_index, counts, bcounts, fobs, sig, b_zero = zb, $
csigma = sigma, bsigma = bsigma
; make the images, for each channel and interval
n_intv = N_ELEMENTS(index1)
IF(N_ELEMENTS(channels) GT 0) THEN chans = channels ELSE chans = [0, 1, 2, 3]
nch = N_ELEMENTS(chans)
nch1 = nch-1
hxi_index = hxi_index(chans, *)
fobs = fobs(*, chans, *) & sig = sig(*, chans, *)
IF(KEYWORD_SET(xxx)) THEN BEGIN ;you need to add any sys. error here
sig = sqrt(sig^2+(hxi_index0.syserr*fobs)^2)
ENDIF
fobs = fobs > 0.0 ;jmm 22-sep-95, this got lost when GETVAR was removed from makimg
sig = sig > 1.0 ;jmm 22-sep-95, this got lost when GETVAR was removed from makimg
ndset = n_intv*nch
;If btot is to be calculated, do it here
IF(KEYWORD_SET(calc_btot)) THEN BEGIN
print, 'calculating btot'
btot0 = hxt_btot(fobs, p(*, *, *, chans))
ENDIF
; Do we output the data? New stuff, we'll output only if the
; outfile keyword is set, and append only if the append
; keyword is set. 22-Jun-95, jmm
IF(KEYWORD_SET(outfile)) THEN BEGIN
outfil = outfile
append = 0
ENDIF ELSE IF(KEYWORD_SET(prompt_outfile)) THEN BEGIN
;prompt for filename
t7 = anytim2ex(index1(0))
t7s = strcompress(string(t7), /remove_all)
lt_10 = where(t7 LT 10)
IF(lt_10(0) NE -1) THEN t7s(lt_10) = '0'+t7s(lt_10)
IF(t7(3) LT 100) THEN t7s(3) = '0'+t7s(3)
exten_wvq = t7s(6)+t7s(5)+t7s(4)+'_'+t7s(0)+t7s(1)+t7s(2)+'.'+t7s(3)
def_filnam = 'hxi' + exten_wvq
input, 'ENTER THE OUTPUT FILE NAME', outfil, def_filnam
append = 0
ENDIF ELSE outfil = 'none'
outfil = strlowcase(outfil)
IF (file_exist(outfil)) THEN BEGIN
print, 'FILE ', outfil, ' EXISTS'
IF(KEYWORD_SET(append_outfile)) THEN BEGIN
print, 'IT WILL BE APPENDED'
append = 1
ENDIF ELSE BEGIN
print, 'IT WILL BE OVERWRITTEN, TO APPEND SET /APPEND_OUTFILE'
append = 0
ENDELSE
ENDIF
;output arrays
data_out = fltarr(pxlnum, pxlnum, ndset)
index_out = replicate(str_merge(index1(0), hxi_index(0)), ndset)
;Pixons and qlook are back, jmm 6/22/95
IF(KEYWORD_SET(qlook)) THEN BEGIN
IF(KEYWORD_SET(smooth)) THEN alpha = smooth ELSE alpha = 0.10
IF(KEYWORD_SET(bin)) THEN bnn = bin ELSE bnn = 3
FOR j = nch1, 0, -1 DO BEGIN
print, ch_ch(chans(j)), 'CHANNEL '
makimg0, ictrl, reform(hxi_index(j, *)), index1, reform(fobs(*, j, *)), $
reform(sig(*, j, *)), p(*, *, *, chans(j)), index_j, data_j, $
alpha, bin = bnn, quiet = quiet, verbose = verbose
FOR k = 0, n_intv-1 DO BEGIN
dset = k*nch+j
str_load, index_out, index = dset, index_j(k)
data_out(*, *, dset) = data_j(*, *, k)
ENDFOR
ENDFOR
IF(KEYWORD_SET(bsmooth)) THEN BEGIN
qsmooth = 3
FOR j = 0, ndset-1 DO data_out(*, *, j) = smooth(data_out(*, *, j), qsmooth)
ENDIF
ENDIF ELSE IF(KEYWORD_SET(pixons)) THEN BEGIN
; The following tags are irrelevant for the pixon image reconstruction
hxi_index(*).lambda = 0
hxi_index(*).chilim = 0
hxi_index(*).laminc = 0
hxi_index(*).gamma0 = 0
hxi_index(*).delta = 0
; This is still a mess...
pixonmap = fltarr(pxlnum, pxlnum, ndset)
IF(nch GT 1) THEN chanorder = sort(chans) $ ; go from lo to hi energy
ELSE chanorder = [0] ;sort doesn't work on scalars, jmm 22-jun-95
FOR jch = 0, nch1 DO BEGIN
j = chanorder(jch)
IF(chans(j) EQ 3 AND hxi_index(0).iway GE 2) THEN $
hgain = 1.10 ELSE hgain = 0.91
FOR k = 0, n_intv-1 DO BEGIN
c = hxi_index(0).area*hgain
hxi_index(j, k).cnts_p_cm2 = total(fobs(0:15, j, k))/c
ENDFOR
print, ch_ch(chans(j)), 'CHANNEL '
IF KEYWORD_SET(lo2hi_guess) AND jch GT 0 THEN BEGIN
guess = last_pixon_image
ENDIF
IF(KEYWORD_SET(calc_btot)) THEN btot_j = reform(btot0(j, *)) ELSE btot_j = 0 hxtpix, ictrl, reform(hxi_index(j, *)), index1, reform(fobs(*, j, *)), $
reform(sig(*, j, *)), p, index_j, data_j, quiet = quiet, verbose = verbose, $
initpatt = initpatt, continued = continued, partan = partan, $
resolution = resolution, snr = snr, guess = guess, $
n_pix = long(pxlnum), pixon_sizes = pixon_sizes, notty = notty, $
poisson = poisson, chisqr = chisqr, pixonerr=error_j , $
pixonmap = pixonmap_j, sensitivity = sensitivity, btot_in = btot_j
IF KEYWORD_SET(lo2hi_guess) AND (jch LT nch1) THEN last_pixon_image = data_j
FOR k = 0, n_intv-1 DO BEGIN
IF KEYWORD_SET(lo2hi_guess) AND (jch LT nch1) THEN BEGIN
; Flip N/S and renormalize to the next filter. The normalization is only approximate.
last_pixon_image(*, *, k) = rotate(last_pixon_image(*, *, k), 7)* $
total(fobs(*, chanorder(jch+1), k))/total(fobs(*, j, k))
ENDIF
dset = k*nch+j
str_load, index_out, index = dset, index_j(k)
data_out(*, *, dset) = data_j(*, *, k)
pixonmap(*, *, dset) = pixonmap_j(*, *, k)
ENDFOR
ENDFOR
IF(KEYWORD_SET(bsmooth)) THEN BEGIN
qsmooth = 3
FOR j = 0, ndset-1 DO data_out(*, *, j) = smooth(data_out(*, *, j), qsmooth)
ENDIF
ENDIF ELSE BEGIN ;MEM is the default
;for each interval, create images for each channel
FOR k = 0, n_intv-1 DO BEGIN
index1_k = index1(k)
print, 'INTERVAL #', k
print, 'INTERVAL START TIME:', fmt_tim(index1_k, /msec)
print, 'ACCUMULATION TIME:', hxi_index(0, k).actim/10.0
IF(KEYWORD_SET(lambda_limit)) THEN lambda_cr = lambda_limit
FOR j = nch1, 0, -1 DO BEGIN
print, ch_ch(chans(j)), 'CHANNEL '
fobs_jk = fobs(*, j, k) & sig_jk = sig(*, j, k)
hxi_jk = hxi_index(j, k)
IF(KEYWORD_SET(calc_btot)) THEN btot_jk = btot0(j, k) ELSE btot_jk = 0 IF(KEYWORD_SET(lambda_limit)) THEN BEGIN
IF(lambda_limit EQ 1) AND (j EQ nch1) THEN BEGIN
makimg, ictrl, hxi_jk, index1_k, fobs_jk, sig_jk, $
p(*, *, *, chans(j)), index_jk, data_jk, bsmooth = bsmooth, $
delta_max = delta_max, iter_max = iter_max, $
chi_limit = chilim, btot_in = btot_jk
; obtain the value of lambda from the highest channel, to use later
lambda_cr = hxi_jk.lambda
IF(lambda_cr LT 8) THEN lambda_cr = 8
ENDIF ELSE makimg_l, ictrl, hxi_jk, index1_k, fobs_jk, sig_jk, $
p(*, *, *, chans(j)), index_jk, data_jk, $
lambda_max = lambda_cr, bsmooth = bsmooth, $
delta_max = delta_max, iter_max = iter_max, btot_in = btot_jk
ENDIF ELSE makimg, ictrl, hxi_jk, index1_k, fobs_jk, sig_jk, $
p(*, *, *, chans(j)), index_jk, data_jk, bsmooth = bsmooth, $
delta_max = delta_max, iter_max = iter_max, $
lambda_max = lambda_max, chi_limit = chilim, btot_in = btot_jk
; Load the output arrays
dset = k*nch+j
str_load, index_out, index = dset, index_jk
data_out(*, *, dset) = data_jk
ENDFOR
ENDFOR
ENDELSE
; Normalize the data, if not asked not to
IF(NOT KEYWORD_SET(nonorm)) THEN data_out = hxi_normap(temporary(data_out))
; Output the data, if necessary
IF(outfil NE 'none') THEN sav_hxi, outfil, index_out, data_out, append = append
RETURN
END