function progver, dummy
;+
;NAME:
;	jitter_subs
;PURPOSE:
;	A suite of routines for MDI jitter mode analysis
;HISTORY:
; V4.1	20-Jan-95 (MDM) - Various cosmetic changes
;			- Broke subroutines into JITTER_SUBS.PRO
;			- Allowed sub-menu list to be selected
; V4.2	20-Jan-95 (MDM) - More modifications and adjustments
;			- Added MF_SAMPLE structure to the summary
;			  structure
;			- Backup up 32 samples at missing 1/16th samples
;			  (when 239 instead of 240)
; V4.3	23-Jan-95 (MDM) - Made freq peak intensity plot a log plot
; V4.4	27-Jan-95 (MDM) - Made patch to handle incomplete dataset
;			  during MK_SUMMARY/SORT_ISS
; V4.5	 2-Feb-95 (MDM) - Patched to handle a mixture of ISS_MODE
;			  and image data in a single file.
; V5.0	 6-Feb-95 (MDM) - Modified to pass an assortment of information
;			  by replacing "ss_blocks" with a structure
;			- Pass and display the input file name
;			- Allow background subtraction
; V5.1	 7-Feb-95 (MDM) - Allowed reading/writing summary files on
;			  non-SGI machine (byte swapping)
; V5.11	 8-Feb-95 (MDM) - Modified to print bkg subraction info on
;			  the PLOT_FFT output
;			- Added menu option to change data directory
; V6.00	17-Feb-95 (MDM) - Added Tarbell code to integrate and calculate
;			  overall jitter
;			- Allowed block selection by clicking on image
;			- Allowed masking of spectra by marking bad
;			  frequencies
;			- Allowed totaling of selected ranges
; V6.01	17-Feb-95 (MDM) - Small changes in plotting format (oplot 0 line)
;			- When removing points, remove only the top 3
; V6.10	 1-Mar-95 (MDM) - Various changes and additions
; V6.20  6-Mar-95 (MDM) - Added capability to add title to the plots/displays
; V6.21  6-Mar-95 (MDM) - Allow using previous background for background
;			  subtraction
; V6.22	 7-Mar-95 (MDM) - Modification to not count the same frequency twice
;			  when integrating to get a total (when marking two
;			  regions which overlap)
; V7.00	 6-Feb-96 (MDM) - Modified to work on the flight system
; V7.01	11-Feb-96 (MDM) - Modified to read the raw TFR files
; V7.02	22-Feb-96 (MDM) - Changed all plot reference to "Power" to "Amplitude"
;			- Added new option to plot X error,y error amplitude
;			- Modified the title creation to be shorter
; V7.03	26-Jun-96 (MDM) - Mod to allow resetting of regions to integrate
;			- Added menu option to save spectra to file
; V7.04 27-Jun-96 (MDM) - Corrected derivation of summary saved spectra file name
;			- Converted integrated total to arcsec RMS
; V7.05	 9-Jul-96 (MDM) - Corrected to save "summary" variable in summary file
; V7.06 16-Jul-96 (MDM) - Made xmenu_sel the default.  Use "qxwin=0" to use
;			  wmenu_sel
;			- Made /data14/jitter the default directory
; V7.07	10-Sep-96 (MDM) - Modification to the way that the raw telemetry file is
;			  read to fix a new error state
; V7.08	11-Nov-96 (MDM) - Modified DISP_FFT option to mark integrated blocks
;			  selected
; V7.09 22-Jan-97 (MDM) - Added OUTDIR option to MK_SUMMARY
;			- Added TFR_DIR to main program
; V7.10  4-Sep-97 (RIB) - Possible fix to subscripting problem
;			  Updated default start and stop times in .jitter_do
; V7.11	25-Sep-97 (MDM) - Fix ss=ss_block to be a LONG, not string
;			  (IDL Ver 5.0 broke this!?!?!)
;-

return, 'JITTER Ver 7.11'
end
;---------------------------------------------------------------
function get_blocks, summary


pzt0 = summary.pzt_peakf(0,0)
pzt1 = summary.pzt_peakf(0,1)
pzt2 = summary.pzt_peakf(0,2)


op_cl_arr = [' Opn ', ' Cls ']
op_cl = (summary.mf.mkarr(3) / 16) mod 2
str1 = op_cl_arr(op_cl)
n = n_elements(summary)
ss = indgen(n/4)*4	;oops - hardwired the n16th number
str1(ss) = strmid(str1(ss), 0, 4) + '*'

fmt = '(f7.1)'
str2 = string(pzt0, format=fmt) +string(pzt1, format=fmt) +string(pzt2, format=fmt) 
menu = fmt_tim(summary) + str1 + str2
	


imenu = wmenu_sel(menu)

return, imenu
end
;---------------------------------------------------------------
function get_item, code, summary, gen_info, item, tit2=tit2, get_tit=get_tit, $
		timrange=timrange, all=all, fid=fid
;
;
;
tit_arr =      ['PZT Voltage 0 (Y)', $
		'PZT Voltage 1', $
		'PZT Voltage 2', $
		'PZT Voltage 2-1 (X)', $
		'X Error', $
		'Y Error', $
		'PZT 0 / Y Error', $
		'PZT 2-1 / X Error', $
		' ', $
		' ', $
		'Diode 0', $
		'Diode 1', $
		'Diode 2', $
		'Diode 3']
tit2_arr = strarr(20)
;
tit1 = tit_arr(code)
tit2 = tit2_arr(code)
if (n_elements(item) ne 0) then case item of
    'avg': tit1 = tit1 + ' Average'
    'dev': tit1 = tit1 + ' Standard Deviation'
    else: 
end
if (keyword_set(get_tit)) then return, tit1
;
;----------------------
;
ss_blocks = long(str2arr(gen_info.ss_blocks))
if (ss_blocks(0) lt 0) then ss_blocks = indgen(n_elements(summary))
gen_info.ss_blocks = arr2str(ss_blocks)
ss = long(ss_blocks)
n = n_elements(ss)
timrange = fmt_tim(summary(ss(0))) + ' to ' + fmt_tim(summary(ss(n-1))) + '  (' + strtrim(n,2) + ' blocks)
fid1 = str_replace(ex2fid( anytim2ex(summary(ss(0))), /sec ), '.', '_')
fid2 = str_replace(ex2fid( anytim2ex(summary(ss(n-1))), /sec ), '.', '_')
fid = fid1 + '-' + strmid(fid2, 7, 6)
;
;----------------------
;
if (keyword_set(item)) then begin
    if (code ge 0) and (code le 3) then begin
	str1 = 'pzt_'
	str2 = strtrim(code,2)
    end
    if (code eq 4) or (code eq 5) then begin
	str1 = 'xy_'
	str2 = strtrim((code-4)<1, 2)
    end
    if (code ge 10) then begin
	str1 = 'diode_'
	str2 = strtrim(code-10,2)
    end
    if (keyword_set(all)) then str2 = '*'
    if (strmid(item,0,4) eq 'peak') then str2 = '*,' + str2	;second dimension
    cmd = 'out = summary(ss).' + str1 + item + '(' + str2 + ')'
    ;;print, cmd
    stat = execute(cmd)
    return, out
