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
;	*** CHanged to work with TRACE ***
;			- Didn't update the program version number or name
; V8.00	 9-Jan-97 (MDM) - Updated prog ver and name
;			- Larger main menu listing
;			- Allowed toggling from PZT_Avg to PZT_Err
;			- New labels
;			- Made integrated peak totals unitless (DN)
;			- Changed summary generation to be B-C (1-2) instead of 2-1
; V8.01	13-Jan-97 (MDM) - Modified MK_SUMMARY to not work on files over 2 meg
; V8.02	24-Jan-97 (MDM) - Patch for case where there are 21 blocks in one sec
; V8.03	27-Jan-97 (MDM) - Fixed max freq plotting
;			- Plot 3 PZT when plotting integ spec
; V8.04	28-Jan-97 (MDM) - Fix OPLOT of single point
;			- Fixed displaying GT image when PZT error is selected
;			- Don't plot phase over time, plot amp in linear scale
;			- On peak finding, use median of first 5 points, not avg
;			- Use "XYESNO" instead of "YESNOX"
; V8.05	 3-Apr-97 (MDM) - Allow B-C option during error plots
; V9.00	10-Feb-98 (MDM) - Modified to work with new TRACE reformatted datafiles
; V9.01	 7-Apr-98 (MDM) - Modified to use TRACE_J1_DIR (not TRACE_I1_DIR)
;			- Fixed 8 hr correction which was needed for ground
;			  testing (because of DEP)
;-

return, 'TJITTER Ver 9.01'
end
;---------------------------------------------------------------
function get_blocks, summary, qxwin


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
	


if (qxwin) then imenu = xmenu_sel(menu) $
	   else 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
;
common tjitter_blk, which_pzt
;
;
tit_arr =      ['PZT Average A (Z)', $
		'PZT Average B', $
		'PZT Average C', $
		'PZT Average B-C (X)', $
		'GT X Error', $
		'GT Z Error', $
		'PZT A / Z Error', $
		'PZT B-C / X Error', $
		' ', $
		' ', $
		'Diode 0', $
		'Diode 1', $
		'Diode 2', $
		'Diode 3']

if (which_pzt) then tit_arr = str_replace(tit_arr, 'Average', 'Error')

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 = 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
	if (which_pzt) then begin
	    str1 = 'perr_'
	    str2 = strtrim(code,2)
	end else begin
	    str1 = 'pzt_'
	    str2 = strtrim(code,2)
	end
    end
    if (code eq 4) or (code eq 5) then begin
	str1 = 'gerr_'
	str2 = strtrim((code-4)<1, 2)
    end
    if (code ge 10) and (code le 13) 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 + which_pzt*20 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(*, 1) - summary(ss).pzt_fft(*, 2)
    4: out = summary(ss).gerr_fft(*,0)
    5: out = summary(ss).gerr_fft(*,1)
    6: out = summary(ss).pzt_fft(*, 0) / summary(ss).gerr_fft(*,1)
    7: out = (summary(ss).pzt_fft(*, 1) - summary(ss).pzt_fft(*, 2)) / summary(ss).gerr_fft(*,0)
    20: out = summary(ss).perr_fft(*, 0)
    21: out = summary(ss).perr_fft(*, 1)
    22: out = summary(ss).perr_fft(*, 2)
    23: out = summary(ss).perr_fft(*, 1) - summary(ss).perr_fft(*, 2)

    24: out = summary(ss).gerr_fft(*,0)
    25: out = summary(ss).gerr_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, /last_dot
