plot_spdyn.pro (Kronos IDL library)

;+
; Contains the plot_spdyn procedure
; and associated procedures
;
; :Author:
;    Philippe Zarka
;
; :History:
;    2004/04/07: Created
;
;    2004/04/07: Last Edit
;-

;+
; remplissage d'un spectre "peigne"
; si hfabc=1, on extrapole le remplissage jusqu'au canal de plus HF de la bande
; x: data entre 0 et xmax, -1000 dans les valeurs a inter/extra-poler
;-
  pro FILL_SPEC, x, hfabc

  nx=n_elements(x)
  test=where(x ge 0)
  if test(0) ne -1 then begin
    if hfabc then xtest=x(test(0):*) else $
				xtest=x(test(0):test(n_elements(test)-1))
    para=xtest ge 0
    ADJUST_LIN2, xtest,para, xnet
    if hfabc then x(test(0):*)=xnet else $
				x(test(0):test(n_elements(test)-1))=xnet
  endif
return
end

;+
; remplissage d'un tableau de spectres "peigne" par appels de FILL_SPEC
; x: data entre 0 et xmax, -1000 dans les valeurs a inter/extra-poler
;-
  pro FILL_TAB, x, hfabc

  if n_elements(hfabc) eq 0 then hfabc=0
  nx=n_elements(x(*,0))
  for i=0,nx-1 do begin
	xx=reform(x(i,*))
	test=where(xx lt 0)
	if test(0) ne -1 then begin
		FILL_SPEC,xx, hfabc
		x(i,*)=xx
	endif
  endfor
return
end

;--------------------------------------------
  pro APPEND_TAB, ex,ez,x1,z1,f1,x2,z2,f2
;--------------------------------------------
  f=[f1,f2] & nf=n_elements(f)
  sf=sort(f) & f=f(sf)
  test=[0b,f(1:nf-1) eq f(0:nf-2)]
  f=f(where(test eq 0)) & nf=n_elements(f)
  nf1=n_elements(f1) & if1=intarr(nf1)
  for i=0,nf1-1 do if1(i)=where(f eq f1(i))
  nf2=n_elements(f2) & if2=intarr(nf2)
  for i=0,nf2-1 do if2(i)=where(f eq f2(i))
  if ex then begin
    nb1=n_elements(x1(*,0)) & nb2=n_elements(x2(*,0))
    x=fltarr(nb1+nb2,nf)-1000
    x(0:nb1-1,if1)=x1 & x(nb1:*,if2)=x2
    x1=x
  endif
  if ez then begin
    nb1=n_elements(z1(*,0)) & nb2=n_elements(z2(*,0))
    z=fltarr(nb1+nb2,nf)-1000
    z(0:nb1-1,if1)=z1 & z(nb1:*,if2)=z2
    z1=z
  endif
  f1=f
return
end

; appelle PREP_SPDYN & SPDYNPS
; PREP_SPDYN, aaaajjj,hd,hf,bande,fmin,fmax,ex,ez,dsec, tabx,tabz,xts,xf
; SPDYNPS, image, xmin,xmax,ymin,ymax, legende_x,legende_y,titre, $
;           grid,fond,seuil,fracmin,fracmax,color,bar_titre, [/log]

;+
; trace de spectres dynamiques [ABC, H1/H2] * [Ex, Ez]
;   pour la periode, les bandes, et la resolution temporelle choisies,
;   avec possibilite d'adjonction d'un texte et de soustraction du background
;
; :Uses:
;    adjust_lin2, fond, prep_spdyn, spdynps
;
; :Params:
;    aaaajjj: in, required, type=sometype
;       date
;    hd:      in, required, type=float
;       heure de debut, peut etre fractionnaire
;    hf:      in, required, type=float
;       heure de debut, peut etre fractionnaire
;    abc:     in, required, type=sometype
;       A parameter named abc
;    h1:      in, required, type=sometype
;       A parameter named h1
;    h2:      in, required, type=sometype
;       A parameter named h2
;    fmin:    in, required, type=sometype
;       debut intervalle de frequence en kHz ou MHz (<=16.125)
;    fmax:    in, required, type=sometype
;       fin intervalle de frequence en kHz ou MHz (<=16.125)
;    ex:      in, required, type=sometype
;       A parameter named ex
;    ez:      in, required, type=sometype
;       A parameter named ez
;    dsec:    in, required, type=sometype
;       doit etre >5 sec
;    b:       in, required, type=sometype
;       b=0 (sans), b=1 (avec), b=2 (sans & avec) (texte,background)
;    texte:   in, required, type=sometype
;       A parameter named texte
;    x:       in, required, type=sometype
;       A parameter named x
;
; :Keywords:
;    nblocks: in, optional, type=sometype
;       creation du spectre dynamique en nblocks contigus
;    log:     in, optional, type=sometype
;       echelle log pour les frequences
;-
pro PLOT_SPDYN, aaaajjj,hd,hf, abc,h1,h2, fmin,fmax, ex,ez, dsec,b, texte, $
      x, nblocks=nblocks, log=log

