prep_spdyn.pro (Kronos IDL library)

;+
; Contains the prep_spdyn procedure
;
; :Author:
;    Philippe Zarka
;
; :History:
;    2006/12/28: Created
;
;    2006/12/28: Last Edit
;-
;
;+
; Kronos  (~NewBench) pour la periode choisie "aaaajjj hd -> hf"
; et les bandes choisies (0=ABC, 1=HF1, 2=HF2, 3=HF1&2)
; + calibre et remplit les tableaux tabx (si ex <> 0) & tabz (si ez <> 0)
; a la resolution temporelle dsec (sec) si dsec > 5, ou dsec*dt_min (sec)
; si dsec <= 5
;
; :Uses:
;    cal, freq_abc, lit_hfr_liste
;
; :Params:
;    aaaajjj: in, required, type=sometype
;       date
;    hd:      in, required, type=sometype
;       heure de debut
;    hf:      in, required, type=sometype
;       heure de fin
;    bande:   in, required, type=sometype
;       bandes choisies (0=ABC, 1=HF1, 2=HF2, 3=HF1&2)
;    fmin:    in, required, type=sometype
;       en kHz (ne sert qu'en H1 et H2)
;    fmax:    in, required, type=sometype
;       en kHz (ne sert qu'en H1 et H2)
;    ex:      in, required, type=sometype
;       calibre et remplit les tableaux tabx (si ex <> 0)
;    ez:      in, required, type=sometype
;       calibre et remplit les tableaux tabz (si ez <> 0)
;    dsec:    in, required, type=sometype
;       resolution temporelle (s) si dsec > 5, ou dsec*dt_min (sec) si dsec <= 5
;    tabx:    out, required, type=sometype
;       tableau rempli et calibré
;    tabz:    in, required, type=sometype
;       tableau rempli et calibré
;    xts:     in, required, type=sometype
;       rampes de temps (sec)
;    xf:      in, required, type=sometype
;       rampes de frequences (kHz [H] ou log(kHz) [ABC])
;-
pro PREP_SPDYN, aaaajjj,hd,hf,bande,fmin,fmax,ex,ez,dsec, tabx,tabz,xts,xf

  common KRONOSCAL, a1h1,a1h2,ai,dbcal_abc,dbcal_h,dbcal

; lecture des donnees

  if bande eq 0 then LIT_HFR_LISTE, aaaajjj, hd, hf, 0, 2, z
  if bande eq 1 or bande eq 2 then $
	LIT_HFR_LISTE, aaaajjj, hd, hf, bande+2, bande+2, z
  if bande eq 3 then LIT_HFR_LISTE, aaaajjj, hd, hf, 3, 4, z
  if bande lt 0 or bande gt 3 then stop,'erreur bande'
  nz = n_elements(z)
  if nz le 1 then begin
      nb = long((hf-hd)*3600./dsec)
      xts = hd*3600.+findgen(nb)*dsec
      xf=[fmin]
      if ex then tabx=fltarr(nb,1) else tabx=0
      if ez then tabz=fltarr(nb,1) else tabz=0
      return
  endif

; rampes de temps & frequences

  time = dblarr(nz)
  for i=0L,nz-1L do time(i)=total(transpose(z(i).time)#[3600.0,60.0,1.0,0.01])
  dt=time(1:nz-1)-time(0:nz-2) & dt_min=min(dt(where(dt gt 0.)))
  print,'Delta t min (sec) = ',dt_min
  ddsec=dsec
;  read, 'Delta t (sec ou N[<=5]xdt_min) ? ', ddsec
  if ddsec le 5 then ddsec=ddsec*dt_min
  print,'dsec = ',ddsec
  nb = long((hf-hd)*3600./ddsec)
  xts = hd*3600.+findgen(nb)*ddsec
  ok=where(time ge xts(0) and time lt xts(nb-1)+ddsec)
; help,ok
  z=z(ok) & time=time(ok) & nz = n_elements(z)

  if bande eq 0 then begin
    nfmax=max(z.filters)
    nf = 3*nfmax
    xf = FREQ_ABC(nfmax)
    xf = alog10(xf)
  endif else begin
    freq=z.channel
    ok=where(freq ge fmin and freq le fmax)
;   help,ok
    z=z(ok) & time=time(ok) & nz = n_elements(z)
    freq=z.channel
;    if bande eq 3 then begin
;      fmaxh1=max(freq(where(z.bande eq 3)))
;      ok=where(z.bande eq 3 or (z.bande eq 4 and freq gt fmaxh1))
;      z=z(ok) & time=time(ok) & nz = n_elements(z)
;    endif
;    freq=z.channel
    ffmin=min(freq) & ffmax=max(freq)
    df=freq(1:nz-1)-freq(0:nz-2) & df=min(df(where(df gt 10)))
    nfmax=max(z.filters)
    print,'fmin, fmax, df, #filters =',ffmin,ffmax,df,nfmax
    nf = ((ffmax-ffmin)/df + 1)*nfmax
    while nf gt 650 do begin
	nf=nf/2 & nfmax=nfmax/2
    endwhile
    ffmax=ffmax+25*(nfmax-1)/nfmax
    xf = ffmin + findgen(nf)*(ffmax-ffmin)/(nf-1)
  endelse
  print,'N (balayages x frequencies) = ',nb,nf
  jb=long((time-xts(0))/ddsec)
  if bande eq 0 then jf=z.bande*nfmax+z.filter*(nfmax/z.filters) else $
	jf=long((z.channel-ffmin)/df)*nfmax+z.filter*(nfmax/z.filters)

; donnees X & Z

  if ex then begin
    tabx=fltarr(nb,nf) & ntabx=intarr(nb,nf)
    xant=((z.df mod 10b)-1b)*(z.agc1 ne 255b)	; +X,-X,D = 0,1,2
    dbx=CAL(z.agc1,z.auto1,z.bande,xant,z.channel,z.filters,z.filter)
    ndbx=(bytarr(nz)+1b) and ((z.agc1 ne 255b) or (z.auto1 ne 255b))
  endif

  if ez then begin
    tabz=fltarr(nb,nf) & ntabz=intarr(nb,nf)
    zant=bytarr(nz)+3b				; Z = 3
    dbz=CAL(z.agc2,z.auto2,z.bande,zant,z.channel,z.filters,z.filter)
    ndbz=(bytarr(nz)+1b) and ((z.agc2 ne 255b) or (z.auto2 ne 255b))
  endif

  for i=0L,nz-1L do begin
;   if i mod 10000L eq 0 then print,i
    if ex then begin
      tabx(jb(i),jf(i))=tabx(jb(i),jf(i))+dbx(i)*ndbx(i)
      ntabx(jb(i),jf(i))=ntabx(jb(i),jf(i))+ndbx(i)
    endif
    if ez then begin
      tabz(jb(i),jf(i))=tabz(jb(i),jf(i))+dbz(i)*ndbz(i)
      ntabz(jb(i),jf(i))=ntabz(jb(i),jf(i))+ndbz(i)
    endif
  endfor

; mise a -1000 des pixels "vides" dans le spectre dynamique

  if ex then begin
    test=where(ntabx eq 0)
    if test(0) ne -1 then begin
	  tabx(test)=-1000 & ntabx(test)=1
    endif
    tabx=tabx/ntabx
  endif

  if ez then begin
    test=where(ntabz eq 0)
    if test(0) ne -1 then begin
	  tabz(test)=-1000 & ntabz(test)=1
    endif
    tabz=tabz/ntabz
  endif

  if ex eq 0 then tabx=0 & if ez eq 0 then tabz=0
;  help,tabx,tabz,xts,xf

return
end