end
;
;----------------------
;
case code of
    0: out = summary(ss).pzt_fft(*, 0)
    1: out = summary(ss).pzt_fft(*, 1)
    2: out = summary(ss).pzt_fft(*, 2)
    3: out = summary(ss).pzt_fft(*, 2) - summary(ss).pzt_fft(*, 1)
    4: out = summary(ss).xy_fft(*,0)
    5: out = summary(ss).xy_fft(*,1)
    6: out = summary(ss).pzt_fft(*, 0) / summary(ss).xy_fft(*,1)
    7: out = (summary(ss).pzt_fft(*, 2) - summary(ss).pzt_fft(*, 1)) / summary(ss).xy_fft(*,1)
    else: out = -1
endcase
;
return, out
end
;---------------------------------------------------------------
function jitter_files, base_dir, st_dattim, en_dattim, exten, qprint=qprint
;
;
;
ff = file_list(base_dir, exten, /cd)
break_file, ff, dsk_log, dir, filnam, ext
ss = where((filnam ge st_dattim) and (filnam le en_dattim), nss)
if (nss ne 0) then out = filnam(ss) else out = 'NONE AVAILABLE'
if (keyword_set(qprint)) then prstr, out

return, out
end
;---------------------------------------------------------------
pro rd_and_sort, infil, out, dep_tim, fast_tim, dsets=dsets
;
;
;
nfil = n_elements(infil)
ndset_arr = intarr(nfil)
for ifil=0,nfil-1 do begin
    infil0 = infil(ifil)
    read_mdi, infil0, xxx, ndset, /nimages
    ndset_arr(ifil) = ndset
    if (keyword_set(dsets)) then ndset_arr(ifil) = n_elements(dsets)
end
ndset_tot = total(ndset_arr)
;
ii = 0
for ifil=0,nfil-1 do begin
    infil0 = infil(ifil)
    ;
    for idset=0,ndset_arr(ifil)-1 do begin
	print, 'Processing dataset ', idset, ' of ', ndset_arr(ifil), ' for file: ', infil0 
	iidset = idset
	if (keyword_set(dsets)) then iidset = dsets(idset)
	read_mdi, infil0, iidset, data
	sort_iss, data, out0, dep_tim0, fast_tim0
	if (ifil eq 0) and (idset eq 0) then begin
	    out = replicate(out0, ndset_tot)
	    fast_tim = dblarr(ndset_tot * n_elements(fast_tim0))
	    dep_tim  = replicate(dep_tim0(0), ndset_tot * n_elements(dep_tim0))
	end
	out(ii) = out0
	fast_tim(ii*n_elements(fast_tim0)) = fast_tim0
	dep_tim(ii*n_elements(dep_tim0)) = dep_tim0
	ii = ii + 1
    end
end
;
end
;---------------------------------------------------------------
pro sort_iss, data, out, dep_tim, fast_tim
;
;
;qdebug = 1
;qstop = 1
;
mf_samp_ref = {struct1c, marker: 0b, lobt: bytarr(6), sample: fix(0), mkarr: bytarr(15)}
;NOTE:  There is an extra byte in between lobt and sample and one of the 15 "mk"
;       fields is missing.   Type:   print, out.mf.sample        or     print, out.mf.mkarr(3) --should be 2 for 
;	going from open to closed loop
;;mf_samp_ref = {struct1b, marker: 0b, lobt: bytarr(6), sample: fix(0), $
;;	mp75lt: 0b, mkpwlt: 0b, mkltloop: 0b, mxam1pos: bytarr(2), mxam2pos: bytarr(2), $
;;	mklt_ofs: bytarr(2), mklt_gn: bytarr(2), mklt_os: bytarr(3), mkltstat: 0b}
	;;mkarr: bytarr(15)}
dep_tick_ref = {struct2, marker: 0b, lobt: bytarr(6), pzt_arr: bytarr(3), $
	xy_err: bytarr(2), diode: bytarr(4)}
fast_samp_ref = {struct3, marker: 0b, pzt_arr: bytarr(3), xy_err: bytarr(2)}
data_ref0 = {struct4,  dep: dep_tick_ref, fast: replicate(fast_samp_ref, 32), nfast: 0}
;
data_ref = {struct5e, mf: mf_samp_ref, one16th: replicate(data_ref0, 240), n16ths: 0}	;every 15 seconds (per dataset)
;
n = n_elements(data(0,*))
out = replicate(data_ref, n)
;
qdebug = keyword_set(qdebug)
qstop  = keyword_set(qstop)
iout = 0
for idset=0,n-1 do begin
    b = byte(data(*,idset), 0, 63882)
    qdone = 0
    ss = where(b eq 'FE'x)	;marks start of major frame sample
    if (ss(0) eq -1) then qdone = 1		;MDM added 8-Feb-96
    while (not qdone) do begin
	ss = where(b eq 'FE'x)	;marks start of major frame sample
	ii = ss(0)
	reads, b(ii:ii+23), mf_samp_ref
	out(iout).mf	= mf_samp_ref
	if (qdebug) then print, 'MF Sample', idset, iout, mf_samp_ref.lobt
	ii = ii + 24
	;
	qdone2 = 0
	iout2 = 0
	while (not qdone2) do begin
	    i16th = 0
	    temp0 = b( indgen(32)*6 + 16 + ii )
	    nn = (16 + 6*32) + 2	;+2 for "nfast"
	    reads, b(ii:ii+nn-1), data_ref0
	    ;
	    temp = data_ref0.fast.marker
	    ss = where((temp eq '00'x) or (temp eq 'FF'x), nss)		;'00' marks start of DEP header
									;'ff' is fill at end of dataset 
									;     -- where dep start would be
	    if (nss eq 1) then begin
		nfast = ss(0)
	    end else begin
		ss2 = where(temp((ss-1)>0) ne '7C'x)
		if (ss2(0) eq -1) then nfast = ss(0) else nfast = ss(ss2(0))
	    end
	    if ((nss eq 2) and (ss(0) ne -1)) then if (max(temp(ss)) eq 0) then nfast = ss(0)  ;RIB fudge 4-Sep-97