;;ss = where((filnam ge st_dattim) and (filnam le en_dattim), nss)
;;if (nss ne 0) then out = filnam(ss) else out = 'NONE AVAILABLE'
out = sel_filetimes(fid2ex(st_dattim), fid2ex(en_dattim), filnam)
if (out(0) eq '') then 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
;
;
;
data = rfits(infil)
sort_iss, data, out, dep_tim, fast_tim
;
end
;---------------------------------------------------------------
pro sort_iss, data, out, dep_tim, fast_tim
;
;
qdebug2 = 0
qdebug = 0
qstop = 0
;
mf_samp_ref = {struct1x, marker: bytarr(2), lobt: bytarr(6), sample: fix(0), mkarr: bytarr(14)}
;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)}
dep_tick_ref = {struct2z, marker: 0b, lobt: bytarr(6), diode: bytarr(4), wedge: bytarr(2), iss_reg: 0b}
;;fast_samp_ref = {struct3, marker: 0b, pzt_arr: bytarr(3), xy_err: bytarr(2)}
fast_samp_ref = {struct3z, marker: 0b, pzt_arr: bytarr(3), gt_err: bytarr(2), pzt_err: bytarr(3)}
;;data_ref0 = {struct4,  dep: dep_tick_ref, fast: replicate(fast_samp_ref, 10), nfast: 0}
data_ref0 = {struct4z,  dep: dep_tick_ref, fast: replicate(fast_samp_ref, 10), nfast: 0}
;
data_ref = {struct5za, mf: mf_samp_ref, one16th: replicate(data_ref0, 20), 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 = bytarr(2048)
    b(0) = data(*,idset)	;embed in pad
    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		;;+23	;temp fix to bypass junk
	;
	qdone2 = 0
	iout2 = 0
	while (not qdone2) do begin
	    i16th = 0
	    temp0 = b( indgen(10)*9 + 14 + ii )		;temp fix (+1)
	    nn = (14 + 9*10 + 1) + 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 'FD'x), nss)		;'00' marks start of DEP header
	    ss = where((temp eq '00'x) or (temp ge 90) or (temp eq 'FD'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 (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
	    ;
	    jj = ii + 14  + 9L*nfast
	    ;if (temp(nfast) lt 90) and (temp(nfast) ne 0) and ((b(jj-1) eq 'FD'x) or (b(jj+1) eq 'FD'x)) then nfast = nfast - 1		;MDM patch 23-Dec-96
	    if (b(jj) eq 90) then nfast = nfast + 1
	    if (b(jj-1-9) eq 'FD'x) then nfast = nfast - 1		;MDM patch 23-Dec-96
	    if (b(jj+1-9) eq 'FD'x) then nfast = nfast - 1		;MDM patch 23-Dec-96
	    if (b(jj-9) eq '5a'x) and (b(jj-9+1) eq 0) then begin
		print, 'Funny end'
		nfast = nfast - 1		;MDM patch 23-Dec-96
		qdone2 = 1
	    end
	    ;
	    if (qdebug) then print, 'DEP Sample', idset, iout, iout2, nfast, ii, data_ref0.dep.lobt
	    if (qdebug2) then print, b(ii:ii+13), format='(30z3)'
	    if (qdebug2) then print, b(ii+14:ii+nn-1), format='(9z3)'
	    ;
	    ;;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(10-nfast:9) 	= data_ref0.fast(0:nfast-1)
		out(iout).one16th(iout2).nfast 			= nfast+1
		;;;nmissing = 10-nfast-1
		;;;for imiss=1,nmissing do out(iout).one16th(iout2).fast(imiss) =  out(iout).one16th(iout2).fast(0) 
		nmissing = 10-nfast
		for imiss=0,nmissing do out(iout).one16th(iout2).fast(imiss) = out(iout).one16th(iout2).fast(10-nfast) 
	    	iout2 = iout2 + 1
	    end else begin
		qdone2 = 1
	    end
	    ii = ii + 14  + 9L*nfast
	    ss3 = where(b(ii:*) eq 'FD'x, nss)		;should be zero (sometimes 1)
	    if (ss3(0) gt 0) then begin
		;;if (ss3(0) gt 5) then ss3 = -1		;23-Dec-96 (MDM) Temp patch
		if (ss3(0) gt 5) then stop
		if (qdebug) then print, 'Re-sync on FD.  ii = ', ii, ' # shifted = ', ss3(0)
		ii = ii + ss3(0)
	    end
	    if (b(ii) eq 'FE'x) or (b(ii) eq '00'x) then qdone2 = 1		;found start of major frame sample
	    if (qstop) then stop
	    if (iout2 eq 20) then qdone2 = 1	;MDM 24-Jan-97/27-Jan-97 patch
	end

	if (keyword_set(qstop2)) then stop

	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*20
nmax_fast = n*20*10
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
if (dt(0) lt 9.4210560e+08) then dt = dt - 8*60.*60	;patch for wrong dep time
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('24-May-68', off=double(dt))
;
;-------------
;
fast_tim = (rebin(dt, nmax_fast))

end
;---------------------------------------------------------------
pro calc_summary, out, dep_tim, fast_tim0, n16th, npeak, out2
;
;
qdebug = 1
;
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*20/n16th
n1 = n16th*10		;approximate number of 200 Hz values
n2 = n1/2 + 1
n3 = n16th/20		;normally = 4
;;mf_samp_ref = {struct1c, marker: 0b, lobt: bytarr(6), sample: fix(0), mkarr: bytarr(15)}
mf_samp_ref = {struct1x}
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), $
		perr_avg: fltarr(3), perr_dev: fltarr(3), perr_fft: complexarr(n2, 3), $
		perr_peaka: fltarr(npeak,3), perr_peakp: fltarr(npeak,3), perr_peakf: fltarr(npeak,3), $
		gerr_avg: fltarr(2), gerr_dev: fltarr(2), gerr_fft: complexarr(n2, 2), $
		gerr_peaka: fltarr(npeak,2), gerr_peakp: fltarr(npeak,2), gerr_peakf: fltarr(npeak,2), $
		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 20) 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(1))(ss) - (out.one16th.fast.pzt_arr(2))(ss), iavg, idev
    out2(iout).pzt_avg(3)  = iavg
    out2(iout).pzt_dev(3)  = idev
    ;
    for i=0,2 do begin
	do_calc, fast_tim(ss), (out.one16th.fast.pzt_err(i))(ss), iavg, idev, spec, npeak, peaka, peakp, peakf
	out2(iout).perr_avg(i)	= iavg
	out2(iout).perr_dev(i)	= idev
	out2(iout).perr_fft(*,i)   = spec
	out2(iout).perr_peaka(*,i) = peaka
	out2(iout).perr_peakp(*,i) = peakp
	out2(iout).perr_peakf(*,i) = peakf
    end
    ;
    for i=0,1 do begin
	do_calc, fast_tim(ss), (out.one16th.fast.gt_err(i))(ss), iavg, idev, spec, npeak, peaka, peakp, peakf
	out2(iout).gerr_avg(i)	= iavg
	out2(iout).gerr_dev(i)	= idev
	out2(iout).gerr_fft(*,i)   = spec
	out2(iout).gerr_peaka(*,i) = peaka
	out2(iout).gerr_peakp(*,i) = peakp
	out2(iout).gerr_peakf(*,i) = peakf
    end

    out2(iout).time	= dep_tim(ss(n1/2)/10).time		;mid time
    out2(iout).day	= dep_tim(ss(n1/2)/10).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 = 40	;how many points to bin (96 for MDI)
    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]	;MDM changed 28-Jan-97
    med_arr = [ median( yy1(0:4)), 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.) * 100.
	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
