pro mk_sdm, wid
;
;+
;NAME:
;	mk_sdm
;PURPOSE:
;	To create the weekly average SXT dark image database
;HISTORY:
;	Written Mar-95 by M.Morrison
;V1.1	11-Jul-95 (MDM) - Modified to write the output to /1d3/yohkoh/sdm
;V1.11	12-Jul-95 (MDM) - Modified to correct minor bug
;V1.12			- Fixed bug when no short exposure is available
;-
;
;
start_time = systime(1)
progverno = 1.120*1000
iver = fix(progverno/1000.)
;
indir = '/ydb/sdml'
;datadir = '/yd3/morrison/sdm'
;outdir = '/yd3/morrison/sdm'
;datadir = [data_paths(), '/0d1/yohkoh_b/MOjob/sdc', '/1d1/yohkoh_b/MOjob/sdc']
datadir = [data_paths(), '/0d1/yohkoh_b/mojob/sdc', '/1d1/yohkoh_b/mojob/sdc']
;outdir = '/yd18/sdm'
outdir = '/1d3/yohkoh/sdm'
;
tarr = weekid2ex(wid)
last_week_tim = anytim2ints(tarr, off=-86400)
last_week = strmid(anytim2weekid(last_week_tim, /str), 0, 5)
;
infil           = concat_dir(indir, 'sdml' + wid + '.genx')
outfil          = concat_dir(outdir, 'sdm' + wid       + 'a.' + string(iver,format='(i2.2)') )
outfil_lastweek = concat_dir(outdir, 'sdm' + last_week + 'a.' + string(iver,format='(i2.2)') )
;
if (not file_exist(infil)) then begin
    print, 'MK_SDM: SDML file not found: ' + infil
    return