;	    if (nss eq 2) and (max(temp(ss)) eq 0) then nfast = ss(0)	;MDM fudge 10-Sep-96
	    if (ss(0) ne -1) then if (temp(ss(0)) eq 'FF'x) then begin
		nfast = ss(0)
		qdone2 = 1
	    end
	    ;
	    if (qdebug) then print, 'DEP Sample', idset, iout, iout2, nfast, data_ref0.dep.lobt
	    ;
	    ;;out(iout).one16th(iout2) = data_ref0
	    if (nfast gt 0) then begin		;MDM 27-Jan-95 - handle case of incomplete dataset
		out(iout).one16th(iout2).dep		= data_ref0.dep
		out(iout).one16th(iout2).fast(0).pzt_arr	= data_ref0.dep.pzt_arr
		out(iout).one16th(iout2).fast(0).xy_err	= data_ref0.dep.xy_err
		;out(iout).one16th(iout2).fast(1:nfast) 	= data_ref0.fast(0:nfast-1)
		;---- 4-Jan-94 4:00 PM (MDM) - put the fast packets in the right slot - copy fast packet "0" in missing
		out(iout).one16th(iout2).fast(32-nfast:31) 	= data_ref0.fast(0:nfast-1)
		out(iout).one16th(iout2).nfast 		= nfast+1
		nmissing = 32-nfast-1
		for imiss=1,nmissing do out(iout).one16th(iout2).fast(imiss) =  out(iout).one16th(iout2).fast(0) 
	    	iout2 = iout2 + 1
	    end else begin
		qdone2 = 1
	    end
	    ii = ii + 16 + 6L*nfast
	    if (b(ii) eq 'FE'x) then qdone2 = 1		;found start of major frame sample
	    if (qstop) then stop
	end

	out(iout).n16ths = iout2

	;Patch for 239 pieces - 5-Jan-95 2:20 PM (MDM)	
	if (iout2 ne 240) then begin
	    for jjj=iout2,239 do begin
		out(iout).one16th(jjj).fast(*) = out(iout).one16th(iout2-1).fast(31)
		out(iout).one16th(jjj).dep.lobt = out(iout).one16th(iout2-1).dep.lobt	;Added 20-Jan-95 18:28
	    end
	end
	qdone = 1
        iout = iout + 1
    end
end
;
;------------
n = n_elements(out)
nmax_dep  = n*240
nmax_fast = n*240*31
ndep  = total(out.n16ths)
nfast = total(out.one16th.nfast)
;
tai_sec = long(out.one16th.dep.lobt(0:3), 0, nmax_dep)
frac_sec = fix(out.one16th.dep.lobt(4:5), 0, nmax_dep)
frac_sec = unsign(frac_sec, 16) / double(2.^16.)
dt = tai_sec + frac_sec
ss_dep = where(dt ne 0)
ss0 = where(dt eq 0)
if (ss0(0) ne -1) then dt(ss0) = dt(ss0-1)
dep_tim = anytim2ints('1-jan-58', off=double(dt))
;
;-------------
;
fast_tim = (rebin(dt, nmax_fast))

end
;---------------------------------------------------------------
pro calc_summary, out, dep_tim, fast_tim0, n16th, npeak, out2
;
;
;
n = n_elements(out)
						;fast and dirty way to do it
ref_tim = fmt_tim(dep_tim(0))
fast_tim = fast_tim0 - int2secarr(ref_tim, '1-jan-58')

nout = n*240/n16th
n1 = n16th*32		;approximate number of 512 Hz values
n2 = n1/2 + 1
n3 = 240/n16th		;normally = 4
mf_samp_ref = {struct1c, marker: 0b, lobt: bytarr(6), sample: fix(0), mkarr: bytarr(15)}
data_ref = {day: 0, time: 0L, mf: mf_samp_ref, $
		pzt_avg: fltarr(4), pzt_dev: fltarr(4), pzt_fft: complexarr(n2, 3), $
		pzt_peaka: fltarr(npeak,3), pzt_peakp: fltarr(npeak,3), pzt_peakf: fltarr(npeak,3), $
		xy_avg: fltarr(2), xy_dev: fltarr(2), xy_fft: complexarr(n2, 2), $
		xy_peaka: fltarr(npeak,2), xy_peakp: fltarr(npeak,2), xy_peakf: fltarr(npeak,2), $
		diode_avg: fltarr(4), diode_dev: fltarr(4), $
		nfast: 0}
out2 = replicate(data_ref, nout)
;
qdebug = 1
for iout=0L,nout-1 do begin
    iii = iout / n3
    out2(iout).mf	= out(iii).mf

    ss = indgen(n1) + iout*n1
    jjj = out(iii).n16ths
    if (jjj ne 240) and ((iout mod n3) eq (n3-1)) then ss = ss-32*(240-jjj)	;back up 32 samples per 1/16 missing

    for i=0,2 do begin
	do_calc, fast_tim(ss), (out.one16th.fast.pzt_arr(i))(ss), iavg, idev, spec, npeak, peaka, peakp, peakf
	out2(iout).pzt_avg(i)	= iavg
	out2(iout).pzt_dev(i)	= idev
	out2(iout).pzt_fft(*,i)	= spec
	out2(iout).pzt_peaka(*,i) = peaka
	out2(iout).pzt_peakp(*,i) = peakp
	out2(iout).pzt_peakf(*,i) = peakf

    end
    do_calc, fast_tim(ss), (out.one16th.fast.pzt_arr(2))(ss) - (out.one16th.fast.pzt_arr(1))(ss), iavg, idev
    out2(iout).pzt_avg(3)  = iavg
    out2(iout).pzt_dev(3)  = idev
    ;
    for i=0,1 do begin
	do_calc, fast_tim(ss), (out.one16th.fast.xy_err(i))(ss), iavg, idev, spec, npeak, peaka, peakp, peakf
	out2(iout).xy_avg(i)	= iavg
	out2(iout).xy_dev(i)	= idev
	out2(iout).xy_fft(*,i)	= spec
	out2(iout).xy_peaka(*,i) = peaka
	out2(iout).xy_peakp(*,i) = peakp
	out2(iout).xy_peakf(*,i) = peakf
    end
    ;
    ss = indgen(n16th) + iout*n16th
    for i=0,3 do begin
	do_calc, fast_tim(ss), (out.one16th.dep.diode(i))(ss), iavg, idev
	out2(iout).diode_avg(i)	= iavg
	out2(iout).diode_dev(i)	= idev
    end

    out2(iout).time	= dep_tim(ss(n16th/2)).time		;mid time
    out2(iout).day	= dep_tim(ss(n16th/2)).day		;mid time

    if (qdebug) then print, iout, ' of ', nout, ss(0), fmt_tim(out2(iout))
end

end		
;----------------------------
pro do_calc, x0, y, iavg, idev, spec, npeak, peak_amp, peak_phase, peak_freq
;
;
;
n = n_elements(y)
idev = stdev(y, iavg)
;

