;
; get 24sep04 data
;
; .compile ../labpix 
;
yymmdd1=20040924
yymmdd2=yymmdd1
nscans=arch_gettbl(yymmdd1,yymmdd1,slar,slfilear)
;
nbeams=7
nrecs=900
ind=where(slar.numrecs eq nrecs,nscans)
;
; hold the total power data for on the 1350 radar and off of it.
;
tpa1=fltarr(nrecs,nscans,nbeams) 
tpa2=fltarr(nrecs,nscans,nbeams) 
tpb1=fltarr(nrecs,nscans,nbeams) 
tpb2=fltarr(nrecs,nscans,nbeams) 
i1_1=608L
i1_2=620L
i2_1=2000L
i2_2=2800l
for iscan=0,nscans-1 do begin &$
	print,'start:',slar[ind[iscan]].scan &$
	istat=arch_getdata(slar,slfilear,ind[iscan],b,type=2) &$
	if iscan eq 0 then b0=b &$
	if (n_elements(b) ne nrecs) then continue &$
	for ibm=0,6 do begin &$
		tpa1[*,iscan,ibm]=total(b.(ibm).d[i1_1:i1_2,0],1)/(i1_2-i1_1+1.) &$
		tpb1[*,iscan,ibm]=total(b.(ibm).d[i1_1:i1_2,1],1)/(i1_2-i1_1+1.) &$
		tpa2[*,iscan,ibm]=total(b.(ibm).d[i2_1:i2_2,0],1)/(i2_2-i2_1+1.) &$
		tpb2[*,iscan,ibm]=total(b.(ibm).d[i2_1:i2_2,1],1)/(i2_2-i2_1+1.) &$
	endfor &$
endfor
;
;  plot the avgerage bandpass for pola, polBdata away from the radar
;
hard=1
bavg=coravg(b0)
frq=corfrq(bavg.b1.h)
;
; average bandpass first scan
;
hard=1
if hard then pscol,'24sep04_avgbp.ps',/full
!p.multi=[0,1,2]
hor,1335,1435
ver,0,5
plot,[0,0],[0,0],/nodata,$
	xtitle='freq [Mhz]',ytitle='spectral density',$
	title ='24sep04 a1946 avg spectrum scan 426802630L polA'
for i=0,6 do oplot,frq,bavg.(i).d[*,0],col=i+1
;flag,[frq[i1_1],frq[i1_2]],linestyle=2
flag,[frq[i2_1],frq[i2_2]],linestyle=2,color=2
;
plot,[0,0],[0,0],/nodata,$
	xtitle='freq [Mhz]',ytitle='spectral density',$
	title ='24sep04 a1946 avg spectrum scan 426802630L polA'
for i=0,6 do oplot,frq,bavg.(i).d[*,1],col=i+1
;flag,[frq[i1_1],frq[i1_2]],linestyle=2
flag,[frq[i2_1],frq[i2_2]],linestyle=2,color=2
labpix
if hard then hardcopy
;
; normalize to median power no cals
;
ya2=tpa2
yb2=tpb2
ya1=tpa1
yb1=tpb1
for ibm=0,6 do begin &$
	for iscan=0,nscans-1 do begin &$
		y=ya1[*,iscan,ibm] &$
		y=y/median(y) - 1. &$
		ind=where(y lt .1) &$
		avg=meanrob(ya1[ind,iscan,ibm]) &$
		ya1[*,iscan,ibm]= ya1[*,iscan,ibm]/avg - 1. &$
		y=yb1[*,iscan,ibm] &$
		y=y/median(y) - 1. &$
		ind=where(y lt .1) &$
		avg=meanrob(yb1[ind,iscan,ibm]) &$
		yb1[*,iscan,ibm]= yb1[*,iscan,ibm]/avg - 1. &$
;
		y=ya2[*,iscan,ibm] &$
		y=y/median(y) - 1. &$
		ind=where(y lt .1) &$
		avg=meanrob(ya2[ind,iscan,ibm]) &$
		ya2[*,iscan,ibm]= ya2[*,iscan,ibm]/avg - 1. &$
		y=yb2[*,iscan,ibm] &$
		y=y/median(y) - 1. &$
		ind=where(y lt .1) &$
		avg=meanrob(yb2[ind,iscan,ibm]) &$
		yb2[*,iscan,ibm]= yb2[*,iscan,ibm]/avg - 1. &$
	endfor &$
