prep_data_for_df.pro (Kronos IDL library)

;+
; prep_data_for_df.pro <file desc here>
;
; :Author:
;    Philippe Zarka
;
; :History:
;    2004/07/19: Created
;
;    2004/07/19: Last Edit
;-
;
; ------------------------------
; program PREP_DATA_FOR_DF.pro
; ------------------------------

  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_DF, aj, hhd, hhf, bande, vvbuf,xtsbuf,xfbuf,msecbuf
    if i eq td then begin
        vv=vvbuf & xt=i*1.d0+xtsbuf/86400.d0 & xf=xfbuf & msec=msecbuf
    endif else begin
        vv=[vv,vvbuf] & xt=[xt,i*1.d0+xtsbuf/86400.d0]
        xf=[xf,xfbuf] & msec=[msec,msecbuf]
    endelse
  endfor
  print
  help, vv,xt,xf,msec
 
; frequency selection
  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,4)

  nom='v'+strtrim(long(aaaajjjd),2)+'_'+strtrim(fix(hd),2)+'_'+ $
        strtrim(long(aaaajjjf),2)+'_'+strtrim(fix(hf),2)+'_'+band(bande)
  path='../idlsave/'
;  path='/user4/zarka/kronos/idlsave/'

; soustraction du fond et estimation du S/B
; FREQFON(n), FON_DB(n,4), SIG_DB(n,4) [4 = +X,Z,-X,Z]

  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
      test=where(abs(f-ftest) le df/2.)
      if test(0) ne -1 then begin
	for j=0,3 do begin
	  vtest=vv(test,j)
	  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)=(vv(test,j)-fon) > 0.
	  sig=fon*(10.^(sig_db(i,j)/10.)-1.)
	  sb(test,j)=10.*alog10(vv(test,j)/sig)
	endfor
      endif
    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
	; ../idlsave/bfile.idlsave
	; contains: FREQFON(n), FON_DB(n,4), SIG_DB(n,4) [4 = +X,Z,-X,Z]
    restore,path+bfile+'.idlsave',/verb
    for i=0,3 do begin
      fon=10.^(interpol(fon_db(*,i),freqfon,xf)/10.)
      vv(*,i)=(vv(*,i)-fon) > 0.
      sig=fon*(10.^(interpol(sig_db(*,i),freqfon,xf)/10.)-1.)
      sb(*,i)=10.*alog10(vv(*,i)/sig)
    endfor
  end
  endcase

; test VV > 0
;  test=where(vv(*,0) gt 0 and vv(*,1) gt 0 and vv(*,2) gt 0 and vv(*,3) gt 0)
;  if test(0) eq -1 then stop,'test (Avv - background) > 0   not passed'
;  vv=vv(test,*) & xt=xt(test) & xf=xf(test) & msec=msec(test) & sb=sb(test,*)
;  print & print,'Avv - background > 0' & help,vv
  test=where(vv(*,0:3) lt 0)
  if test(0) ne -1 then begin
      (vv(*,0:3))(test)=0.
      sb(test)=-1.
  endif

; test S/B >= N dB
;  ant=['DF','Cal. X','Cal. Z']
;  if ib ne 0 then begin
;    print & read,'DF (0) or Antenna Calibration (1=X, 2=Z) ? ',iant
;seuil:
;    read,'S/N threshold in dB ? ',ndb
;    case iant of
; 0: test=where((sb(*,0) ge ndb or sb(*,1) ge ndb) and $
;					(sb(*,2) ge ndb and sb(*,3) ge ndb))
; 1: test=where(sb(*,1) ge ndb and sb(*,3) ge ndb)
; 2: test=where(sb(*,0) ge ndb and sb(*,2) 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,*) & xt=xt(test) & xf=xf(test) & msec=msec(test) & sb=sb(test,*)
;    print & print,ant(iant),' with S/N >= ',ndb & help,vv
;  endif

; test |cross-correlations| <= 1
;  test=where(abs(vv(*,4)) le 1 and abs(vv(*,5)) le 1 and abs(vv(*,6)) le 1 $
;	and abs(vv(*,7)) le 1)
;  if test(0) eq -1 then stop,'test |cross-correlations| <= 1   not passed'
;  vv=vv(test,*) & xt=xt(test) & xf=xf(test) & msec=msec(test) & sb=sb(test,*)
;  print & print,'|cross-correlations| <= 1' & help,vv

; test delta(ZZ) <= N%
;delta_zz:
;  print
;  read,'Maximum allowed relative difference between consecutive Azz (0.1) ? ',dz
;  test=where(abs((vv(*,1)-vv(*,3))/(vv(*,1)+vv(*,3))) le dz)
;  help,test
;  read,'Ok (1=yes, 0=no) ? ',Ok
;  if ok ne 1 then goto,delta_zz
;  if test(0) eq -1 then stop,'test delta(ZZ)   not passed'
;  vv=vv(test,*) & xt=xt(test) & xf=xf(test) & msec=msec(test) & sb=sb(test,*)
;  print & print,'delta(Azz) <= ',fix(dz*100),' %' & help,vv

; sauvegarde des donnees
  ti=lonarr(n_elements(xt))+100L & fi=ti
  data=replicate({data_n25a},n_elements(xt))
  data.ti=ti
  data.fi=fi
  data.vv=transpose(vv)
  data.sn=transpose(sb)
;  save,vv,xt,xf,msec,sb,freqfon,fon_db,sig_db, $
;	file=path+nom+'.idlsave',/verb,/xdr
  write_data_binary,'$DATA_RPWS/2001_001_090/temp/N25a_2001001.02',data

end