if (n_params() ge 5) then begin		;only calculate if asked for
    spec = fft(y, 1)
    spec = spec(0:n/2)
    ;;if (min(abs(spec)) eq 0) then stop, 'Minimum of ABS(SPEC) is zero
    ;;spec = abs(spec(1:n/2))
    ;
    yy0 = spec(1:*)
    yy1 = alog10( abs(yy0)>0.1 )
    ;
    nn = 96
    navg = (n/2) / nn
    med_arr = fltarr(navg)
    for i=0,navg-1 do med_arr(i) = median( yy1(i*nn:(i+1)*nn-1) )
    xx1 = (findgen(navg)+0.5)*nn
    xx1 = [0, xx1]
    med_arr = [ total(yy1(0:4))/5, med_arr]
    coeff = poly_fit(xx1, med_arr, 5, yfit)
    xx3 = findgen(n/2)
    yy3 = poly(xx3, coeff)

    ;;plot, yy1
    ;;oplot, xx1, med_arr, thick=3
    ;;oplot, xx3, yy3

    yy4 = yy1 - yy3
    ;;plot, yy4

    ;;npeak = 10
    peak_amp = fltarr(npeak)
    peak_phase = fltarr(npeak)
    peak_freq  = fltarr(npeak)
    ;256 HZ full range / 960 data points ==> 0.266 Hz per point
    mask_width = 7	;mask out 15*0.266 window (4 Hz window, +/- 2 HZ)
    yy5 = yy4
    dmin = min(yy5)
    for ipeak=0,npeak-1 do begin
	dmax = max(yy5, imax)
	;peak_amp(ipeak) = abs(yy0(imax))
	peak_amp(ipeak) = abs(yy0(imax)) - 10.^yy3(imax)
	peak_phase(ipeak) = atan( imaginary(yy0(imax)), float(yy0(imax)) ) * 180/!pi
	peak_freq(ipeak) = imax / (n/2.) * 256.
	ii0 = (imax - mask_width) > 0
	ii1 = (imax + mask_width) < (n_elements(yy5)-1)
	yy5( ii0:ii1 ) = dmin
    end

    ;;plot, yy5
    ;;oplot, peak_freq,fltarr(npeak)+1,psym=2

    ;;yy = smooth(yy0, 10)
    ;;dyy = deriv_arr(yy)
    ;find + to - locations (zero derivative)
    ;;ss = where( (dyy(1:*) le 0) and (dyy ge 0) )
end
;
end
;---------------------------------------------------------------
pro iss_plot1, x, y, n, nx, ny, tit, ref_tim, timerange, gen_info
;
tv2, 800, 600, /land, /init, /already
;
savep = !p
if (n eq 3) then !p.charsize = 1.7 else !p.charsize = 0
!p.multi=[0,nx,ny]
for i=0,n-1 do begin
    y0 = y(i,*,*,*,*)
    utplot, x, y0, ref_tim, psym=psym, tit=tit(i), /ynoz, timerange=timerange, /xstyle
end
;
plottime, 0, .97, gen_info.infil
plottime, 0, 0, progver()
plottime
pprint
!p = savep
end
;----------------------------
pro plot_jit1, out2, gen_info

tv2, 800, 600, /land, /init, /already

save_p = !p
!p.multi=[0,4,4]
!p.psym = 10
!p.charsize = 1.4
for i=0,3 do begin
    ii = strtrim(i,2)
    ii2 = strtrim(i-2,2)

    if (i le 3) then utplot, out2, out2.pzt_avg(i), /ynoz, tit='Average ' + get_item(i,/get_tit)	$
		else !p.multi(0) = !p.multi(0) - 1
    if (i le 3) then utplot, out2, out2.pzt_dev(i), /ynoz, tit='STDDEV ' + get_item(i,/get_tit) 	$
		else !p.multi(0) = !p.multi(0) - 1
    if (i le 1) then utplot, out2, out2.xy_avg(i), /ynoz, tit='Average ' + get_item(i+4,/get_tit)	$
		else utplot, out2, out2.diode_avg(i-2), /ynoz, tit='Diode Average #' +ii2
    if (i le 1) then utplot, out2, out2.xy_dev(i), /ynoz, tit='STDDEV ' + get_item(i+4,/get_tit)	$
		else utplot, out2, out2.diode_avg(i), /ynoz, tit='Diode Average #'+ii	
end
!p = save_p
;
plottime, 0, .97, gen_info.infil, 1.2
plottime, 0, 0, progver()
plottime
pprint
end		
;----------------------------
pro plot_jit4, out2, gen_info, codes, item, qpause=qpause

ss_blocks = long(str2arr(gen_info.ss_blocks))
if (n_elements(ss_blocks) eq 1) and (ss_blocks(0) ne -1) then begin
    tbeep, 3
    print, 'More than one data point required for time plots.  Returning'
    return
end

if (n_elements(item) ne 4) then item = strarr(4) + item
tv2, 800, 600, /land, /init, /already

save_p = !p
!p.multi=[0,2,2]
!p.psym = 10
!p.charsize = 1.4

x = anytim2ints(out2)
for i=0,3 do begin
    y = get_item(codes(i), out2, gen_info, item(i), timrange=timrange)
    tit = get_item(codes(i), out2, gen_info, item(i), /get_tit)
    utplot, x, y, /ynozero, tit=tit, /xstyle
end
;
!p = save_p
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
plottime, 0, .97, gen_info.infil, 1.2
plottime, .02, .02, timrange, 1.2, 90
plottime, 0, 0, progver()
plottime
if (keyword_set(qpause)) then pause
pprint
end
;----------------------------
pro disp_fft, out2, code, gen_info, nn=nn, qpause=qpause
;
;
;  for i=0,2 do disp_fft, out02, 1, i
;  for i=0,1 do disp_fft, out02, 0, i

common disp_fft_blk, x0, y0, nx, ny

y = out2(0).pzt_fft(*,0)
n = n_elements(y)
n2 = n_elements(out2)
xx = (findgen(n)+1)/n * 512. *0.5
;;plot, xx, y

;;tvscl, congrid( out2.pzt_fft(*,0), 512, 512)
;;tvscl, congrid( out2.pzt_fft(0:200,1), 512, 512)

x0 = 70
y0 = 150
nx = 512
ny = 512

if (n_elements(ii) eq 0) then ii = 0
if (n_elements(nn) eq 0) then nn = 200

tv2, 650, 750, /init, /already, win=0
if (!d.name eq 'X') then wshow

tit = get_item(code, /get_tit)
yyy = get_item(code, out2, gen_info, timrange=timrange)
yyy = yyy(1:nn,*)
nblock = n_elements(yyy(0,*))

amp = abs(yyy)
if (gen_info.qbkg) and (code le 5) then for i=0,nblock-1 do amp(*,i) = amp(*,i) - gen_info.bkg_amp(1:*,code)


if (min(amp) le 0) then begin	;MDM patch 10-Jan-95
    ss = where(amp le 0)
    ss2 = where(amp ne 0)
    ;amp(ss) = min(amp(ss2))
    if (not gen_info.qbkg) then amp = amp > 1
    ;;stop
