pro plot_fft, timarr, data, out, ff2, freq, period=period, smoo=smoo, ytit=ytit, hc=hc, step=step, $
		tit=tit, xrange=xrange, yrange=yrange, xtype=xtype, ytype=ytype, $
		mark_freq=mark_freq, mark_period=mark_period, $
		random=random
;+
;NAME:
;	plot_fft
;PURPOSE:
;	To plot facilitate generic plotting of an fft
;SAMPLE CALLILNG SEQUENCE:
;	plot_fft, timarr, data
;
;	quick_hkplot,'23-apr 20:00','24-apr 24:00',style=3,info=info
;	plot_fft, info.mdi.daytime, info.mdi.value(*,0),ytit='Oven temperature (mtopts4)',/xtype,xrange=[1e-4,.1],step=7.5
;
;	plot_fft, timarr, data, out, ff2, freq
;INPUTS:
;	timarr	- time array in seconds (OR) time structure array
;	data	- The data array
;OUTPUTS:
;	out	- abs(ff2(1:n2/2))^2
;	ff2	- The FFT results (complex)
;	freq	- The frequency array (which goes with "out")
;OPTIONAL KEYWORD INPUT:
;	period	- Plot versus the period, not the frequency
;	smoo	- Smooth the FFT spectra by this value
;	ytit	- The raw data plot Y axis label
;	hc	- Make a hardcopy
;	step	- The sample frequency (default is to find the minimum
;		  separation)
;	tit	- Main title
;	xrange	- The FFT frequency (or period) range to plot
;	yrange	- The FFT amplitude to plot
;	xtype	- xtype=0 will make X axis linear (default is log)
;	ytype	- ytype=0 will make Y axis linear (default is log)
;	mark_freq- A list of frequencies to mark on the plot
;	mark_period- A list of periods to mark on the plot
;	random	- Instead of interpolating across missing data, insert
;		  random samples of the input data to fill the gap
;HISTORY:
;	Written 25-Apr-96 by M.Morrison
;V2.1	25-Apr-96 (MDM) - Added /HC option
;			- changed FFT plot symbol to 10
;V2.2	25-Apr-96 (MDM) - Added passing in list of periods or freq to mark
;			- Made LOG/LOG the default
;V2.3	16-Dec-96 (MDM) - Made the FINDGENs into DINDGENs
;V2.4	17-Dec-96 (MDM) - Added /RANDOM option
;			- Added documentation header
;-
;
progver = 'PLOT_FFT  Ver 2.4'
;
if (keyword_set(hc)) then begin
    set_plot,'ps
    device,/land
end
;
if (data_type(timarr) eq 8) then xx1 = int2secarr(timarr) $
				else xx1 = timarr
n1 = n_elements(timarr)	
yy1 = data
dxx1 = deriv_arr(xx1)
if (n_elements(step) eq 0) then step = min(dxx1)	;sample step size
print, 'Sample step size used is: ', step
;
n2 = max(xx1)/step
xx2 = dindgen(n2)*step
if (not keyword_set(random)) then begin
    yy2 = interpol(yy1, xx1, xx2)
    filltype = '(interpolated)'
end else begin
    ss2 = interpol( lindgen(n2), xx2, xx1 )         ;get subscripts of input in output array
    yy2 = fltarr(n2)
    yy2(ss2) = 1
    ss3 = where(yy2 eq 0, nss3)     ;where are the holes

    yy2(ss2) = yy1
    if (nss3 ne 0) then begin
	ss4 = randomu(xx, nss3) * (n1-1)
	yy2(ss3) = yy1(ss4)         ;put random values into
    end

    filltype = '(random fill)'
end
;
ff2 = fft(yy2,1)
freq = (dindgen(n2/2)+1) / (n2/2) * (0.5 * 1/step)
peri = 1/freq
;
out = abs(ff2(1:n2/2))^2
ytit2 = 'ABS(FFT)^2'
if (keyword_set(smoo)) then begin
    out = smooth(out, smoo)
    ytit2 = ytit2 + '  Smooth ' + strtrim(smoo,2)
end
;
xplot = freq
xtit = 'Frequency (Hz)'
if (keyword_set(mark_period)) then begin
    mark = 1./mark_period
    mark_lab = string(1/mark,format='(f8.2)') + ' sec'
end
if (keyword_set(mark_freq)) then begin
    mark = mark_freq 
    mark_lab = string(mark,format='(f8.4)') + ' Hz'
end
if (keyword_set(period)) then begin
    xplot = peri
    xtit = 'Period (sec)'
    mark = 1./mark
end
;
idev = stdev(data, iavg)
if (n_elements(ytit) eq 0) then 	ytit = ''
;
if (n_elements(tit) eq 0) then 		tit = ''
if (n_elements(xrange) eq 0) then	xrange = [min(xplot), max(xplot)]
if (n_elements(yrange) eq 0) then	yrange = [min(out), max(out)]
if (n_elements(xtype) eq 0) then	xtype = 1
if (n_elements(ytype) eq 0) then	ytype = 1
;
!p.multi = [0,1,2]
if (data_type(timarr) eq 8) then utplot, timarr, data, ytit=ytit, /ynozero $
			else plot, timarr, data, ytit=ytit, /ynozero
plottimes, .7, .95, 'Average = ' + string(iavg)
plottimes, .7, .90, 'Std Dev = ' + string(idev)
;
plot, xplot, out, ytit=ytit2, xtit=xtit, psym=10, $
		tit=tit, xrange=xrange, yrange=yrange, xtype=xtype, ytype=ytype
plottimes, .7, .95, strtrim(long(n2),2) + ' data points'
plottimes, .7, .90, string( (n2-n1)/n2*100, format='(f4.1)') + '% of data filled ' + filltype
;
for i=0,n_elements(mark)-1 do begin
    yyy = !y.crange
    if (keyword_set(ytype)) then yyy = 10.^yyy
    oplot, [1,1]*mark(i), yyy
    xyouts, mark(i), yyy(0), orient=90, mark_lab(i)
end
;
plottime, 0, 0, progver
plottime
if (keyword_set(hc)) then pprint, /reset
end