if not(keyword_set(nblocks)) then nblocks = 1
if not(keyword_set(log)) then log=0

exname=['','X'] & ezname=['','Z']
textebande=''
aaaa=long(aaaajjj/1000L)
jjj=long(aaaajjj) mod 1000
ffmin=fmin & ffmax=fmax
if fmax le 16.125 then begin
  ffmin=ffmin*1000. & ffmax=ffmax*1000.
endif

hl = fix(hf)+(fix(hf) ne hf) - fix(hd)
if hl lt nblocks then nblocks = hl


; Bande ABC
if abc then begin
  print & print,'******************** ABC ********************' & print
  textebande=textebande+'ABC'
  PREP_SPDYN,aaaajjj,fix(hd),fix(hd)+fix(hl/nblocks),0,ffmin,ffmax,ex,ez,dsec, $
     x,z,t,f
  if nblocks gt 1 then for ib = 1,nblocks-1 do begin
    hdb = fix(hd) + fix(hl*ib/nblocks)
    hfb = fix(hd) + fix(hl*(ib+1)/nblocks)
    PREP_SPDYN,aaaajjj,hdb,hfb,0,ffmin,ffmax,ex,ez,dsec,xtmp,ztmp,ttmp,ftmp
;   x = [x, xtmp]
;   z = [z, ztmp]
    t = [t, ttmp]
    APPEND_TAB, ex,ez,x,z,f,xtmp,ztmp,ftmp
endfor
  ok=where(t ge hd*3600. and t le hf*3600.)
  t=t(ok)
  if ex then begin
    x=x(ok,*) & help,x
    test=where(x ge -999)
    if test(0) ne -1 then x(test)=x(test)-min(x(test))
    test=where(x lt 0)
    if test(0) ne -1 then FILL_TAB,x,1
  endif
  if ez then begin
    z=z(ok,*) & help,z
    test=where(z ge -999)
    if test(0) ne -1 then z(test)=z(test)-min(z(test))
    test=where(z lt 0)
    if test(0) ne -1 then FILL_TAB,z,1
  endif
  xabc=x & zabc=z
  if log then fabc=f else fabc=10.^f
  ffmin=max([ffmin,325.])		; si ABC, alors H1&2 >= 325 kHz
endif

; Bandes HF1 / HF2
if h1 or h2 then begin
  if h1 and not h2 then begin
    print & print,'******************** HF1 ********************' & print
    textebande=textebande+'H1'
    PREP_SPDYN,aaaajjj,fix(hd),fix(hd)+fix(hl/nblocks),1,ffmin,ffmax,ex,ez, $
	dsec, x,z,t,f 
    if nblocks gt 1 then for ib = 1,nblocks-1 do begin
      hdb = fix(hd) + fix(hl*ib/nblocks)
      hfb = fix(hd) + fix(hl*(ib+1)/nblocks)
      PREP_SPDYN,aaaajjj,hdb,hfb,1,ffmin,ffmax,ex,ez,dsec,xtmp,ztmp,ttmp,ftmp
;      APPEND_TAB, ex,ez,x,z,f,xtmp,ztmp,ftmp
      x = [x, xtmp]
      z = [z, ztmp]
      t = [t, ttmp]
    endfor
  endif
  if h2 and not h1 then begin
    print & print,'******************** HF2 ********************' & print
    textebande=textebande+'H2'
    PREP_SPDYN,aaaajjj,fix(hd),fix(hd)+fix(hl/nblocks),2,ffmin,ffmax,ex,ez, $
	dsec, x,z,t,f
    if nblocks gt 1 then for ib = 1,nblocks-1 do begin
      hdb = fix(hd) + fix(hl*ib/nblocks)
      hfb = fix(hd) + fix(hl*(ib+1)/nblocks)
      PREP_SPDYN,aaaajjj,hdb,hfb,2,ffmin,ffmax,ex,ez,dsec,xtmp,ztmp,ttmp,ftmp
;      APPEND_TAB, ex,ez,x,z,f,xtmp,ztmp,ftmp
      x = [x, xtmp]
      z = [z, ztmp]
      t = [t, ttmp]
    endfor
  endif
  if h1 and h2 then begin
    print & print,'******************** HF1&2 ********************' & print
    textebande=textebande+'H1H2'
    PREP_SPDYN,aaaajjj,fix(hd),fix(hd)+fix(hl/nblocks),3,ffmin,ffmax,ex,ez, $
	dsec, x,z,t,f 
    if nblocks gt 1 then for ib = 1,nblocks-1 do begin
      hdb = fix(hd) + fix(hl*ib/nblocks)
      hfb = fix(hd) + fix(hl*(ib+1)/nblocks)
      PREP_SPDYN,aaaajjj,hdb,hfb,3,ffmin,ffmax,ex,ez,dsec,xtmp,ztmp,ttmp,ftmp