end
amin = min(amp)>1	;MDM added >1 1-Mar-95
amax = max(amp)
if (gen_info.yrange(0) ne gen_info.yrange(1)) then begin
    amin = gen_info.yrange(0)
    amax = gen_info.yrange(1)
end
if (not gen_info.qbkg) then bimg = congrid(  bytscl( alog10(amp)), nx, ny) $
			else bimg = congrid(  bytscl(amp, amin, amax), nx, ny) 
erase
xyouts2, 650/2, 750- 50, tit, align=0.5, siz=2.5
tv2, bimg, x0, y0
ss_blocks = long(str2arr(gen_info.ss_blocks))
time_list = anytim2ints( out2( ss_blocks( indgen(9)/8. * (nblock-1) + 0.5 ) ) )
xyouts2, x0+nx, y0+ny+10, gt_day( time_list(0), /str), align=0.5, siz=1.5
;for i=0,8 do xyouts2, x0+nx+5, y0+i/8.*ny, gt_time( time_list(i), /str )
for i=0,8 do begin
    xyouts2, x0+nx+15, y0+i/8.*ny, gt_time( time_list(i), /str )
    plots2, [x0+nx+1, x0+nx+10], [1,1]*(y0+i/8.*ny), /dev
end
imax = xx(nn)
step = fix(imax/8)

;if (gen_info.ss_integrate ne '') then begin
;    ss_marks = long(str2arr(gen_info.ss_blocks))
;    print, ss_marks
;    for i=0,n_elements(ss_marks)-1 do plots2, x0+[-7,-3], y0+ss_marks(i)/float(n2)*ny, /dev 
;end
if (gen_info.ss_data ne '') then begin
    ss_marks = long(str2arr(gen_info.ss_data))
    print, ss_marks
    for i=0,n_elements(ss_marks)-1 do plots2, x0+[-15,-8], y0+ss_marks(i)/float(n2)*ny, /dev 
end

if (nblock gt 1) then yavg = total( amp, 2) / nblock else yavg = amp
plot, xx(0:nn), yavg, xtit='Frequency (Hz)', ytit='Amplitude', tit = 'Average over Time', $
	position = [x0/650., .08, (x0+nx)/650., .16], /noerase, /xstyle, ytype=(not gen_info.qbkg), psym=10, $
	yrange = [amin, amax]

plottime, 0, .97, gen_info.infil, 1.2
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
if (gen_info.qbkg) and (code le 5) then plottime, 0, .93, 'BKG Subtracted with: ' + gen_info.bkg_info
plottime, 0, 0, progver()
plottime, .02, .02, timrange, 1.2, 90
plottime
if (keyword_set(qpause)) then pause
pprint
end
;----------------------------
pro plot_fft, out2, code, gen_info, nn=nn, qpause=qpause
;
;
;
y = out2(0).pzt_fft(*,0)
n = n_elements(y)
n2 = n_elements(out2)
xx = (findgen(n)+1)/n * 512. *0.5

if (n_elements(nn) eq 0) then nn = 200

tv2, 650, 900, /init, /already, win=0
if (!d.name eq 'X') then wshow

yyy = get_item(code, out2, gen_info, timrange=timrange)
yyy = yyy(1:nn,*)
nblock = n_elements(yyy(0,*))
tit = get_item(code, /get_tit)
;;tit = tit + '!c ' + timrange

amp = abs(yyy)
if (gen_info.qbkg) and (code le 5) then for i=0,nblock-1 do amp(*,i) = amp(*,i) - gen_info.bkg_amp(1:*,code)
phase = atan(yyy)*180/!pi

if (nblock gt 1) then begin
   amp = total(amp, 2) / nblock
   phase = total(phase, 2) / nblock
end

amin = min(amp)>1	;MDM added >1 1-Mar-95
amax = max(amp)
if (gen_info.yrange(0) ne gen_info.yrange(1)) then begin
    amin = gen_info.yrange(0)
    amax = gen_info.yrange(1)
end
!p.multi=[0,1,2]
plot, xx(0:nn), amp,   xtit='Frequency (Hz)', /xstyle, psym=10, ytit='Amplitude', tit=tit, $
		ytype=(not gen_info.qbkg), yrange=[amin, amax]
plot, xx(0:nn), phase, xtit='Frequency (Hz)', /xstyle, psym=10, ytit='Phase', tit=tit

plottime, 0, .97, gen_info.infil, 1.2
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
if (gen_info.qbkg) and (code le 5) then plottime, 0, .93, 'BKG Subtracted with: ' + gen_info.bkg_info
plottime, 0, 0, progver()
plottime, .02, .02, timrange, 1.2, 90
plottime
if (keyword_set(qpause)) then pause
pprint
end
;----------------------------
pro plot_xy_amp, out2, gen_info, nn=nn, qpause=qpause
;
;
;
y = out2(0).pzt_fft(*,0)
n = n_elements(y)
n2 = n_elements(out2)
xx = (findgen(n)+1)/n * 512. *0.5

if (n_elements(nn) eq 0) then nn = 200

tv2, 650, 900, /init, /already, win=0
if (!d.name eq 'X') then wshow

code = 4		;------------------- X error
yyy = get_item(code, out2, gen_info, timrange=timrange)
yyy = yyy(1:nn,*)
nblock = n_elements(yyy(0,*))
tit1 = get_item(code, /get_tit)
amp1 = abs(yyy)
if (gen_info.qbkg) and (code le 5) then for i=0,nblock-1 do amp1(*,i) = amp1(*,i) - gen_info.bkg_amp(1:*,code)
if (nblock gt 1) then amp1 = total(amp1, 2) / nblock

code = 5		;------------------- X error
yyy = get_item(code, out2, gen_info, timrange=timrange)
yyy = yyy(1:nn,*)
nblock = n_elements(yyy(0,*))
tit2 = get_item(code, /get_tit)
amp2 = abs(yyy)
if (gen_info.qbkg) and (code le 5) then for i=0,nblock-1 do amp2(*,i) = amp2(*,i) - gen_info.bkg_amp(1:*,code)
if (nblock gt 1) then amp2 = total(amp2, 2) / nblock

!p.multi=[0,1,2]
amin = min(amp1)>1	;MDM added >1 1-Mar-95
amax = max(amp1)
if (gen_info.yrange(0) ne gen_info.yrange(1)) then begin
    amin = gen_info.yrange(0)
    amax = gen_info.yrange(1)
end
plot, xx(0:nn), amp1,   xtit='Frequency (Hz)', /xstyle, psym=10, ytit='Amplitude', tit=tit1, $
		ytype=(not gen_info.qbkg), yrange=[amin, amax]
;
amin = min(amp2)>1	;MDM added >1 1-Mar-95
amax = max(amp2)
if (gen_info.yrange(0) ne gen_info.yrange(1)) then begin
    amin = gen_info.yrange(0)
    amax = gen_info.yrange(1)