;
if (data_type(gen_info) eq 8) then 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.gerr_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.gerr_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 = 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 * 200. *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 = 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 * 200. *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 * 200. *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 * 200. *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 * 200. *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 = 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, y2(0,*), /xstyle, psym=10, ytit='Peak Amplitude', tit=tit, yrange=yrange2
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
;
;	mk_summary, file_list('/td02/log/rawimage/961221','*jitter*.fits')
;
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 = 20*4		;4 sec
npeak = 10		;find location of 10 peaks
head = fltarr(20)
;
nfil = n_elements(infil)
for ifil=0,nfil-1 do begin
    infil0 = infil(ifil)
    print, 'Now PRocessing: ', infil0
    ;
    break_file, infil0, dsk_log, dir, filnam, ext
    ;
    ndp = 1
    qnew = strmid(filnam, 0, 3) eq 'trj'
    if (qnew) then begin
	img = rfits(infil0, h=h, /nodata)
	ndp = sxpar(h, 'EXT_NROW')	;number of data products
    end

    for i=0,ndp-1 do begin
	if (qnew) then mxf_read_data, infil0, i, index, data	$ ;read_trace, infil0, i, index, data $
		else data = rfits(infil0)
	ndset = n_elements(data(0,*))
	;
	;outfil = dsk_log + dir + filnam + '.iss_sum'
	;outfil = concat_dir('/td10/jitter', filnam + '.iss_sum')
	if (n_elements(outdir) eq 0) then outdir = concat_dir(getenv('GSE_DATA_DIR3'), 'jitter')
	if (qnew) then outfil = concat_dir(outdir, filnam +ext + '_'+ string(i+1,format='(i3.3)') + '.iss_sum') $
		else outfil = concat_dir(outdir, filnam + '.iss_sum')
	file_delete, outfil
	openw, lun, outfil, /get_lun
	head(0) = 1.0		;version
	head(1) = ndset
	head(2) = n16th
	head(3) = ndset * (64*20/n16th)
	head(4) = npeak
	if (qswap) then dec2sun, head
	writeu, lun, head
	;
	idl_speed
	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
	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 + fix(head(1)*(20/head(2)))
    free_lun, lun
