;
; Computes and plots the frequency response of various filters.
;
; We want lots of overlapping band-passed filters, whose width increases
; linearly with time scale (from a width of 10 for the 0-10 year timescale,
; to a width of 50 for the 50-100 year timescale, noting that only an 80-yr
; run of data is analysed).
;
; Define filters to apply
;
allmin=[findgen(26),findgen(8)*2+26,replicate(40,12)]
allmax=[findgen(26)*2+10,findgen(8)*5+65,110.,replicate(120,11)]
allts=[replicate(0.,26+8+1),findgen(11)*0.1]
nband=n_elements(allmin)
print,allmin
print
print,allmax
print
print,allts
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=8,ncol=3,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=850,xsize=700
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=11
endelse
def_1color,10,color='blue'
def_1color,11,color='red'
def_1color,12,color='mlgrey'
;
perlist=findgen(148)+3.
nper=n_elements(perlist)
nt=1000
;
; Consider each filter separately
;
for iband = 0 , nband-1 do begin
  t1=allmin(iband)
  t2=allmax(iband)
  ts=allts(iband)
  ;
  ; Generate sine waves at different frequencies and filter these
  ;
  fresp=fltarr(nper)
  for iper = 0 , nper-1 do begin
    ;
    y=sin(findgen(nt)*2*!pi/perlist(iper))
    ;
    if t1 gt 0 then begin
      filter_cru,t1,tsin=y,tshigh=yh1,tslow=yl1
    endif else begin
      yh1=y*0.
      yl1=y
    endelse
    ;
    if t2 ne 999 then begin
      filter_cru,t2,tsin=y,tshigh=yh2,tslow=yl2
    endif else begin
      yh2=y
      yl2=y*0.
    endelse
    ;
    yfilt=yl1-yl2+ts*yl2
    ;
    dummy=moment(y,sdev=sd1)
    dummy=moment(yfilt,sdev=sd2)
    fresp(iper)=sd2/sd1
    ;
  endfor
  ;
  pause
  plot,perlist,fresp,$
    /xstyle,xrange=[0,150],xtitle='Period  (years)',$
    /ystyle,yrange=[-0.1,1.1],ytitle='Filter response',$
    title=string([t1,t2,ts],format='(2I5,F5.1)')
  oplot,!x.crange,[0,0],linestyle=1
  oplot,!x.crange,[1,1],linestyle=1
  oplot,[t1,t1],!y.crange,linestyle=1
  oplot,[t2,t2],!y.crange,linestyle=1
  ;
endfor
;
end