end
plot, xx(0:nn), amp2,   xtit='Frequency (Hz)', /xstyle, psym=10, ytit='Amplitude', tit=tit2, $
		ytype=(not gen_info.qbkg), yrange=[amin, amax]
;

plottime, 0, .97, gen_info.infil, 1.2
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
if (gen_info.qbkg) and (code le 5) then plottime, 0, .93, 'BKG Subtracted with: ' + gen_info.bkg_info
plottime, 0, 0, progver()
plottime, .02, .02, timrange, 1.2, 90
plottime
if (keyword_set(qpause)) then pause
pprint
end
;----------------------------
pro plot4fft, out2, gen_info, nn=nn, qpause=qpause
;
;
;
y = out2(0).pzt_fft(*,0)
n = n_elements(y)
xx = (findgen(n)+1)/n * 512. *0.5
if (n_elements(nn) eq 0) then nn = 200

tv2, 650, 900, /init, /already, win=0
if (!d.name eq 'X') then wshow

codelist = [3,0,4,5]
;
!p.multi = [0,1,4]
!p.charsize = 1.5
for icode=0,3 do begin
    code = codelist(icode)
    yyy = get_item(code, out2, gen_info, timrange=timrange)
    yyy = yyy(1:nn,*)
    nblock = n_elements(yyy(0,*))
    tit = get_item(code, /get_tit)
    ;;tit = tit + '!c ' + timrange

    amp = abs(yyy)
    if (gen_info.qbkg) and (code le 5) then for i=0,nblock-1 do amp(*,i) = amp(*,i) - gen_info.bkg_amp(*,code)

    if (nblock gt 1) then amp = total(amp, 2) / nblock

    plot, xx(0:nn), amp,   xtit='Frequency (Hz)', /xstyle, /ytype, psym=10, ytit='Amplitude', tit=tit
end

plottime, 0, .97, gen_info.infil, 1.2
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
if (gen_info.qbkg) and (code le 5) then plottime, 0, .93, 'BKG Subtracted with: ' + gen_info.bkg_info
plottime, 0, 0, progver()
plottime, .02, .02, timrange, 1.2, 90
plottime
if (keyword_set(qpause)) then pause
pprint
end
;----------------------------
pro plot_transf, out2, gen_info, nn=nn, qpause=qpause
;
;
;
y = out2(0).pzt_fft(*,0)
n = n_elements(y)
n2 = n_elements(out2)
xx = (findgen(n)+1)/n * 512. *0.5

if (n_elements(nn) eq 0) then nn = 200
tv2, 650, 900, /init, /already, win=0
if (!d.name eq 'X') then wshow

codelist = [6,7]
;
!p.multi = [0,1,4]
for icode=0,1 do begin
    code = codelist(icode)
    yyy = get_item(code, out2, gen_info, timrange=timrange)
    yyy = yyy(1:nn,*)
    nblock = n_elements(yyy(0,*))
    tit = get_item(code, /get_tit)
    ;;tit = tit + '!c ' + timrange

    amp = abs(yyy)
    phase = atan(yyy)*180/!pi

    if (nblock gt 1) then begin
       amp = total(amp, 2) / nblock
       phase = total(phase, 2) / nblock
    end

    plot, xx(0:nn), amp,   xtit='Frequency (Hz)', /xstyle, /ytype, psym=10, ytit='Amplitude', tit=tit
    plot, xx(0:nn), phase, xtit='Frequency (Hz)', /xstyle        , psym=10, ytit='Phase', tit=tit
end

plottime, 0, .97, gen_info.infil, 1.2
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
plottime, 0, 0, progver()
plottime, .02, .02, timrange, 1.2, 90
plottime
if (keyword_set(qpause)) then pause
pprint
end
;----------------------------
pro plot_peak, out2, code, gen_info, qpause=qpause
;
;
;
savep = !p
y1 = get_item(code, out2, gen_info, 'peakf')
y2 = get_item(code, out2, gen_info, 'peaka')
y3 = get_item(code, out2, gen_info, 'peakp', timrange=timrange)
ss_blocks = long(str2arr(gen_info.ss_blocks))
x = out2(ss_blocks)
;
n = 3	;3 top peaks
yrange1 = [0, max(y1(0:n-1,*))]
yrange2 = [1, max(y2(0:n-1,*))]
yrange3 = [-200,200]
;
tit = get_item(code, /get_tit)
;;tit = tit + '!c ' + timrange
;
tv2, 650, 900, /init, /already, win=0
if (!d.name eq 'X') then wshow
;
!p.charsize = 1.4
!p.multi=[0,1,3]
utplot, x, y1(0,*), /xstyle, psym=10, ytit='Peak Freq', tit=tit, yrange=yrange1
for i=1,n-1 do outplot, x, y1(i,*), psym=10, linestyle=i
utplot, x, y2(0,*), /xstyle, psym=10, ytit='Peak Amplitude', tit=tit, yrange=yrange2, /ytype
for i=1,n-1 do outplot, x, y2(i,*), psym=10, linestyle=i
utplot, x, y3(0,*), /xstyle, psym=10, ytit='Peak Phase', tit=tit, yrange=yrange3
for i=1,n-1 do outplot, x, y3(i,*), psym=10, linestyle=i

plottime, 0, .97, gen_info.infil, 1.2
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
plottime, 0, 0, progver()
plottime, .02, .02, timrange, 1.2, 90
plottime
if (keyword_set(qpause)) then pause
!p = savep
pprint
end
;------------------------------
pro mk_summary, infil, dsets, outdir=outdir
;
;HISTORY:
;	 1-Mar-95 (MDM) - Added "dsets" optional input
;
;
;
swap_os=['vms','ultrix','OSF']                  ; list of swap machines
chk=where(!version.os eq swap_os,nscount)       ; current OS in noswap list?
qswap = nscount ne 0                            ; assign boolean
;
n16th = 60		;3.75 sec
npeak = 10		;find location of 10 peaks
head = fltarr(20)
;
dpc = 'ffffff30'x
dpc = '6a000000'x
dpc_str = string(dpc, format='(z8.8)')
;
nfil = n_elements(infil)
for ifil=0,nfil-1 do begin
    infil0 = infil(ifil)
    read_mdi, infil0, xxx, ndset0, /nimages
    ;
    break_file, infil0, dsk_log, dir, filnam, ext
    ;
    tfd = file_list(dir, filnam + '*.tfd')
    spawn, 'grep ' + dpc_str + ' ' + arr2str(tfd, delim=' '), result
    if (result(0) eq '') then begin
	print, 'No ' + dpc_str + ' datasets found in ' + tfd
	tbeep, 3
    end else begin
	list_tfi, infil0, info, /quiet			;--- MDM added 2-Feb-95
	ss = where(info.dpc eq dpc, ndset)		;--- MDM added 2-Feb-95
	print, ndset, ' of ', ndset0, ' are FFFFFF30'
	;
	extra = ''
	if (n_elements(dsets) ne 0) then begin
	    ndset = n_elements(dsets)
	    print, 'Manually overriding:  Using', ndset, ' of ', ndset0
	    ss = dsets
	    extra = '_' + string(ss(0), format='(i3.3)')
	end
	;
	;outfil = dsk_log + dir + strmid(filnam, 0, 13) + '.iss_sum'
	if (keyword_set(outdir)) then dir = outdir
	;outfil = dsk_log + dir + filnam + extra + '.iss_sum'
	outfil = concat_dir(dir, filnam + extra + '.iss_sum')
	ssw_file_delete, outfil
	openw, lun, outfil, /get_lun
	head(0) = 1.0		;version
	head(1) = ndset
	head(2) = n16th
	head(3) = ndset * (240/n16th)
	head(4) = npeak
	if (qswap) then dec2sun, head
	writeu, lun, head
	;
	;for idset=0,ndset-1 do begin
	for i=0,ndset-1 do begin
		idset = ss(i)
		print, 'Processing dataset ', idset-ss(0), ' of ', ndset, ' for file: ', infil0 
		idl_speed
		read_mdi, infil0, idset, data
		sort_iss, data, out, dep_tim, fast_tim
		calc_summary, out, dep_tim, fast_tim, n16th, npeak, out2
		if (qswap) then dec2sun, out2
		writeu, lun, out2
		idl_speed, 1
	end
	free_lun, lun
    end