end
;
restgen, list, file=infil
nimg = max(list.sdm_imgnum)+1
;
recomp_list = ''
for iout=0,nimg-1 do begin
    ss = where(list.sdm_imgnum eq iout, nss)
    nss_out = [nss, 0, 0]
    list1 = list(ss)
    res = list1(0).res
    dummy = min(abs(list1.tim2week_cen), iclosest)
    ;
    avg_slope = total(list1.slope) / nss
    avg_int   = total(list1.intercept) / nss
    ;
    params = [1.9883, -7.1686, 10.1224, -6.9061, 2.2587, -0.28179]
    xxx = alog10( list1.fms )
    orb_corr = 10.^poly(xxx, params)
    ;
    tot_time_weight = total(list1.time_weight)
    if (res eq 0) then begin
	corn = [list1.corn]
	ucorn = corn(uniq(corn, sort(corn)))
	ncorn = n_elements(ucorn)
	tot_time_weight = fltarr(ncorn)
	for j=0,ncorn-1 do begin
	    ss2 = where(list1.corn eq ucorn(j), nss2)
	    tot_time_weight(j) = total(list1(ss2).time_weight)
	    nss_out(j) = nss2
	end
    end
    ;
    nx = 1024 / (2^res)
    out = fltarr(nx, nx)
    if (res ne 0) then wedge = (lindgen(nx, nx)/nx) / float(nx-1) $
		else wedge = (lindgen(nx, 512)/nx) / float(512-1)
    ;
    dpe = list1(0).dpe
    if (dpe gt 13) then begin
	ss3 = -1
	if (file_exist(outfil)) then begin
	    rd_roadmap, outfil, rmap_out
	    ss3 = where((gt_dpe(rmap_out) le 13) and (gt_res(rmap_out) eq res))
	    if (ss3(0) ne -1) then rd_xda, outfil, ss3, short_index, short
	end
	if (ss3(0) eq -1) then if (file_exist(outfil_lastweek)) then begin
	    print, 'MK_SDM: Checking ' + outfil_lastweek + ' for a short exposure
	    rd_roadmap, outfil_lastweek, rmap_out
	    ss3 = where((gt_dpe(rmap_out) le 13) and (gt_res(rmap_out) eq res))
	    if (ss3(0) ne -1) then rd_xda, outfil_lastweek, ss3, short_index, short
	end
	if (ss3(0) eq -1) then begin
	    short = out		;work around to not crash - bad situation - needs fixing
	    print, 'MK_SDM: No short exposure to use for orbit correction.  None to be used
	    tbeep, 3
	end
    end
    ;
    for i=0,nss-1 do begin
	list0 = list1(i)

	if (i eq 0) then begin
            str = string(list0.res, list0.dpe, nss, format="('----------- Res:', i1, ' DPE:',i2, "+$
                                        "' #Img:', i3, ' ----------')")
            print, str
	end

        str = fmt_tim(list0) + string(list0.tim2week_cen/86400., list0.time_weight, list0.avg, list0.dev, $
                                                format='(4f8.3)') + '  ' + string(list0.st$filename)
        print, str


	infil_raw = string(list0.st$filename)
	infil0 = file_list(datadir, infil_raw)
	infil0 = infil0(0)
	if (infil0 eq '') then begin
	    infil0 = file_list(datadir, infil_raw+'.Z')
	    infil0 = infil0(0)
	    if (infil0 eq '') then stop
	    file_uncompress, infil0
	    infil0 = str_replace(infil0, '.Z', '')
	    recomp_list = [recomp_list, infil0]
	end
	rd_roadmap, infil0, rmap
	idset = tim2dset(rmap, list0, del=del)
	if (del(0) ne 0) then stop, 'Time match delta is not zero' 
	rd_xda, infil0, idset, index, data
	data = sxt_decomp(data)
	;
	;---- Ajust the slope by the orbit factor
	data = (data-list0.intercept)/orb_corr(i) + list0.intercept
	;  xx=findgen(nx) & yy=total(data,1)/nx & coeff=poly_fit(xx,yy,1)
	;
	;---- Ajust the integrated signal by the orbit factor
	icorn = list0.corn
	if (dpe gt 13) then begin
	    if (res ne 0) then data = short + (data-short)/orb_corr(i) $
		else data = short(*,icorn:icorn+511) + (data-short(*,icorn:icorn+511))/orb_corr(i)
	end
	;
	;---- Build up the average
	if (res ne 0) then begin
	    out = out + data * list0.time_weight/tot_time_weight
	end else begin
	    ss3 = where(ucorn eq icorn)
	    ss3 = ss3(0)
	    out(*,icorn:icorn+511) = out(*,icorn:icorn+511) + data * list0.time_weight/tot_time_weight(ss3)
	end
	;
	if (i eq iclosest) then begin
	    index_out = index
	    his_index, /enable
	    his_index, index_out
	    his_index, index_out, 0, 'Q_COMPOSITE', progverno
	    his_index, index_out, 0, 'DAY_COMPOS', nss_out		;save number of images
	    his_index, index_out, 0, 'TIME_COMPOS', [avg_slope, avg_int, 0]*1000
	end
    end
    ;
    if (res eq 0) then if (max(ucorn) eq 512) then begin
	;---- Shift the second half ramp up
	xx=findgen(nx)
	yy=total(out,1)/nx
	coeff1 = poly_fit(xx(10:511),yy(10:511),1)
	coeff2 = poly_fit(xx(512:*), yy(512:*) ,1)
	offset = coeff1(0) - coeff2(0)
	out(*,512:*) = out(*,512:*) + offset	
    end
    ;
    sav_sda, outfil, index_out, out, append=(iout ne 0), /compress
end
;
print, '--------------- Finished ... Recompressing files now -------------'
for i=1,n_elements(recomp_list)-1 do file_compress, recomp_list(i)
end_time = systime(1)
run_time = (end_time-start_time)/60.
print, 'MK_SDM Took ', run_time, ' minutes for week: ' + wid
end
;--------------------------------------------------------
;for i= 1,10 do mk_sdm, '95_'+string(i,format='(i2.2)')		;95_09 data and later not available - need to get off MO disk
;for i= 1,53 do mk_sdm, '94_'+string(i,format='(i2.2)')
;for i=23,53 do mk_sdm, '94_'+string(i,format='(i2.2)')
;for i= 1,53 do mk_sdm, '93_'+string(i,format='(i2.2)')
 for i=14,53 do mk_sdm, '93_'+string(i,format='(i2.2)')
 for i= 1,53 do mk_sdm, '92_'+string(i,format='(i2.2)')
 for i=35,53 do mk_sdm, '91_'+string(i,format='(i2.2)')
end