end
;
n2 = (head(2) * 10) / 2 + 1
npeak = head(4)
mf_samp_ref = {struct1x, marker: bytarr(2), lobt: bytarr(6), sample: fix(0), mkarr: bytarr(14)}
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), $
		perr_avg: fltarr(3), perr_dev: fltarr(3), perr_fft: complexarr(n2, 3), $
		perr_peaka: fltarr(npeak,3), perr_peakp: fltarr(npeak,3), perr_peakf: fltarr(npeak,3), $
		gerr_avg: fltarr(2), gerr_dev: fltarr(2), gerr_fft: complexarr(n2, 2), $
		gerr_peaka: fltarr(npeak,2), gerr_peakp: fltarr(npeak,2), gerr_peakf: fltarr(npeak,2), $
		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 = fix(head(1)*(20/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,PERR_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
;	27-Jan-97 (MDM) - Added PERR_PSD

NF = N_ELEMENTS(SUMMARY(0).GERR_FFT(*,0))
XY_PSD = 0.*FLTARR(NF,2)
PZT_PSD = 0.*FLTARR(NF,4)
PERR_PSD = 0.*FLTARR(NF,3)
IF (NOT KEYWORD_SET(BL)) THEN BL = INDGEN(N_ELEMENTS(SUMMARY))
NBL = N_ELEMENTS(BL)
XY = SUMMARY(BL).GERR_FFT
PZT = SUMMARY(BL).PZT_FFT
PERR = SUMMARY(BL).PERR_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    
	PERR_PSD = PERR_PSD + ABS(PERR(*,*,I))^2
ENDFOR

FAC = 1./FLOAT(NBL)
XY_PSD = FAC*XY_PSD
PZT_PSD = FAC*PZT_PSD
PERR_PSD = FAC*PERR_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, group=group
;
;
common tjitter_blk, which_pzt
tv2, 600, 800, /init, /already
;
if (n_elements(qpzt) eq 0) then qpzt = 0
case qpzt of
    0: begin
	    tits0 = get_item([4,5], /get_tit)
	    y = gen_info.gerr_amp
	end
    1: begin
	    which_pzt = 0
	    tits0 = get_item([0,1,2],/get_tit)
	    y = gen_info.pzt_amp
	end
    2: begin
	    which_pzt = 1
	    tits0 = get_item([0,1,2],/get_tit)
	    y = gen_info.perr_amp
	end
endcase
;
np = n_elements(tits0)
!p.multi = [0,1,np]
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 * 200. *0.5
if (n_elements(nn) eq 0) then nn = 400
xrange = [0, xx(nn)]
;
if (gen_info.ss_mask ne '') then ss_mask = [str2arr(gen_info.ss_mask)] else ss_mask = -1
if (gen_info.ss_integrate ne '') then ss_total = [str2arr(gen_info.ss_integrate)] else ss_total = -1

for i=0,np-1 do begin
    tit = tits0(i)
    if (ss_total(0) ne -1) then tit = tit + '  Integrated total: ' + string(PSD_RMS2(y(*,i), ss_total)) + ' DN 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
	    case qpzt of
		0: y = gen_info.gerr_amp
		1: y = gen_info.pzt_amp
		2: y = gen_info.perr_amp
	    endcase
	    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,np-1 do gen_info.gerr_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'
	xyesno, 'Mark another range? ', ans, group=group
	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