end
;
end
;------------------------------
pro rd_summary, infil, out
;
;
;
swap_os=['vms','ultrix','OSF']                  ; list of swap machines
chk=where(!version.os eq swap_os,nscount)       ; current OS in noswap list?
qswap = nscount ne 0                            ; assign boolean
;
nfil = n_elements(infil)
head = fltarr(20)
;;    head(0) = 1.0		;version
;;    head(1) = ndset
;;    head(2) = n16th
;;    head(4) = npeak
;
nout = 0
for ifil=0,nfil-1 do begin
    infil0 = infil(ifil) + '.iss_sum'
    openr, lun, infil0, /get_lun
    readu, lun, head
    if (qswap) then dec2sun, head
    nout = nout + head(1)*(240/head(2))
    free_lun, lun
end
;
n2 = (head(2) * 32) / 2 + 1
npeak = head(4)
mf_samp_ref = {struct1c, marker: 0b, lobt: bytarr(6), sample: fix(0), mkarr: bytarr(15)}
data_ref = {day: 0, time: 0L, mf: mf_samp_ref, $
		pzt_avg: fltarr(4), pzt_dev: fltarr(4), pzt_fft: complexarr(n2, 3), $
		pzt_peaka: fltarr(npeak,3), pzt_peakp: fltarr(npeak,3), pzt_peakf: fltarr(npeak,3), $
		xy_avg: fltarr(2), xy_dev: fltarr(2), xy_fft: complexarr(n2, 2), $
		xy_peaka: fltarr(npeak,2), xy_peakp: fltarr(npeak,2), xy_peakf: fltarr(npeak,2), $
		diode_avg: fltarr(4), diode_dev: fltarr(4), $
		nfast: 0}
out = replicate(data_ref, nout)
;
ii = 0
for ifil=0,nfil-1 do begin
    infil0 = infil(ifil) + '.iss_sum'
    openr, lun, infil0, /get_lun
    readu, lun, head
    if (qswap) then dec2sun, head
    nout0 = head(1)*(240/head(2))
    out0 = replicate(data_ref, nout0)
    readu, lun, out0
    if (qswap) then dec2sun, out0
    out(ii:ii+nout0-1) = out0
    ii = ii + nout0
    free_lun, lun
end
;
end
;------------------------------
PRO AVE_PSD,SUMMARY,XY_PSD,PZT_PSD,BLOCKS=BL
; AVERAGES POWER SPECTRA OF SELECTED BLOCKS OF AN ISS DATA SUMMARY STRUCTURE
;  INPUTS:
;    SUMMARY = INPUT SUMMARY STRUCTURE
;    BL = ARRAY OF INDICES OF WHICH BLOCKS TO AVERAGE (DEFAULT = ALL)
;  OUTPUTS:
;    XY_PSD  = OUTPUT ARRAY OF X & Y ERROR POWER SPECTRA, FLTARR(961,2)
;    PZT_PSD = OUTPUT ARRAY OF PZT POWER SPECTRA, FLTARR(961,4)
;    PZT_PSD(*,3) = POWER SPECTRUM OF PZT X = PZT 2-1
; TDT  2/15/95

NF = N_ELEMENTS(SUMMARY(0).XY_FFT(*,0))
XY_PSD = 0.*FLTARR(NF,2)
PZT_PSD = 0.*FLTARR(NF,4)
IF (NOT KEYWORD_SET(BL)) THEN BL = INDGEN(N_ELEMENTS(SUMMARY))
NBL = N_ELEMENTS(BL)
XY = SUMMARY(BL).XY_FFT
PZT = SUMMARY(BL).PZT_FFT

FOR I=0,NBL-1 DO BEGIN
	XY_PSD = XY_PSD + ABS(XY(*,*,I))^2
	PZT_PSD(*,0:2) = PZT_PSD(*,0:2) + ABS(PZT(*,*,I))^2
	PZT_PSD(*,3) = PZT_PSD(*,3) + ABS(PZT(*,2,I)-PZT(*,1,I))^2    
ENDFOR

FAC = 1./FLOAT(NBL)
XY_PSD = FAC*XY_PSD
PZT_PSD = FAC*PZT_PSD
END
;------------------------------
FUNCTION PSD_RMS2, PSD, ss_freq
; OUTPUT:  PSD_RMS = RMS DEVIATION OF SIGNAL OVER SPECIFIED BAND
; INPUTS:
;   PSD = POWER SPECTRAL DENSITY FROM 0 TO NYQUIST FREQUENCY INCLUSIVE
; TDT  2/15/95
;HISTORY
; TDT  2/15/95
;	17-Feb-95 (MDM) - Changed to accept a list of indicies to sum

NP = N_ELEMENTS(PSD)
if (n_elements(ss_freq) eq 0) then ss_freq = indgen(np-1) + 1	;skip (0) element 

SUM = TOTAL(PSD(ss_freq))
out = SQRT(2.*SUM)/(2.*(NP-1))
;
return, out
END
;------------------------------
FUNCTION PSD_FIX,PSD,NFIX
; REPLACES SPECIFIED INTERVALS OF BAD DATA IN ARRAY PSD BY INTERPOLATION
; PSD = DATA ARRAY
; NFIX = INTARRAY(2,N) OF N INTERVALS (LOWER & UPPER INDEX) TO REPLACE
; PSD_FIX = NEW CORRECTED ARRAY