;     x = [x, xtmp]
;     z = [z, ztmp]
      t = [t, ttmp]
      APPEND_TAB, ex,ez,x,z,f,xtmp,ztmp,ftmp
    endfor
  endif
  ok=where(t ge hd*3600. and t le hf*3600.)
  t=t(ok)
  if ex then begin
    x=x(ok,*) & help,x
    test=where(x ge -999)
    if test(0) ne -1 then x(test)=x(test)-min(x(test))
    test=where(x lt 0)
    if test(0) ne -1 then FILL_TAB,x
  endif
  if ez then begin
    z=z(ok,*) & help,z
    test=where(z ge -999)
    if test(0) ne -1 then z(test)=z(test)-min(z(test))
    test=where(z lt 0)
    if test(0) ne -1 then FILL_TAB,z
  endif
  xh=x & zh=z
  if log then fh=alog10(f*1.) else fh=f*1.
endif

help, fabc,fh, xabc,xh,zh
ff=[0.]
if abc then ff=[ff,fabc]
if h1 or h2 then ff=[ff,fh]
ff=ff[1:*] & nff=n_elements(ff)
nf=fix(600./max([!p.multi(2),1]))
if log then f=alog10(fmin*1.)+findgen(nf)*alog10(fmax*1./fmin)/(nf-1.) else $
  f=fmin+findgen(nf)*(fmax-fmin)/(nf-1.)
textefreq='kHz' & if fmax le 16.125 then textefreq='MHz'
if fmax le 16.125 then if log then f=f+3. else f=f*1000.
nt=n_elements(t)
aj=strtrim(aaaajjj,2)

if ex then begin
  x=fltarr(nt,nf)
  for i=0,nt-1 do begin
    xx=[0.]
    if abc then xx=[xx,reform(xabc(i,*))]
    if h1 or h2 then xx=[xx,reform(xh(i,*))]
    xx=xx[1:*]
    if nff gt 1 then x(i,*)=interpol(xx,ff,f)
  endfor
  if nt lt 1000./max([!p.multi(1),1]) then $
	x=rebin(x,nt*ceil(1000./max([!p.multi(1),1])/nt),nf)
  test=where(x lt 0)
  if test(0) ne -1 then x(test)=-1000
  xnorm=fltarr(nf)+1.
  for i=0,nf-1 do begin
    test=where(x(*,i) gt 0)
    if test(0) ne -1 then xnorm(i)=total(x(test,i))/n_elements(test)
  endfor
  xnorm=rebin(reform(xnorm,1,nf),n_elements(x(*,0)),nf)
  xn=x/xnorm
  test=where(xn lt 0)
  if test(0) ne -1 then xn(test)=-1000
  if (b mod 2) eq 0 then SPDYNPS,xn,hd,hf,fmin,fmax, $
	'Hour of '+strmid(aj,0,4)+' / '+strmid(aj,4,3), $
	'Frequency ('+textefreq+')','RPWS ['+textebande+'/Ex] '+texte, $
	0,0,0,0.01,0.93,0,'.',log=log
  if b ge 1 then SPDYNPS,x,hd,hf,fmin,fmax, $
	'Hour of '+strmid(aj,0,4)+' / '+strmid(aj,4,3), $
	'Frequency ('+textefreq+')','RPWS ['+textebande+'/Ex] (-B) '+texte, $
	0,1,-0.5,0.,0.93,0,' ',log=log
;	0,1,-0.1,0.01,0.95,0,' '
endif

if ez then begin
  z=fltarr(nt,nf)
  for i=0,nt-1 do begin
    zz=[0.]
    if abc then zz=[zz,reform(zabc(i,*))]
    if h1 or h2 then zz=[zz,reform(zh(i,*))]
    zz=zz[1:*]
    if nff gt 1 then z(i,*)=interpol(zz,ff,f)
  endfor
  if nt lt 1000./max([!p.multi(1),1]) then $
	z=rebin(z,nt*ceil(1000./max([!p.multi(1),1])/nt),nf)
  test=where(z lt 0)
  if test(0) ne -1 then z(test)=-1000
  znorm=fltarr(nf)+1.
  for i=0,nf-1 do begin
    test=where(z(*,i) gt 0)
    if test(0) ne -1 then znorm(i)=total(z(test,i))/n_elements(test)
  endfor
  znorm=rebin(reform(znorm,1,nf),n_elements(z(*,0)),nf)
  zn=z/znorm
  test=where(zn lt 0)
  if test(0) ne -1 then zn(test)=-1000
  if (b mod 2) eq 0 then SPDYNPS,zn,hd,hf,fmin,fmax, $
	'Hour of '+strmid(aj,0,4)+' / '+strmid(aj,4,3), $
	'Frequency ('+textefreq+')','RPWS ['+textebande+'/Ez] '+texte, $
	0,0,0,0.01,0.93,0,'.',log=log
  if b ge 1 then SPDYNPS,z,hd,hf,fmin,fmax, $
	'Hour of '+strmid(aj,0,4)+' / '+strmid(aj,4,3), $
	'Frequency ('+textefreq+')','RPWS ['+textebande+'/Ez] (-B) '+texte, $
	0,1,-0.5,0.,0.93,0,' ',log=log
;	0,1,-0.1,0.01,0.95,0,' '
endif

RETURN
END