auto_prep_data_for_df.pro (Kronos IDL library)

;+
; auto_prep_data_for_df.pro <file desc here>
;
; :Author:
;    Baptiste Cecconi, Philippe Zarka
;
; :History:
;    2003/12/08: Created
;
;    2003/12/08: Last Edit
;-
;
;+
; auto_prep_data_for_df is a procedure that <behavior desc here>
;
; :Uses:
;    dynmin, fond, prep_df
;
; :Params:
;    aaaajjj: in, required, type=sometype
;       A parameter named aaaajjj
;    bande:   in, required, type=sometype
;       A parameter named bande
;    iant:    in, required, type=sometype
;       A parameter named iant
;    ndb:     in, required, type=sometype
;       A parameter named ndb
;-
pro AUTO_PREP_DATA_FOR_DF, aaaajjj,bande,iant,ndb
; automated version of PREP_DATA_FOR_DF
; day per day (aaaajjj)
; 1 band (0=A, 1=B, 2=C, 3=HF1, 4=HF2), no frequency selection
; substracting Jupiter background
; with S/N >= ndb (dB) for 
; DF measurements (0) or Antenna Calibration (1=X, 2=Z)
; with Dzz <= 10%

  band=['A','B','C','HF1','HF2']

;  print
;  read,'yyyyddd_beg, hr_beg, yyyyddd_end, hr_end ? ',aaaajjjd,hd,aaaajjjf,hf
  aaaajjjd=aaaajjj & aaaajjjf=aaaajjj & hd=0 & hf=24
;  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
  ffmin=0 & ffmax=0
  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/'

; 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
  ib=2
;  print

  case ib of
0:begin
    print,'No S/N selection'
    freqfon=[-1] & fon_db=[0.] & sig_db=[0.]
  end
1:begin
    ib1 = 10
    ;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]
    bfile='fond_jupiter_'+band(bande)
    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 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
;    iant=0
seuil:
;    read,'S/N threshold in dB ? ',ndb
;  ndb = 10.
    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
  dz=0.1
  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
  save,vv,xt,xf,msec,sb,freqfon,fon_db,sig_db, $
	file=path+nom+'.idlsave',/verb,/xdr
end