endfor
v1=.2
v2=1.
var=[[1,2],[1,1],[1,1],[1,1],[1,1],[1,1],[1,1]]
hard=1
if hard then pscol,'24sep04_faaradar.ps',/full
x=findgen(nrecs)
hor
for ibm=0,6 do begin &$
tit1='pwr 1395 Mhz band for 9 scans ' &$
!p.multi=[0,1,2] &$
print,ibm &$
v=(var[0,ibm] eq 1)? v1:v2 &$
ver,0,v &$
inc=v*.1 &$
lab=string(format='("Pixel ",i1)',ibm) &$
stripsxy,x,ya2[*,*,ibm],inc/2,inc,/step,$
	xtitle='time [seconds for 1 scan]',$
	ytitle='pwr [Tsys]',$
	title=tit1+lab+' polA' &$
v=(var[1,ibm] eq 1)? v1:v2 &$
ver,0,v &$
inc=v*.1 &$
stripsxy,x,yb2[*,*,ibm],inc/2,inc,/step,$
	xtitle='time [seconds for 1 scan]',$
	ytitle='pwr [Tsys]',$
	title=tit1+lab+' polB' &$
endfor
if hard then hardcopy
x
;---------------------------------------------------------------------
; using a blowup show:
; alignment of radar and out of band
; blowup of the recovery time 
; flat spectra.
;
if hard then pscol,'24sep04_pwrjumps.ps',/full
hor,0,60
cs=1.5
!p.multi=[0,1,3]
ver,-.05,.8
sym=10
plot ,yb2[*,0,0],psym=sym,charsize=cs,$
	xtitle='Sample [secs]',ytitle='power [Tsys]',$
	title='Jumps align in time with 1350 radar (pwr vs time)'
oplot,yb1[*,0,0],color=2,psym=sym
ln=2
xp=.04
scl=.6
note,ln      ,'150Khz about 1350 Radar',xp=xp,color=2
note,ln+1*scl,'20Mhz away from 1350 Radar',xp=xp
flag,findgen(5)*12+8.5,linestyle=1
;
ver,-.03,.05
plot ,yb2[*,0,0],psym=sym,charsize=cs,$
	xtitle='Sample [secs]',ytitle='power [Tsys]',$
	title='pwr vs time. blowup showing recovery time for jump'
oplot,yb1[*,0,0],color=2,psym=sym
;
ver,0,.2
hor
y=b0[452].b1.d[*,1]/b0[451].b1.d[*,1]-1.
tosmo=13
plot,frq,smooth(y,tosmo),charsize=cs,$
	xtitle='freq [Mhz]',ytitle='spectral density [Tsys]',$
	title=$
'SpectrumAtJump/spectrumBeforeJump -1. Jump is flat in freq (smp 452)'
if hard then hardcopy
x
;---------------------------------------------------------------------
;
; show correlation jump and strength of radar
;
if hard then pscol,'24sep04_corAmpl.ps',/full 
hor,.4,1.5
ver,-.02,.15 &$
!p.multi=[0,1,2]
clip1=.4
clip2=.2
tit1='Strength outside radar band vs strength radar '
for ibm=0,6 do begin &$
	y1=ya1[*,*,ibm] &$
	y2=ya2[*,*,ibm] &$
	ind=where((y1 gt clip1) and (y2 lt clip2)) &$
	lab=string(format='("Pixel ",i1," polA")',ibm) &$
	plot,y1[ind],y2[ind],psym=2,$
	xtitle='pwr in radar [Tsys]',$
	ytitle='pwr outside  radar [Tsys]',$
	title=tit1+ lab &$
 	y1=yb1[*,*,ibm] &$
 	y2=yb2[*,*,ibm] &$
	ind=where((y1 gt clip1) and (y2 lt clip2)) &$
	lab=string(format='("Pixel ",i1," polB")',ibm) &$
	plot,y1[ind],y2[ind],psym=2, $
	xtitle='pwr in radar [Tsys]',$
	ytitle='pwr outside  radar [Tsys]',title=tit1 + lab &$
endfor
if hard then hardcopy
x
 
end