T=PSD
TT=[0.,PSD,0.]
NFX = N_ELEMENTS(NFIX(0,*))
FOR I=0,NFX-1 DO BEGIN 
	NL=NFIX(0,I)  &  NU=NFIX(1,I)
	NBAD=1+NU-NL
	XBAD=NL+(1+FINDGEN(NBAD))/FLOAT(NBAD+1)
;	print,xbad
;	print,[TT(0:NL),TT(NU+2:*)]
	T(NL:NU)=INTERPOLATE([TT(0:NL),TT(NU+2:*)],XBAD)
END

RETURN,T
END
;------------------------------
function get_mark_ranges, xx, nss0, x0, x1

wshow
print, 'Mark the starting freqency (right button to reset and abort)'
cursor, x0, y0, /data, /down
if (!err eq 4) then return, -99		;MDM added 26-Jun-96
print, 'Mark the ending frequency'
cursor, x1, y1, /data, /down
ss0 = where((xx ge x0) and (xx le x1), nss0)

return, ss0
end
;------------------------------
pro plot_integ, gen_info, nn=nn, qmask=qmask, qtotal=qtotal, qpzt=qpzt, qstop=qstop
;
;
tv2, 600, 800, /init, /already
;
if (keyword_set(qpzt)) then begin
    tits0 = get_item([3,0],/get_tit)
    y = gen_info.pzt_amp
end else begin
    tits0 = get_item([4,5], /get_tit)
    y = gen_info.xy_amp
end
!p.multi = [0,1,2]
yrange = gen_info.yrange
if (total(yrange) eq 0) and (max(abs(y)) gt 10000) then yrange = [-500, 1500]
ymark = (yrange(1)-yrange(0))*.96 + yrange(0)
;
n = n_elements(y(*,0))
xx = (findgen(n)+1)/n * 512. *0.5
if (n_elements(nn) eq 0) then nn = 960
xrange = [0, xx(nn)]
;
if (gen_info.ss_mask ne '') then ss_mask = long(str2arr(gen_info.ss_mask)) else ss_mask = -1
if (gen_info.ss_integrate ne '') then ss_total = long(str2arr(gen_info.ss_integrate)) else ss_total = -1

for i=0,1 do begin
    tit = tits0(i)
    if (ss_total(0) ne -1) then tit = tit + '  Integrated total: ' + string(0.05*PSD_RMS2(y(*,i), ss_total)) + ' arcsec RMS'
    yy = y(1:*,i)
    if (gen_info.smooth ge 2) then yy = smooth(yy, gen_info.smooth)
    plot, xx(1:*), yy, xrange=xrange, yrange=yrange, psym=10, /xstyle, tit=tit, xtit='Frequence (Hz)'	;TODO - ytit='Power or Amp'
    if (ss_mask(0) ne -1) then oplot, xx(ss_mask), yy(ss_mask), psym=2
    if (ss_total(0) ne -1) then oplot, xx(ss_total), fltarr(1000)+ymark, psym=3
    oplot, !x.crange, [0,0]
end
if (gen_info.quse_bkg and (gen_info.ss_bkg ne '')) then begin
    out_bkg = 'Bkg  Blocks:' + strcompress(gen_info.ss_bkg,/remove)
    if (gen_info.infil ne gen_info.bkg_infil) then out_bkg = out_bkg + ' (' + gen_info.bkg_infil + ')'
end else begin
    out_bkg = 'Bkg  Blocks: None used'
end

plottime, 0, .97, gen_info.infil, 1.2
plottime, 0.5, .98, gen_info.title, 1.4, align=0.5
plottime, 0, 0, progver()
plottime, .02, .02, 'Data Blocks:' + strcompress(gen_info.ss_data,/remove), 1., 90
plottime, .04, .02, out_bkg, 1., 90
if (gen_info.smooth ge 2) then plottime, .06, .02, 'Smooth:' + strtrim(gen_info.smooth,2), 1., 90
plottime
if (keyword_set(qstop)) then stop
if (keyword_set(qpause)) then pause
pprint

if (keyword_set(qmask)) then begin
    ss = -1
    qdone = 0
    while (not qdone) do begin
	wshow
	print, 'Mark the freqency to remove with the left button.  Use center button to exit'
	cursor, x0, y0, /data, /down
	button = !err
	dummy = min(abs(xx-x0), imin)
	ss0 = indgen(15)-7 + imin	;look at 15 points around where marked for the peak
	;;ss0 = get_mark_ranges(xx, nss0, x0, x1)
	;;print, 'Masking ', nss0<3, ' points between ', x0, ' and ', x1
	qdone = button eq 2
	if (not qdone) then begin
	    if (keyword_set(qpzt)) then y = gen_info.pzt_amp $
				else y = gen_info.xy_amp
	    temp = abs(y(*,0)) + abs(y(*,1))
	    dummy = max(temp(ss0), itemp)
	    ss00 = [ss0(itemp)-1, ss0(itemp), ss0(itemp)+1] >0<(n_elements(xx)-1)
	    ;;ss00 = [min(ss0), max(ss0)]
	    for i=0,1 do gen_info.xy_amp(*,i) = psd_fix(y(*,i), ss00)
	    if (ss(0) eq -1) then ss = ss00 else ss = [ss, ss00]
	    gen_info.ss_mask = arr2str(ss)
	    plot_integ, gen_info, nn=nn, qpzt=qpzt
	end
	;;yesnox, 'Mark another range? ', ans, 'Yes'
	;;qdone = 1-ans
    end
end
;
if (keyword_set(qtotal)) then begin
    ss = -1
    qdone = 0
    while (not qdone) do begin
	ss0 = get_mark_ranges(xx, nss0, x0, x1)
	if (ss0(0) eq -1) then begin
	    tbeep, 3
	    print, 'No values selected - try again'
	end else begin
	    if (ss0(0) eq -99) then begin
	        print, 'Resetting all selected ranges'
	        ss0 = -1
 	    end
	    if (ss0(0) ne -1) then print, 'Selecting ', nss0, ' points between ', x0, ' and ', x1
	    if (ss(0) eq -1) then ss = ss0 else ss = [ss, ss0]
	    ss = ss( uniq(ss, sort(ss)) )		;MDM added 7-Mar-95
	    gen_info.ss_integrate = arr2str(ss)
	    plot_integ, gen_info, nn=nn, qpzt=qpzt
	end
	yesnox, 'Mark another range? ', ans, 'Yes'
	qdone = 1-ans
    end
end
;
end
;------------------------------
pro jitt_save_gen, gen_info, summary, outdir
;
;
;
filnam = strmid(gen_info.infil, 0, 13)
ff = file_list(outdir, filnam+'*.genx', file=file)
if (ff(0) eq '') then last_ver = 0 else last_ver = max(fix(strmid(file, 15, 3)))
outfil = concat_dir(outdir, filnam+'_V' + string(last_ver+1, format='(i3.3)'))
savegen, gen_info, summary, file=outfil
;;savegen, gen_info, file=outfil
;
end