band=['A','B','C','HF1','HF2']
print
read,'yyyyddd_beg, hr_beg, yyyyddd_end, hr_end ? ',aaaajjjd,hd,aaaajjjf,hf
read,'HFR band (0=A, 1=B, 2=C, 3=HF1, 4=HF2) ? ',bande
print,'min,max frequency (kHz) within ',band(bande), $
' band ? [0,0 = whole band] '
read,ffmin,ffmax
td=long(aj_t97(aaaajjjd)) & tf=long(aj_t97(aaaajjjf))
for i=td,tf do begin
if i eq td then hhd=hd else hhd=0
if i eq tf then hhf=hf else hhf=24
aj=long(t97_aj(i))
PREP_DF2, aj, hhd, hhf, bande, vvbuf,xantbuf,xtsbuf,xfbuf,msecbuf
test=where(xantbuf eq 3)
if test(0) ne -1 then vvbuf(test,0)=vvbuf(test,0)*2.
if i eq td then begin
vv=vvbuf & xant=xantbuf
xt=i*1.d0+xtsbuf/86400.d0 & xf=xfbuf & msec=msecbuf
endif else begin
vv=[vv,vvbuf] & xant=[xant,xantbuf]
xt=[xt,i*1.d0+xtsbuf/86400.d0] & xf=[xf,xfbuf] & msec=[msec,msecbuf]
endelse
endfor
print
help, vv,xant,xt,xf,msec
if ffmin ne 0. and ffmax ne 0. then begin
test=where(xf ge ffmin and xf le ffmax)
if test(0) eq -1 then stop,'test frequency range not passed'
vv=vv(test,*) & xt=xt(test) & xf=xf(test) & msec=msec(test)
print,'Frequency range selection' & help,vv & print
endif
n=n_elements(xf) & sb=fltarr(n,2)
nom='v'+strtrim(long(aaaajjjd),2)+'_'+strtrim(fix(hd),2)+'_'+ $
strtrim(long(aaaajjjf),2)+'_'+strtrim(fix(hf),2)+'_'+band(bande)
path='../idlsave/'
print & read,'substract background (0) None, (1) Auto, (2) File ? ',ib
print
case ib of
0:begin
print,'No S/N selection'
freqfon=[-1] & fon_db=[0.] & sig_db=[0.]
end
1:begin
read,'(0) gaussian or (N) lower N% of intensity histogram ? ',ib1
fracmin=ib1/100.
f=xf & df=3.
if bande le 2 then begin
f=alog10(xf) & df=0.02
endif
fmin=min(f) & fmax=max(f) & pas=long((fmax-fmin)/df)
freqfon=fmin+findgen(pas+2)*df
fon_db=fltarr(pas+2,4) & sig_db=fon_db
for i=0,pas+1 do begin
ftest=fmin+i*df
for j=0,3 do begin
if j eq 0 then test=where(abs(f-ftest) le df/2)
if j ge 1 then test=where((abs(f-ftest) le df/2) and (xant eq j))
if test(0) ne -1 then begin
vtest=vv(test,j eq 0)
if ib1 eq 0 then begin
vtest=10.*alog10(vtest)
vtestmin=min(vtest) & vtest=vtest-vtestmin+1.
FOND,vtest,vfon,vsig
vmin=vfon-1.+vtestmin
fon_db(i,j)=vmin & sig_db(i,j)=vsig
fon=10.^(vmin/10.)
endif else begin
DYNMIN,vtest,fracmin,fon
fon_db(i,j)=10.*alog10(fon) & sig_db(i,j)=1.
endelse
vv(test,j eq 0)=(vv(test,j eq 0)-fon) > 0.
sig=fon*(10.^(sig_db(i,j)/10.)-1.)
sb(test,j eq 0)=10.*alog10(vv(test,j eq 0)/sig)
endif
endfor
endfor
test=where(rebin(fon_db,pas+2,1) ne 0.)
freqfon=freqfon(test) & fon_db=fon_db(test,*) & sig_db=sig_db(test,*)
if bande le 2 then freqfon=10.^freqfon
end
2:begin
bfile=''
read,'Background file (without path / extension) ? ',bfile
restore,path+bfile+'.idlsave',/verb
for i=0,3 do begin
if i eq 0 then test=lindgen(n_elements(xf))
if i ge 1 then test=where(xant eq i)
if test(0) ne -1 then begin
fon=10.^(interpol(fon_db(*,i),freqfon,xf(test))/10.)
vv(test,i eq 0)=(vv(test,i eq 0)-fon) > 0.
sig=fon*(10.^(interpol(sig_db(*,i),freqfon,xf(test))/10.)-1.)
sb(test,i eq 0)=10.*alog10(vv(test,i eq 0)/sig)
endif
endfor
end
endcase
test=where(vv(*,0) gt 0 and vv(*,1) gt 0)
if test(0) eq -1 then stop,'test (Avv - background) > 0 not passed'
vv=vv(test,*) & xant=xant(test)
xt=xt(test) & xf=xf(test) & msec=msec(test) & sb=sb(test,*)
print & print,'Avv - background > 0' & help,vv
if ib ne 0 then begin
print & read,'S/N test on (1) or (2) antennas ? ',iant
seuil:
read,'S/N threshold in dB ? ',ndb
case iant of
1: test=where(sb(*,0) ge ndb or sb(*,1) ge ndb)
2: test=where(sb(*,0) ge ndb and sb(*,1) ge ndb)
endcase
help,test
read,'Ok (1=yes, 0=no) ? ',Ok
if ok ne 1 then goto,seuil
if test(0) eq -1 then stop,'test S/N not passed'
vv=vv(test,*) & xant=xant(test)
xt=xt(test) & xf=xf(test) & msec=msec(test) & sb=sb(test,*)
print & print,'S/N >= ',ndb,' on ',iant,' antenna(s)' & help,vv
endif
test=where(abs(vv(*,2)) le 1 and abs(vv(*,3)) le 1)
if test(0) eq -1 then stop,'test |cross-correlations| <= 1 not passed'
vv=vv(test,*) & xant=xant(test)
xt=xt(test) & xf=xf(test) & msec=msec(test) & sb=sb(test,*)
print & print,'|cross-correlations| <= 1' & help,vv
save,vv,xant,xt,xf,msec,sb,freqfon,fon_db,sig_db, $
file=path+nom+'.idlsave2',/verb,/xdr
end