;-----------  RESONANCE/ELMAVAN tests KIA -------------
; reads and processes tm files
; history
; 	2011-11-02  TM ver 14, w/o HK packets, used for LEMI interface tests
;   2011-11-11  TM ver 15, used for SUSPI interface tests
;   2012-03-15  TM ver 20
;   2012-09-01  processing of ELMAVAN Ground added
;   2012-09-20  spectrograms added
;   2012-09-29  overview spectrograms, log added
;   2012-10-01  log changed to elm.txt
;   2012-10-06 bug in insertion of data gaps fixed in elmspg.ps
;
function ver
  return,'2012-10-06'
end

;--------------------------------- kia procedures ----------------------------

pro one, name,sub=sub
  ; runs one measurement using kia.exe (must be placed in the actual dir)
  ; commands are read from comm.txt (must be placed in the actual dir)
  ; generates output files
  ;     elm'yyyymmdd_hhmmss_name'_dat.bin
  ;     elm'yyyymmdd_hhmmss_name'_log.txt
  ;     xxx__elm'yyyymmdd_hhmmss_name'_dat.ps, if waveforms are measured, xxx=001,017,
  ;
  ; ___ run kia
  if not keyword_set(name) then name=''
  spawn,'kyja com23 com24 comm.txt '+name
  openr,1,'output_file_name.txt' & fnam='' & readf,1,fnam & close,1
  ; ___ process data
  f=file_search('elm'+fnam+'*.bin', count=n)
  for i=0,n-1 do process, f(i)
end




pro processall, ground=ground
  ;plot_init, 'all.ps'
  files=file_search('*.bin',count=nf)
  for i=0,nf-1 do process, files(i),ground=ground
  ;plot_close
end



pro process,tmfile,ground=ground
;
; ELMAVAN TM file processing
; reads the telemetry stream
; writes files according to the ID bytes in separate
;  readable outputs
;
  parse,tmfile, nid,ids,nrecs,types,fss,sources,files,t0wf
  print,'---- ELMAVAN TM Processing ----'
 for i=0,nid-1 do begin
    case types(i) of
    0: begin ;______HK
         rdhk, files(i),hk,ok=ok
         if ok then writehk,hk,files(i)
       end
    1: begin; ______WF
    	 rdwf, files(i),wf,ok=ok
    	 if ok then $
    	   if keyword_set(ground) then begin
    	     tit='3D WAVEFORMS ('+decodeid(wf(0).id)+') '+$
    	         strtrim(string(nrecs(i),form='(I15)'),2)+' packets of 144 samples from '+t0wf+'.'
    	     analwf,wf,fss(i),files(i),tit,ip,np,t0wf
    	     plotg,wf,fss(i),files(i),tit,ip,np,t0wf,/logf,sub=10;,donotclose=donotclose,sub=sub,logf=logf
    	   end else begin
    	     tit='7D WAVEFORMS ('+decodeid(wf(0).id)+') '+$
    	         strtrim(string(nrecs(i),form='(I15)'),2)+' packets of 144 samples from '+t0wf+'.'
    	     analwf,wf,fss(i),files(i),tit,ip,np,t0wf
    	     plotwf,wf,fss(i),files(i),tit,ip,np,t0wf,/logf,sub=10;,donotclose=donotclose,sub=sub,logf=logf
    	   end $
         else begin
           print,' ---- Reading '+files[i]+' FAILED.'
         end
       end
    2: begin ;_______SM
       end
    3: begin ;_______SP
       end
    7: if ids(i) eq 255 then begin ;_______link test
         rdlink, files(i),link,ok=ok
         if ok then analink,link,files(i)
       end else print,' ---- UNKNOWN link test ID:',ids(i)
    else: begin
           print,' ---- UNKNOWN data type:',types(i)
       end
    endcase
  end
  print,'---- Processing done. '
  print
end



pro parse,tmfile, nid,ids,nrecs,types,fss,sources,files,t0wf
;
; ELMAVAN TM file parser
; reads the telemetry stream
; writes files according to the ID bytes in separate records
;
; nid... number of differrent ids found in the tmfile
; ids ... array of id
; types... first 3bits of id - determines the processing method: HK,WF,SM,SP
; fss  ...sampling frequency, per id: 5e4/2.5e4 (Hz)
; sources...test data bit, per id: source=0 measured signal, =1 test data
; nrecs ... number of records per id
; files ... filename per id
;
; ID
  stypes=['-HK-','-WF-','-SM-','-SP-','????','????','????','LINK']
; TM structure:
  pac={id:0b,d:bytarr(2031)}
  ;
  print,'---- ELMAVAN TM Parser ----'
  print,'---- TM file: ',tmfile
 openr,1,tmfile
  s=fstat(1)
  ndat= s.size /2032L
  print,'---- Number of TM packets: ',ndat
  if ndat gt 0 then begin
    dat=replicate(pac,ndat)
    readu,1,dat
    ;
   if ndat gt 1 then h=histogram(dat.id,bin=1b,min=0b,max=255b) $
     else begin
      h=lonarr(256) & h(dat.id)=1
     end
   ids=where(h ne 0,nid)
   nrecs=h(ids)
   files=strarr(nid)
   types=bytarr(nid)
   fss=fltarr(nid)
   sources=bytarr(nid)
   for i=0,nid-1 do begin
    if (ids(i) and 8b) eq 8b then fss(i)=2.5e4 else  fss(i)=5e4
    if (ids(i) and 128b) eq 128b then sources(i)=1b else  sources(i)=0b
    irec=where(dat.id eq ids(i),nrec2)
    if nrec2 ne nrecs(i) then print,'------- !!!!',nrecs(i),'<>',nrec2
    types(i)= ids(i) and 7b
    files(i)=string(ids(i),form='(i3.3)')+'__'+tmfile
    mes='----    ID='+string(ids(i),form='(i3.3)')+' ('+stypes(types(i))
    if (types(i) ge 1)and(types(i) le 3) then mes=mes+decodeid(ids(i)) ; decode only for WF,SM and SP
    mes=mes+'): '+string(nrec2)+' packets -> '+files(i)
    print,mes
    openw,2,files(i)
    writeu,2,dat(irec)
    close,2
   end
  end else nid=0
  close,1
  ;
   logfile=strmid(tmfile,0,strpos(/reverse_s,tmfile,'.bin')-5)+'log.txt'
   print,'---- Log file: ',logfile
   openr,1,logfile & t0s='' & a='' & addt0=0d0
   readf,1,t0s,form='(///21x,a26)'
   if strmid(t0s,0,1) eq 'e' then t0s=strmid(t0s,3,23) $ ;old version, form='(///24x,a23)'
                             else t0s=strmid(t0s,0,23)
   while (strpos(a,'data=  20') eq -1) and not eof(1) do readf,1,a
    if a ne '' then reads,a,addt0, form='(g12)'
   close,1
   str_tns, t0s+'.000.000',tns & tns=tns+ulong64(addt0*1d9)
   tns_str,tns,t0wf & t0wf=strmid(t0wf,0,27)
   ;
   print,'---- Parsing done.'
   print
  end


function decodeid, id
  tst=['Measured data','Artificial data']
  fs=['50 kHz','25 kHz']
  fr=[ ['781-Hz', '781-Hz'] , ['195-Hz','97-Hz'] ]
  art=(id and 128b) eq 128b
  sid=tst(art)
  if art then begin
     if (id and 64b) eq 64b then sid=sid+', 781-Hz squares to antenna'
     if (id and 16b) eq 16b $
         then sid=sid+', '+fr[(id and 8b) eq 8b,(id and 32b) eq 32b]+' ramp'$
         else sid=sid+', channel constants'
     if (id and 32b) eq 32b then sid=sid+', without decimation' else sid=sid+', with decimation'
  end
  sid=sid+'  sampled at '+fs((id and 8b) eq 8b)
  return, sid
end


; --------------------- HK reading & diagnostics --------------------


pro rdhk, file, hk, ok=ok
;
; reading the telemetry stream of ELMAVAN HK packets
;
;-----------------
; structure of settings inside HK packets
  set={gdig:0b, wf_sram:0b, wf_g:0b, wf_i:0b, s_n:0b, s_tab:0b, sp_ant:0b,$
     ic7:0b,ic5:0b,ic4:0b,sc:0b, reg:0u, other:bytarr(3)}
  crc={ crckHz200h:0u, crckHz200:0ul, addr:0ul, datah:0b, data:0ul, crc:0b}
; HK packets
  hkpac={id:0b,  timh:0b, tim:0ul, kHz200h:0u, kHz200:0ul,$
       trigger:0b,nosuspi:0u, calkHz200h:0u, calkHz200:0ul, cal1Hz:0ul,$
       i33:0u,i15:0u, v33:0u, v15:0u, v50:0u, tmp:0u, rad:0ul,$
       curset:set, newset:set, nuks:0b, npatch:0b, ncrcm:0b, ncrcb:0b, $
       uks:bytarr(16),patch:bytarr(66),$
       crcm:replicate(crc,87), crcb:replicate(crc,30), ee:0b }
;-----------------
; reading TMfile, extracting HK
;
  print,'----    Reading housekeeping TM file: ',file
  openr,1,file, /SWAP_IF_LITTLE_ENDIAN
   s=fstat(1)
  ndat= s.size /2032L
  if ndat gt 0 then begin
    hk=replicate(hkpac,ndat)
    readu,1,hk
    ok=1
  end else ok=0
  close,1
end


pro writehk,hk,file
; writing HK information to an ASCII outfile
; d .. HK structure
;
  outfile=strmid(file,0,strpos(file,'.',/reverse_s))+'.txt'
  print,'----    Writing housekeeping ASCII file: ',outfile
  nhk=n_elements(hk)
  openw,1,outfile
    printf,1,'-------------------------------------------'
    printf,1,'----- Resonance/ELMAVAN HK ASCII file -----'
    printf,1,'-- Generated on '          +systime()+ ' --'
    printf,1,'-------------------------------------------'
    printf,1,' ' &    printf,1,'   Source: '+file
    printf,1, nhk,' HK packets. '
  for i=0,nhk-1 do begin
    printf,1,' '&  printf,1,'------------- ',i+1,' --------------'
     printf,1,'  s/c time:', strtrim(string(hk(i).timh*'100000000'XULL+hk(i).tim),2)
     printf,1,'  200 kHz counter: ', strtrim(string(hk(i).kHz200h*'100000000'XULL+hk(i).kHz200),2)
     printf,1,'  Trigger:',hk(i).trigger
     printf,1,'  No. of SUSPI packets:',hk(i).nosuspi
     if i gt 0 then difcal=hk(i).cal1Hz-hk(i-1).cal1Hz else difcal=0
     if difcal ne 0 then rat='    200kHz/1Hz: '+$
          strtrim(string(( hk(i).calkHz200h*'100000000'XULL+hk(i).calkHz200  -  $
                hk(i-1).calkHz200h*'100000000'XULL-hk(i-1).calkHz200)/2e3/$
              difcal),2)+'%'$
      else rat=''
     printf,1,'  CAL-200kHz: '+strtrim(string(hk(i).calkHz200h*'100000000'XULL+hk(i).calkHz200),2)+$
              '    CAL-1Hz: '+strtrim(string(hk(i).cal1Hz),2)+rat
     printf,1,'  Current 3.3V (mA-TBU): ',strtrim(string(hk(i).i33*0.6495,form='(i20)'),2)
     printf,1,'  Current 1.5V (mA-TBU): ',strtrim(string(hk(i).i15*0.6495,form='(i20)'),2)
     printf,1,'  Voltage 3.3V (mV): ',strtrim(string(hk(i).v33* 1.3431,form='(i20)' ),2)
     printf,1,'  Voltage 1.5V (mV): ',strtrim(string(hk(i).v15* 0.6105,form='(i20)'),2)
     printf,1,'  Voltage 5.0V (mV): ',strtrim(string(hk(i).v50* 2.0350,form='(i20)'),2)
     printf,1,'  Temperature (deg C) : ',strtrim(string(hk(i).tmp*0.12989 - 273.15,form='(i10)'),2)
     printf,1,'  Radiation count: ',strtrim(string(hk(i).rad),2)
     printf,1,'  Previous regular STOP (FFFF): '+string(hk(i).curset.reg,form='(z4.4)')
     wrset,1,hk(i).curset,'  Current settings -> '
     wrset,1,hk(i).newset,'  Waiting settings -> '
     if hk(i).nuks le 0 then uks='' else begin
       uks='  UKS data: '
       for j=0,hk(i).nuks-1 < 7 do $
          uks=uks+strtrim(string(hk(i).uks(2*j),form='(i3)'),2)+$
                  '-'+strtrim(string(hk(i).uks(2*j+1),form='(i3)'),2)+' '
     end
     printf,1,'  Number of UKS packets: '+strtrim(string(hk(i).nuks,form='(i3)'),2)+uks
     if hk(i).npatch le 0 then patch='' else begin
       patch='  Patch headers: '
       for j=0,hk(i).npatch-1 < 65 do patch=patch+strtrim(string(hk(i).patch(j),form='(i3)'),2)+' '
     end
     printf,1,'  Number of Patch packets: '+strtrim(string(hk(i).npatch,form='(i3)'),2)+patch
     if hk(i).ncrcm gt 0 then begin
        crc='  CRC errors (200 kHz counter, address, data, crc) : '
        for j=0,hk(i).ncrcm-1 <86 do crc=crc+$
        strtrim(string(hk(i).crcm(j).crckHz200h*'100000000'XULL+hk(i).crcm(j).crckHz200),2)+$
        ', '+string(hk(i).crcm(j).addr,form='(z8.8)')+', '+$
        string(hk(i).crcm(j).datah*'1000000'XULL+hk(i).crcm(j).data,form='(z10.10)')+', '+$
        string(hk(i).crcm(j).crc,form='(z2.2)')+';   '
      end else crc=''
     printf,1,'  Number of CRC errors in MRAM: '+strtrim(string(hk(i).ncrcm,form='(i3)'),2)+crc
     if hk(i).ncrcb gt 0 then begin
        crc='  CRC errors (200 kHz counter, address, data, crc) : '
        for j=0,hk(i).ncrcb-1 <29 do crc=crc+$
        strtrim(string(hk(i).crcb(j).crckHz200h*'100000000'XULL+hk(i).crcb(j).crckHz200),2)+$
        ', '+string(hk(i).crcb(j).addr,form='(z8.8)')+', '+$
        string(hk(i).crcb(j).datah*'1000000'XULL+hk(i).crcb(j).data,form='(z10.10)')+', '+$
        string(hk(i).crcb(j).crc,form='(z2.2)')+';   '
      end else crc=''
     printf,1,'  Number of CRC errors in BRAM: '+strtrim(string(hk(i).ncrcb,form='(i3)'),2)+crc
     printf,1,'  EE (238):',hk(i).ee
  end
  close,1
end



pro wrset,u,s,tit
; write settings s to unit u
  printf,u, tit+'  Data:'+decodeid(s.gdig)
  printf,u,'                      SRAM:'+string(s.wf_sram,form='(z2.2)')+' G:'+string(s.wf_g,form='(i3)')+$
     ' I:'+string(s.wf_i,form='(i3)')+$
  ;' SM/SP:',strtrim(string(s.smsp),2)+
  ' Analog(IC7:'+string(s.ic7,form='(z2.2)')+' IC5:'+string(s.ic5,form='(z2.2)')+$
  ' IC4:'+string(s.ic4,form='(z2.2)')+' SC:'+string(s.sc,form='(z2.2)')+')
end




; --------------------- WF reading & diagnostics --------------------



pro rdwf, file,wf,ok=ok
;
; reads the telemetry stream of WF ELMAVAN data packets
;
; TM structure of WF data:
  datpac={id:0b,  timh:0b, tim:0ul, kHz200h:0u, kHz200:0ul,$
       d:intarr(7,144), adr:0u,crc:0u}
;--------------------------
  print,'----    Reading waveform TM file:     ',file
   openr,1,file, /SWAP_IF_LITTLE_ENDIAN
   s=fstat(1)
  nwf= s.size /2032L
  if nwf gt 0 then begin
    wf=replicate(datpac,nwf)
    readu,1,wf
    ok=1
  end else ok=0
  close,1
end


pro analwf,wf,fs,file,title,ip,np,t0wf
  ; WF timing analysis
  ;
  nwf=n_elements(wf)
  txtfile=strmid(file,0,strpos(file,'.',/reverse_s))+'.txt'
  print,'----    Writing waveform analysis file: ',txtfile
  openw,1,txtfile
      printf,1,'-------------------------------------------'
      printf,1,'-----  Resonance/ELMAVAN WF analysis  -----'
      printf,1,'-- Generated on '          +systime()+ ' --'
      printf,1,'-------------------------------------------'
      printf,1,' ' &    printf,1,'   Source: '+file
      printf,1, '   '+title
    ;
    if fs eq 2.5e4 then pak200=1152UL else  pak200=576UL ;(144*200/50)
;     pak200=144
    dif=wf(1:nwf-1).kHz200 - wf(0:nwf-2).kHz200
    mdif=median(dif)
    printf,1,'   Median difference of the 200-kHz counter in neigbouring packets: '+string(mdif,form='(i15)')
    printf,1,'   Expected difference of the 200-kHz counter for continuous data : '+string(pak200,form='(i15)')
    ;
    printf,1
    h=histogram((wf.adr and 'f000'xu)/'1000'xu,bin=1,min=0,max=15)
    ih=where(h ne 0,nh) & sh='' & for i=0,nh-1 do sh=sh+string(ih(i),form='(i3)')
    printf,1,'   SRAM used (0-7): '+sh
    h=histogram(long(wf(1:nwf-1).adr)- long(wf(0:nwf-2).adr),bin=1l,min=-65535l,max=65535l)
    ih=where(h ne 0,nh) & sh='' & for i=0,nh-1 do sh=sh+string(ih(i)-65535l,form='(i7)')
    printf,1,'   SRAM address increments: '+sh
    ;
    if (wf(0).id eq 177) or (wf(0).id eq 185) then $
      for i=0,6 do begin
         h=histogram((long(wf.d(i,*)))(1:nwf*144l-1)- (long(wf.d(i,*)))(0:nwf*144l-2),bin=1l,min=-65535l,max=65535l)
         ih=where(h ne 0,nh) & sh='' & for j=0,nh-1 do sh=sh+string(ih(j)-65535l,form='(i7)')
         printf,1,'   Channel '+string(i,form='(i1)')+' data increments: '+sh
      end
    ;
    ip=where(  dif ne pak200   ,np)
    if np gt 0 then ip=[-1,ip,nwf-1] else  ip=[-1,nwf-1]
    lp=ip(1:np+1)-ip(0:np)
    ;stop
    printf,1
    printf,1,'   WF capture start time: '+t0wf
    printf,1,'   Number of continuous data intervals: '+strtrim(string(np+1),2)
    printf,1,'   List of continuous data intervals (DI): '
    printf,1,'   DI No.  Duration (ms)  Packets   Gap to the next (ms)'+$
            '  Aver. 1 (TMU)  Aver. 2 (TMU)  Aver. 3 (TMU)  Aver. 4 (TMU)  Aver. 5 (TMU)  Aver. 6 (TMU)  Aver. 7 (TMU)'+$
            '  SDev. 1 (TMU)  SDev. 2 (TMU)  SDev. 3 (TMU)  SDev. 4 (TMU)  SDev. 5 (TMU)  SDev. 6 (TMU)  SDev. 7 (TMU)'
    for i=0,np do begin
       std=fltarr(7) & for j=0,6 do std(j)= stddev( wf(ip(i)+1:ip(i+1)).d(j,*) )
       ofs=fltarr(7) & for j=0,6 do ofs(j)= total(wf(ip(i)+1:ip(i+1)).d(j,*))/lp(i)/144.
       minv=intarr(7) & for j=0,6 do minv(j)= min( wf(ip(i)+1:ip(i+1)).d(j,*) )
       maxv=intarr(7) & for j=0,6 do maxv(j)= max( wf(ip(i)+1:ip(i+1)).d(j,*) )
       if (np gt 0)and(i lt np) then gap=(dif(ip(i+1))-pak200)/200. else gap=0.
       printf,1,i, lp(i)*pak200/200., lp(i),gap , ofs,std,$
                  form='(i9,f15.3,i9,f23.3,14f15.3)'
       printf,1
       printf,1,i, minv, maxv, form='(14i6)'
    end
    ;
    printf,1
    err=0
    openw,2,'crcfile.txt'
    printf,2,file
    close,2
    spawn,'crc'
    openr,2,'crcfile.txt'
    readf,2,file1
    crc=replicate(0u,nwf)
    readf,2,crc,form='(z4.4)'
    close,2
    for i=0l,nwf-1 do $
;      if crcwf(wf(i)) ne wf(i).crc then begin
      if crc(i) ne wf(i).crc then begin
        err=1
        printf,1,'   SRAM CRC error in packet no. '+strtrim(string(i),2)+': stored CRC='+$
                  string(wf(i).crc,form='(z4.4)')+' actual CRC='+string(crcwf(wf(i)),form='(z4.4)')+$
                  ' in SRAM no. '+string((wf(i).adr and 'f000'xu)/'1000'xu, form='(i2)')+$
                  ' at address '+string((wf(i).adr and 'fff'xu), form='(z3.3)')
      end
    if err eq 0 then printf,1,'   No SRAM CRC errors found. '
 close,1
end



pro plotwf,wf,fs,file,title, ip,np, t0wf,donotclose=donotclose,sub=sub,logf=logf
;
; plots E+B 7 channel waveforms+spectra from wf
; sub plots also subintervals of 5000 points
; if sub=-1 plot all subintervals
;
  psfile=strmid(file,0,strpos(file,'.',/reverse_s))+'.ps'
  print,'----    Plotting waveform figure file: ',psfile
  if not keyword_set(donotclose) then plot_init, psfile
    nt=n_elements(wf.d)/7
    wfm=reform(wf.d,7,nt)
  next_page
  plot_wf_mult,7,wfm,fs, 0.1,0.75,   ip=ip,np=np,tit='ELMAVAN  '+title+ $
           '!c !c!u Processed '+systime()+' by ELM ver.'+ver()+' from '+file+'!n',$
           xtit=' Time (s) from '+t0wf
  plot_spect_mult,7,wfm,fs,  0.82,0.995, tit='Power spectra',logf=logf
;
;---- subintervals
  subl=50000; length in samples
if (nt gt subl) and keyword_set(sub) then begin
  if sub eq -1  then maxsub=nt/subl else maxsub=sub-1
  for i=0l,maxsub do begin
     next_page
     plot_wf_mult,7,wfm,fs,  0.1,0.75, i*subl, (i+1)*subl-1,  ip=ip,np=np,$
      tit='ELMAVAN  '+title+'  Part '+$
           strtrim(string(i+1),2)+'/'+strtrim(string(nt/subl+1),2)+ $
           '!c !c!u Processed '+systime()+' by ELM ver.'+ver()+' from '+file+'!n',$
           xtit=' Time (s) from '+t0wf
     if (i+1)*subl-1 le nt then $
       plot_spect_mult,7,wfm, fs, 0.82,0.995, i*subl,(i+1)*subl-1, $
         tit='Power spectra',logf=logf
  end
end
  if not keyword_set(donotclose) then plot_close
end





pro plotg,wf,fs,file,title, ip,np, t0wf, donotclose=donotclose,sub=sub,logf=logf
;
; plots elmavan-g E+B 3 channel waveforms+spectra from wf
; sub plots also subintervals of 5000 points
; if sub=-1 plot all subintervals
;
;
  nt=n_elements(wf.d)/7
  wfm=reform(wf.d,7,nt)
;
; time -> UTC, t0,t1
  str_tns,'2012-10-01T08:00:00.000.000.000',tns_utc2 ;___ UTC+2 used beore this date
  str_tns,t0wf+'.000',tns0
  if tns0 lt tns_utc2 then begin
    tns0=ulong64(tns0-7200000000000ll)
    tns_str,tns0,t0wf
  end
  tns_str,tns0+ulong64(nt*1000000000ull/fs),t1wf
  ;
  psfile=strmid(file,0,strpos(file,'.',/reverse_s))+'.ps'
  print,'----    Plotting waveform figure file: ',psfile
  if not keyword_set(donotclose) then plot_init, psfile
  next_page
  plot_wf_mult,3,wfm,fs, 0.08,0.8,   ip=ip,np=np,tit='ELMAVAN  '+title+ $
           '!c !c!u Processed '+systime()+' by ELM ver.'+ver()+' from '+file+'!n',$
           /noparams,xtit=' Time (s) from '+strmid(t0wf,0,27),yexp=0.4
  plot_spect_mult,3,wfm,fs,  0.87,0.995, tit='Power spectra',logf=logf,yexp=0.4
  plot_spgram_mult,3,wfm,fs,  0.08,0.8 ,xtit='Time (s) from '+strmid(t0wf,0,27),$
    tit='',yexp=0.54,ytop=0.54,outspec=spec,outfreq=freq,npoi=4096,shft=2048;, logf=logf
  plot_params_mult,3,wfm,0.89,yexp=0.54,ytop=0.54, 0,nt-1,  avewf,sdewf,minwf,maxwf
;
;---- subintervals
  subl=50000; length in samples
if (nt gt subl) and keyword_set(sub) then begin
    wfmax=fltarr(nt/subl)
    for i=0l,nt/subl-1 do wfmax[i]=max(abs(wfm[*,i*subl:(i+1)*subl-1]))
  imax=rotate(sort(wfmax),2)
  if sub eq -1 then maxsub=nt/subl-1 else maxsub=sub-1
  for iii=0l,maxsub do begin
     if iii eq maxsub then i=imax[nt/subl-1] else i=imax[iii]
     next_page
     plot_wf_mult,3,wfm,fs,  0.08,0.8, i*subl, (i+1)*subl-1,  ip=ip,np=np,$
           /noparams,tit='ELMAVAN  '+title+'  Part '+$
           strtrim(string(i+1),2)+'/'+strtrim(string(nt/subl+1),2)+ $
           '!c !c!u Processed '+systime()+' by ELM ver.'+ver()+' from '+file+'!n',$
           xtit=' Time (s) from '+strmid(t0wf,0,27),yexp=0.4
     if (i+1)*subl-1 le nt then $
       plot_spect_mult,3,wfm, fs, 0.87,0.995, i*subl,(i+1)*subl-1, tit='Power spectra',logf=logf,yexp=0.4
     plot_spgram_mult,3,wfm,fs,  0.08,0.8, i*subl,(i+1)*subl-1 ,xtit='Time (s) from '+strmid(t0wf,0,27),$
       tit='',yexp=0.54,ytop=0.54;, logf=logf
     plot_params_mult,3,wfm,0.89,yexp=0.54,ytop=0.54, i*subl,(i+1)*subl-1
  end
end
if not keyword_set(donotclose) then plot_close
;
; web files
psfile='el.ps'
  print,'----    Plotting waveform figure file: ',psfile
  plot_init, psfile
  next_page
  plot_wf_mult,3,wfm,fs, 0.08,0.8,   ip=ip,np=np,tit='ELMAVAN  '+title+ $
           '!c !c!u Processed '+systime()+' by ELM ver.'+ver()+' from '+file+'!n',$
           /noparams,xtit=' Time (s) from '+strmid(t0wf,0,27),yexp=0.4
  plot_spect_mult,3,wfm,fs,  0.87,0.995, tit='Power spectra',logf=logf,yexp=0.4
  plot_spgram_mult,3,wfm,fs,  0.08,0.8 ,xtit='Time (s) from '+strmid(t0wf,0,27),$
    tit='',yexp=0.54,ytop=0.54,outspec=spec,outfreq=freq,npoi=4096,shft=2048;, logf=logf
  plot_params_mult,3,wfm,0.89,yexp=0.54,ytop=0.54, 0,nt-1,  avewf,sdewf,minwf,maxwf
plot_close
;---- subintervals
psfiles=['el1.ps','el10.ps']
maxsub=1
  for iii=0l,maxsub do begin
    print,'----    Plotting waveform figure file: ',psfiles[iii]
    plot_init, psfiles[iii]
     if iii eq maxsub then i=imax[nt/subl-1] else i=imax[iii]
     next_page
     plot_wf_mult,3,wfm,fs,  0.08,0.8, i*subl, (i+1)*subl-1,  ip=ip,np=np,$
           /noparams,tit='ELMAVAN  '+title+'  Part '+$
           strtrim(string(i+1),2)+'/'+strtrim(string(nt/subl+1),2)+ $
           '!c !c!u Processed '+systime()+' by ELM ver.'+ver()+' from '+file+'!n',$
           xtit=' Time (s) from '+strmid(t0wf,0,27),yexp=0.4
     if (i+1)*subl-1 le nt then $
       plot_spect_mult,3,wfm, fs, 0.87,0.995, i*subl,(i+1)*subl-1, tit='Power spectra',logf=logf,yexp=0.4
     plot_spgram_mult,3,wfm,fs,  0.08,0.8, i*subl,(i+1)*subl-1 ,xtit='Time (s) from '+strmid(t0wf,0,27),$
       tit='',yexp=0.54,ytop=0.54;, logf=logf
     plot_params_mult,3,wfm,0.89,yexp=0.54,ytop=0.54, i*subl,(i+1)*subl-1
    plot_close
  end
  ;
;__elm.txt -> blesklog.txt file--------------------
      openw,2,'elm.txt'
      printf,2,'   '+strmid(t0wf,11,8)+'-'+strmid(t1wf,11,8)+$
                   string(sdewf[0],sdewf[1],sdewf[2],form='(3f9.2)')
      close,2
    ;
;__accumulated overview spectrogram-----------------
  overviewspect_m,file, t0wf,t1wf,nt/fs/3600.,freq,spec, $
     savfile='elmavan_spectrogram.idl', psfile='elspg.ps',logy=0
end




pro overviewspect_m,file, ts0,ts1,durh,freq,spec, savfile=savfile, psfile=psfile,logy=logy
   nsp=3
;__accumulated overview figure-----------------
   maxtim=4ull*86400000000000ull ; 4 days,
   if keyword_set(ts0) then str_tns,ts0+'.000.000',tns0
   ;__ data need to exist
   prev= file_test(savfile,/read) and not file_test(savfile,/zero)
   ;__ data need to be a 2D matrix
   if prev then begin ;-------
       restore,savfile
       siz=size(spgram)
       if siz[0] ne 3 then prev=0
       if n_elements(tim) eq 0 then prev=0
       if siz[1] ne n_elements(tim) then prev=0
       if siz[2] ne n_elements(frq) then prev=0
   end
   ;__spectrogram data need to have the same number of frequencies as the current spectrum
   if prev and keyword_set(ts0) then begin
      if siz[2] ne n_elements(freq) then prev=0
      if n_elements(frq) ne n_elements(freq) then prev=0
   end
      ;__spectrogram data need to have the same frequencies as the current spectrum
   if prev and keyword_set(ts0) then begin
     for ifr=0,n_elements(freq)-1 do if freq[ifr] ne frq[ifr] then prev=0
   end
    ;__ not more time than maxtim
   if prev then begin
      if keyword_set(ts0) then tim=[tim,tns0]
      ntim=n_elements(tim)
      outt=0
      while tim[ntim-1]-tim[outt] gt maxtim do outt=outt+1
      tim=tim[outt:ntim-1] & ntim=ntim-outt
      if ntim lt 2 then prev=0
   end
   if prev then begin
      spgramn=fltarr(ntim,siz[2],nsp)
      spgramn[0:ntim-2,*,*]=spgram[outt:ntim+outt-2,*,*]
      if keyword_set(ts0) then spgramn[ntim-1,*,*]=spec $
                 else spgramn[ntim-1,*,*]=spgram[ntim+outt-1,*,*]
      spgram=spgramn
    ;__ no backward steps in time
     dtim=long64(tim[1:ntim-1]) - long64(tim[0:ntim-2])
     ibad=where(dtim le 0 ,nbad)
     if nbad gt 0 then prev=0
   end
   if prev then begin
      medtim=median(dtim,/doub)
      igap=where(dtim gt 2*medtim,ngap)
      if ngap gt 0 then begin
        gap=lonarr(ngap+1)
        for ig=1,ngap do gap[ig]=round(dtim[igap[ig-1]]/medtim)
        igap=[-1,igap,ntim-1]
      end else begin
        gap=[0]
        igap=[-1,ntim-1]
      end
      ;
      spgramn=replicate(!values.f_nan,siz[1]+1+total(gap),siz[2],nsp)
      for ig=0,ngap do $
        spgramn[igap[ig]+1+total(gap[0:ig]):igap[ig+1]+total(gap[0:ig]),*,*]=$
         spgram[igap[ig]+1:igap[ig+1],*,*]
      if keyword_set(ts0) then print,'----    Added spectrum at '+strmid(ts0,0,19)+' to the overview data.'
   end else begin
       if not keyword_set(ts0) then begin
         print,'----    Bad file '+savfile
         return
       end
       print,'----    First value for the overview data at '+strmid(ts0,0,19)
       print,'   --> '+savfile+' saved.'
       spgram=reform(spec,1,n_elements(freq),3) & spgramn=spgram & tim=tns0 & frq=freq
   end
; _saving data
   save,file=savfile,tim,frq,spgram
   if not keyword_set(ts0) then begin
     tns0=tim[ntim-1]
     tns_str,tns0,ts0 & ts1=ts0 & durh=0
   end
; _plotting
   spgfile=strmid(file,0,strlen(file)-4)+'_'+psfile
   plot_init,spgfile,yfrac=0.9
   !p.multi=0 &  !p.thick=2 & !x.thick=2 & !y.thick=2
   !p.charsize=1
   plot_spgram_mult,3,fltarr(3,3),1,  0.1,0.92, $
         tit='Overview power spectrograms till '+strmid(ts1,0,19),$
         xtit='Time (hours) from '+strmid(ts0,0,19),inspgram=spgramn,infreq=frq,$
         intint=[-float((tns0-tim[0])/1000000000ull)/3600.,durh]
   plot_close
   file_copy,spgfile,psfile,/overwrite
end






function crcwf,wf
 b=bytarr(2030)
 b(0)=wf.id & b(1)=wf.timh
 b(2)=wf.tim/'1000000'xul & b(3)=(wf.tim and 'ff0000'xul)/'10000'xul
 b(4)=(wf.tim and 'ff00'xul)/'100'xul & b(5)=(wf.tim and 'ff'xul)
 b(6)=wf.khz200H/'100'xul & b(7)=(wf.khz200H and 'ff'xul)
 b(8)=wf.khz200/'1000000'xul & b(9)=(wf.khz200 and 'ff0000'xul)/'10000'xul
 b(10)=(wf.khz200 and 'ff00'xul)/'100'xul & b(11)=(wf.khz200 and 'ff'xul)
 b(12 + 2*indgen(1008))= (wf.d and 'FF00'xul)/'100'xul
 b(12 + 2*indgen(1008)+1)= (wf.d and 'FF'xul)
 b(2028)=wf.adr/'100'xul & b(2029)=(wf.adr and 'ff'xul)
 return,crc16(b)
end

;{calculates CRC}
function CRC16, data ; 16-bit byte data array
  CRC='FFFF'xu;
  for i=0, n_elements(data)-1 do begin
    CRC= CRC XOR (data(i)*'100'xu)
    for j=1, 8 do begin
      if (CRC AND '8000'xu) eq '8000'xu $
      then CRC = ((CRC AND '7FFF'xu)*'2'xu) XOR '1021'xu $
      else CRC = CRC*'2'xu
    end
  end
  return,crc
end



; --------------------- Link test reading & diagnostics --------------------



pro rdlink, file,link,ok=ok
;
; reads the telemetry stream of Link test ELMAVAN data packets
;
; TM structure of link test data:
  datpac={id:0b, d:bytarr(2031)}
;--------------------------
  print,'----    Reading link test TM file: ',file
   openr,1,file
   s=fstat(1)
   n= s.size /2032L
  if n gt 0 then begin
    link=replicate(datpac,n)
    readu,1,link
    ok=1
  end else ok=0
  close,1
end


pro analink,link,file

  ; link test analysis
  ;
  n=n_elements(link)
  txtfile=strmid(file,0,strpos(file,'.',/reverse_s))+'.txt'
  print,'----    Writing link test analysis file: ',txtfile
  openw,1,txtfile
      printf,1,'--------------------------------------------------'
      printf,1,'-----  Resonance/ELMAVAN link test analysis  -----'
      printf,1,'-----  Generated on '          +systime()+ ' -----'
      printf,1,'--------------------------------------------------'
      printf,1,' ' &    printf,1,'   Source: '+file
      printf,1, n,' link test packets. '
    ;
    ok=1
    lastfirst=255b
    for i=0,n-1 do begin
      first=link(i).d(0)
      if first ne lastfirst+1b then $
          printf,1,'      Link test discontinuity: packet no.'+string(i+1)+' - 1st byte is '+$
                     string(first,form='(i3)')+' but should be '+string(lastfirst+1b,form='(i3)')
      dif=link(i).d(1:2030)-link(i).d(0:2029)
      ier=where(dif ne 1,ner)
      if ner gt 0 then  begin
        for k=0,ner-1 do $
          printf,1,'      Link test error: packet no.'+string(i+1)+' byte '+string(ier(k)+1)+' is '+$
                     string(link(i).d(ier(k)+1),form='(i3)')+' but should be '+string(link(i).d(ier(k))+1b,form='(i3)')
        ok=0
      end
      lastfirst=first
    end
       printf,1
       printf,1,'-------------------------------------------------'
    if ok then begin
       printf,1,'-----              Link test OK             -----'
    end else begin
       printf,1,'-----          Link test NOT OK             -----'
    end
       printf,1,'-------------------------------------------------'
  close,1
end





;-------------- plotting subroutines-------------------------------------



pro plot_params_mult,mult,wfm,x0,yexp=yexp,ytop=ytop, i1,i2,  avewf,sdewf,minwf,maxwf
  avewf=fltarr(mult) & sdewf=fltarr(mult) & minwf=fltarr(mult) & maxwf=fltarr(mult)
  for i=0,mult-1 do begin
    avewf[i]=(moment(wfm(i,i1:i2),sdev=sdewfi))(0) & sdewf[i]=sdewfi
    maxwf[i]=max(wfm(i,i1:i2),min=minwfi) & minwf[i]=minwfi
    xyouts,x0,ytop-0.005-yexp*i*0.91/mult,/norm,color=255,$
      '!c !c Channel '+strtrim(string(i),2)+$
      '!c !c mn!3: '+strtrim(string(minwf[i],form='(i10)'),2)+$
      '!c !c mx!3: '+strtrim(string(maxwf[i],form='(i10)'),2)+$
      '!c !c !Mm!3: '+strtrim(string(avewf[i],form='(f10.1)'),2)+$
      '!c !c !Ms!3: '+strtrim(string(sdewf[i],form='(f10.1)'),2)
  end
end

pro plot_wf_mult,mult,wfm, fs, x0,x1, i1,i2, ip=ip,np=np,tit=tit,xtit=xtit,$
                 yexp=yexp,ytop=ytop,noparams=noparams
;
;  plots a set of waveforms of mult channels
;
    if not keyword_set(yexp) then yexp=1 ; expansion/compression factor to make smaller plots
    if not keyword_set(ytop) then ytop=0.95; top vertical coordinate to place the plots
    noisescale=35. ; min y axis range
  ;
  ; plot all data if not selected
  n=n_elements(wfm(0,*))
  if not keyword_set(i2) then i2=n-1
  if not keyword_set(i1) then i1=0l

  for chan=0,mult-1 do begin
    if keyword_set(noparams) then params='' else begin
      std=stddev(wfm(chan,i1:i2))
      ave=total(wfm(chan,i1:i2))/float(i2-i1+1)
      params='!c !c !Mm!3: '+strtrim(string(ave,form='(f10.1)'),2)+$
             '!c !c !Ms!3: '+strtrim(string(std,form='(f10.1)'),2)
    end
    if chan eq 0 then toptit=tit else tit=''
    plot_wf,wfm(chan,*),fs,i1,i2,ytit='Channel '+strtrim(string(chan),2)+$
                 params+'!c !c!3TM units',$
            pos=[x0,ytop-yexp*(chan+1)*0.91/mult,x1,ytop-0.005-yexp*chan*0.91/mult],$
             tit=tit,noisescale=noisescale, nox=(chan lt mult-1),xtit=xtit ; yran=yran
  end
   if keyword_set(np) and keyword_set(ip) then begin
     for i=1,np do if (ip(i) gt i1) and (ip(i) lt i2) $
       then plots,/norm,replicate(x0+(x1-x0)/float(i2-i1)*((ip(i)+1l)*144.-0.5-i1),2),$
                  ytop-[yexp*0.91,0.005],line=1,col=255
   end
end


pro plot_spect_mult,mult,wfm,fs,  x0,x1, i1,i2,tit=tit,logf=logf,yexp=yexp,ytop=ytop
;
; plots a set of spectra of mult channels
;
    if not keyword_set(yexp) then yexp=1 ; expansion/compression factor to make smaller plots
    if not keyword_set(ytop) then ytop=0.95; top vertical coordinate to place the plots
   for chan=0,mult-1 do begin
    if chan eq 0 then toptit=tit else toptit=''
    plot_spect,wfm(chan,*),fs,i1,i2,tit=toptit,$ ;yran=yran
            pos=[x0,ytop-yexp*(chan+1)*0.91/mult,x1,ytop-0.005-yexp*chan*0.91/mult],$
            nox=(chan lt mult-1),logf=logf
  end
end


pro plot_spgram_mult,mult,wfm,fs,  x0,x1, i1,i2,tit=tit,logf=logf,yexp=yexp,ytop=ytop, xtit=xtit,$
                     outspec=spec,outfreq=freq,inspgram=inspgram,infreq=infreq,intint=tint,npoi=npoi,shft=shft
;
; plots a set of spectra of mult channels
;
    if not keyword_set(yexp) then yexp=0.95 ; expansion/compression factor to make smaller plots
    if not keyword_set(ytop) then ytop=0.95; top vertical coordinate to place the plots
   for chan=0,mult-1 do begin
    if not keyword_set(inspgram) then $
       get_spectrogram,wfm(chan,*),fs,i1,i2,logf=logf,spgram,freq,tint,npoi=npoi,shft=shft $
     else begin
        spgram=inspgram(*,*,chan) & freq=infreq
     end
    if chan eq 0 then begin
      toptit=tit
      spec=fltarr(n_elements(freq),mult)
    end  else toptit=''
    spec(*,chan)=total(spgram,1)/n_elements(spgram[*,0])
    plot_spectrogram,spgram,tint,freq/1e3,tit=toptit,logy=logf,xtit=xtit,$
            ytit='Channel '+strtrim(string(chan),2)+'!c !c Frequency (kHz)',$
            pos=[x0,ytop-yexp*(chan+1)*0.91/mult,x1,ytop-0.005-yexp*chan*0.91/mult],$
            nox=(chan lt mult-1)
  end
end


;---------------------------------------------------------------
;--------     general time routines, nanosecond MJD    ---------
;---------------------------------------------------------------



function into_tns,nanosecond,second,minute,hour,day,month,year
 ;
 ; composition of modified jul. date in nanoseconds,
 ;   8byte unsigned integer, nanoseconds from 17nov1858 till 6jun2443
 ;
  return,nanosecond +1000000000ull* (second+60ull*(minute+60ull*(hour + $
                  24ull*(julday(Month, Day, Year)- 2400001ull) ) ) )
end


pro from_tns, tns, nanosecond,second,minute,hour,day,month,year
 ;
 ; decomposition of modified jul. date in nanoseconds,
 ;   8byte unsigned integer, nanoseconds from 17nov1858 till 6jun2443
 ;
  CALDAT, long(tns    /86400000000000ull +2400001ull), Month , Day, Year
  nanosecond=ulong( tns mod 1000000000ull)
  second    =byte( (tns/    1000000000ull) mod 60ull)
  minute    =byte( (tns/   60000000000ull) mod 60ull)
  hour      =byte( (tns/ 3600000000000ull) mod 24ull)
end


pro str_tns, str,tns
 ;
 ; string -> modified jul. date in nanoseconds in 8byte unsigned integer
 ;
  n=n_elements(str)
  Year=uintarr(n) & Month=bytarr(n) & Day=bytarr(n)
  Hour=bytarr(n) & Minute=bytarr(n) & Second=bytarr(n) & nanosecond=ulonarr(n)
  Yeari=0ul &  Monthi=0ul &  Dayi=0ul &  Houri=0ul &  Minutei=0ul
  Secondi=0ul &  millisecondi=0ul &  microsecondi=0ul & nanosecondi=0ul
  for i=0,n-1 do begin
     reads,str(i),Yeari, Monthi, Dayi, Houri, Minutei, Secondi, millisecondi, microsecondi,nanosecondi, $
          form='(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2,1x,i3,1x,i3,1x,i3)'
     Year(i)=Yeari & Month(i)=Monthi & Day(i)=Dayi
     Hour(i)=Houri & Minute(i)=Minutei & Second(i)= Secondi
     nanosecond(i)=nanosecondi+1000ul*microsecondi+1000000ul*millisecondi
  end
  tns=into_tns(nanosecond,second,minute,hour,day,month,year)
end


pro tns_str,  tns,str
 ;
 ;  modified jul. date in nanoseconds in 8byte unsigned integer -> string
 ;
   from_tns, tns, nanosecond,second,minute,hour,day,month,year
  n=n_elements(tns)
  str=strarr(n)
  for i=0l,n-1 do $
    str(i)=string(Year(i), Month(i), Day(i), Hour(i), Minute(i), Second(i), $
                  nanosecond(i)/1000000ul, (nanosecond(i) mod 1000000ul)/1000ul,$
                  nanosecond(i) mod 1000ul,$
         form='(i4.4,"-",i2.2,"-",i2.2,"T",i2.2,":",i2.2,":",i2.2,".",i3.3,".",i3.3,".",i3.3)')
end





;---------------------------------------------------------------
;--------------     general plotting routines    ---------------
;---------------------------------------------------------------







pro plot_init,file,yfrac=yfrac,portrait=portrait
 ;
 ; open lanscape ps plot
 ; close by device,/close
 ;
 set_plot,'ps'
 if not keyword_set(yfrac) then yfrac=1.
  if keyword_set(portrait) then landscape=0 else landscape=1
 device,file=file,xsize=25.,ysize=18.*yfrac,landscape=landscape,xoffset=1.,yoffset=28.,$
             bits_per_pixel=8,/color, /helvetica
 !p.font=0
 !p.charsize=0.6
 !p.thick=3
 !x.thick=3
 !y.thick=3
 bw=0
 color_palette, bw, rp,gp,bp, colors
 tvlct,rp,gp,bp
end


pro plot_close
;
; closes ps file, restores normal crt plotting
 !p.multi=0
 device,/close
 if !VERSION.OS_FAMILY eq 'Windows' then set_plot,'win' else set_plot,'x'
 !p.font=1
 !p.thick=1
 !x.thick=1
 !y.thick=1

end



pro next_page
  plot,[0,1],/nodata,xst=4,yst=4
end



function linticks,axis,index,value
 return,strtrim(string(value,form='(i10)'),2)
end


pro plot_wf,wf,fs,i1,i2, ytit=ytit, tit=tit, pos=pos,yran=yran, noisescale=noisescale, nox=nox,xtit=xtitle
;
; plots a  waveform wf
;  selects data from index i1 to index i2 (if absent plots all data)
; keywords: ytit - y title
;           tit - top title
;           yran - y axis range, autoscaled if absent
;           nox -  if not set:  x axis in ms, annotated
;                  by 't (milliseconds) from '+time
;           noisescale - if set, this minimum scale is used when autoscale
;           pos
;
  n=n_elements(wf)
  y=wf
  ;
  ; plot all data if not selected
  if not keyword_set(i2) then i2=n-1
  if not keyword_set(i1) then i1=0l
  ;
  ; plot empty space if wanting more
    if i1 lt 0 then begin
       y=[replicate(!values.f_nan,-i1),reform(y)]
       i2=i2-i1 & n=n-i1 & i1=0
      end
    if i2 gt n-1 then y=[reform(y),replicate(!values.f_nan,i2-n+1)]
  ;
  ; axis and plot annotations
  !y.ticklen=-0.005/(pos(2)-pos(0))
  !x.ticklen=-0.005/(pos(3)-pos(1))
  x=(dindgen(i2-i1+1)+i1)/fs
  if nox then begin
      xtit=''
      xtickn=replicate(' ',30)
  end else begin
      if not keyword_set(xtitle) then xtit='t (seconds)' else xtit=xtitle
      xtickn=''
  end
  if not keyword_set(ytit) then ytit='y (a.u.)'
  if not keyword_set(tit) then tit=''
  ;
  y= y(i1:i2)  ;  copy the input data
  ;
  ; compress if more points than resolution
   n=n_elements(y)
   maxpoi=2000L ; must be even and close to total number of readable pixels
   if n gt maxpoi*3 then begin ; compress for ratios > 6 -> shorter data gaps on the sides
     com=n/maxpoi*2; compression factor
     cent=(n-com*maxpoi/2)/2 ; centering the compressed waveform
     maxy=max(reform(y(cent:com*maxpoi/2-1+cent),com,maxpoi/2),dim=1,min=miny,/nan)
     xc=rebin(x(cent:com*maxpoi/2-1+cent),maxpoi,/sample)
     y=fltarr(maxpoi) & y(lindgen(maxpoi/2)*2)=miny & y(lindgen(maxpoi/2)*2+1)=maxy
     if com*maxpoi/2 lt n then begin ; data gaps on the sides
       y=[!values.f_nan, y,!values.f_nan ]
       xc=[x(0),xc,x(n-1)]
       n=maxpoi+2
     end else n=maxpoi
     x=xc
     compr='!u   .!n'
   end else compr=''
  ;
  ; plotting
  if (not keyword_set(yran)) and keyword_set(noisescale) then begin
   totmaxy=max(y, min=totminy,/nan)
   yran=[(totminy+totmaxy)/2.-noisescale<totminy*1.05,$
         (totminy+totmaxy)/2.+noisescale>totmaxy*1.05 ]
  end
   plot,x,y,/noerase, ystyle=17, yrange=yran,ytitle=ytit+compr,xtitle=xtit,$
        ytickform='linticks',xtickn=xtickn,title=tit,xstyle=1,thick=2,pos=pos,col=255
end


pro logav,x,y,n,xl,yl,nl
 ; average A*q^i points
 ;
  q=(x(n-1)/x(0))^(1./(nl-1)) ;
  i1=long(alog10(1./n/(1-q)*(1-q^nl))/alog10(q)) ; leave as is for indices below i1
  nk=round((n-i1)*(1-q)/(1-q^(nl-i1))*q^lindgen(nl-i1))  ; numbers of data to average for indices above i1
  nk(nl-i1-1)=nk(nl-i1-1)-total(nk)+n-i1 ; adjustement for rounding
  ni=[replicate(1,i1),nk] ;  numbers of data to average
  ia=[0,total(ni,/cumul)] ; indices to data for averaging
  xl=fltarr(nl) & yl=xl
  for i=0,nl-1 do begin
    xl(i)=exp(total(alog(x(ia(i):ia(i+1)-1)))/ni(i))
    yl(i)=total(y(ia(i):ia(i+1)-1))/ni(i)
  end
end


function logticks,axis,index,value
 return,'10!u'+strtrim(string(alog10(value),form='(i10)'),2)+'!n'
end



pro plot_spect,wf,fs,i1,i2, ytit=ytit,tit=tit,yran=yran,xran=xran,psd=psd, $
         nox=nox,pos=pos,auxtit=auxtit,logf=logf
;
; plots a power spectrum of a waveform wf
;  selects data from index i1 to index i2 (if absent plots all data)
; keywords: ytit - y title
;           tit - top title
;           yran - y axis range, autoscaled if absent
;           nox -  if not set:  x axis  annotated
;           pos
  ;
  ; plot all data if not selected
  if not keyword_set(i2) then i2=n_elements(wf)-1
  if not keyword_set(i1) then i1=0
   n=i2-i1+1

  ; axis and plot annotations
  fscale=1 ; 1 for Hz, 1e-3 for kHz or 1e-6 for MHz
  if nox then begin
      xtit=''
      xtickn=replicate(' ',30)
  end else begin
      xtit='Frequency (Hz)'
      xtickn=''
  end
   if not keyword_set(ytit) then psdtit='Power (dB)' else psdtit=ytit

  if keyword_set(auxtit) then psdtit=auxtit+'!c!c'+psdtit
  if not keyword_set(tit) then tit=''
  !y.ticklen=-0.005/(pos(2)-pos(0))
  !x.ticklen=-0.005/(pos(3)-pos(1))
  ;
  ; substraction of the DC offset
  y= wf(i1:i2)-total(wf(i1:i2))/float(n)
  ;
  ;spectrum params
  if keyword_set(logf) then begin;  (Hanning window, max fft length, overlap forced to fit data with 3 ffts)
   nsmooth=1; moving average
   n2= 2L^long(alog(n)/alog(2))   ; final no of samples in a fft - power of 2
   nspec=n2/2    ; final number of points in the averaged spectrum
   nn=3 ; final no of averaged ffts
  end else begin
   npoi=1024 ; final frequency average of the spectrum to max npoi points (must be power of 2) +  moving average
   nfft=16 ;   no of averaged ffts will be from nfft to 2*nfft,
           ; nfft must be power of 2  (Hanning window, approx 50% overlap)
   nsmooth=1; moving average
   ;
   n2= 2L^long(alog(n)/alog(2))/nfft*2
           ; final no of samples in a fft - power of 2
   nspec=min([npoi, n2/2])
           ; final number of points in the averaged spectrum
   nn=ceil(float(n)/n2)*2-1 ; final no of averaged ffts
  end
   ovl=double(n2*nn - n)/double(nn-1) ; no of samples to overlap
   ;
     ; spectrum
   s=fltarr(n2/2)
   for i=0L,nn-1L do s=s+ 16./3.*$ ; correction factor for hanning window
               n2/fs/float(nn)*abs( ( fft(hanning(n2)*$
                y(i*(n2-ovl):i*(n2-ovl)+n2-1)) )(0:n2/2-1) )^2
   s=smooth( (rebin(s,nspec))(1:nspec-1) ,nsmooth, /EDGE_TRUNCATE, /NAN)
   x=fscale*(findgen(nspec-1)+1)/float(nspec-1)*fs/2.
  ;
   if keyword_set(logf) then begin
     nl=200 < nspec/2
     logav,x,s,nspec-1,xl,sl,nl
     x=xl &s=sl
   end
   if keyword_set(logf) and not nox then xform='logticks' else xform=''
   ;
    ierr=where(s le 0.,nerr)
  if nerr gt 0 then s(ierr)=1e-16
 ;
  if not keyword_set(xran) then xran=[min(x),max(x)]
  i=where( (x ge xran[0])and(x le xran[1]),ni)
  if (ni gt 0) then begin
     freq=x[i]
     if keyword_set(psd) then begin
       spec=s[i]
       yform='logticks'
     end else begin
       spec=10.*alog10(s[i]*fs/2.)+3.01
       yform=''
     end

  ; plotting
   plot,freq,spec, /noerase,yrange=yran,xtickform=xform,$
     ylog=keyword_set(psd),ytickform=yform,$
     ytitle=psdtit,xtitle=xtit,xtickn=xtickn, $
     pos=pos,title=tit,xstyle=1,xlog=logf,col=255
 end
end

pro get_spectrogram,wf,fs,i1,i2,logf=logf,spgram,freq,tint,psd=psd,npoi=npoi,shft=shft
;
; calculates a power spectrogram of a waveform wf
;  selects data from index i1 to index i2 (if absent plots all data)
; keywords:
  ;
  ; plot all data if not selected
  if not keyword_set(i2) then i2=n_elements(wf)-1
  if not keyword_set(i1) then i1=0
   n=i2-i1+1

   ; substraction of the DC offset
  y= wf(i1:i2)-total(wf(i1:i2))/float(n)
  ;
  ;spectrum params
  if keyword_set(logf) then begin;
   if not keyword_set(npoi) then npoi=1024
   nn=1
  end else begin
    if not keyword_set(npoi) then npoi=256
   nn=1
  end
  if not keyword_set(shft) then shft=16 ; no of samples to shift the fft window - overlap=npoi-shft
   ;
     ; spectrum
  nt=(n-npoi)/(shft*nn)
  spgram=fltarr(nt,npoi/2)
  for it=0L,nt-1 do $
   for i=0L,nn-1L do $
      spgram[it,*]=spgram[it,*]+ 16./3.*$ ; correction factor for hanning window
               npoi/fs/float(nn)*abs( ( fft(hanning(npoi)*$
                y[(it*nn+i)*shft:(it*nn+i)*shft+npoi-1]) )(0:npoi/2-1) )^2
  freq=(findgen(npoi/2)+1)/float(npoi/2)*fs/2.
  tint=[0.,((nt*nn-1)*shft+npoi)/fs]+i1/fs
  ;
  if keyword_set(logf) then begin
     nl=200 < npoi/4
     spgraml=fltarr(nt,nl)
     for it=0L,nt-1 do begin
        logav,x,spgram[it,*],npoi/2,xl,sl,nl
        spgraml[it,*]=sl
     end
     spgram=spgraml
     freq=xl
  end
;
  if not keyword_set(psd) then spgram=10.*alog10(spgram*fs/2.)+3.01
end

pro plot_spectrogram,sp,x,y,xtit=xtit,ytit=ytit,ztit=ztit,tit=tit,$
   zran=zran,logy=logy,logz=logz,nox=nox,pos=pos
;
; plots a power spectrogram of a matrix sp
;           ytit - y title
;           tit - top title
;           yran - y axis range
  ; axis and plot annotations
  if not keyword_set(tit) then tit=''
  ;
  if not keyword_set(ztit) then psdtit='Power (dB)' else psdtit=ztit
  if not keyword_set(logz) then logz=0 else logz=1
;
  if not keyword_set(zran) then begin
    hhh=histogram(sp,nbins=200,loc=zzz,/nan)
    hcum=total(hhh,/cumulative)/float(total(hhh))
    zran= [zzz[(where(hcum gt 0.1))(0)],zzz[(where(hcum gt 0.995))(0)] ]
  end
  ;
  if not keyword_set(logy) then logy=0 else logy=1
  if not keyword_set(ytit) then ytit='Frequency (Hz)'
  ;
 bottom=64
 if not keyword_set(pos) then pos=[0.2,0.21,0.9,0.95]
 posi=[pos[0],pos[1],pos[2],pos[3]]
 poss=[pos[2]+0.065,pos[1],pos[2]+0.07,pos[3]]
 ;
 ;
 if logy then begin
    v=10^(alog10(min(y))+alog10(max(y)/min(y))/(n_elements(y)-1)*findgen(n_elements(y)))
   yticks=floor(alog10(max(y)))-ceil(alog10(min(y)))
    u=ceil(alog10(min(y)))+findgen(yticks+1)
   ytickv=interpol(v,y,10^u)
   ytickname='10!u'+strtrim(string(u,form='(i10)'),2)+'!n'
   um=0
   for i=u(0)-1,u(yticks) do for j=2,9 do begin
     mt=j*10^i
     if (min(y) lt mt)and(max(y) gt mt) then um=[um,mt]
   end
   mticks=interpol(v,y,um[1:n_elements(um)-1])
 end else begin
   yticks=0
   ytickv=0
   ytickname=replicate('',30)
   mticks=0
 end
 ;
 if logz then begin
   scalez=alog10(zran)
   zform='logticks'
 end  else begin
   scalez=zran
   zform=''
 end
 !x.ticklen=-0.001 & !y.ticklen=-0.5
 plot,[0,1],zran, /noerase,/nodata,ylog=logz,$                        ;vykresli stupnici
      ytit=psdtit,ytickform=zform,xstyle=4,$
      pos=poss,ythick=3, xthick=3,ystyle=1,color=255
 TVImage,rotate(scalez(0)+findgen(255)/254.*(scalez(1)-scalez(0)),1),poss, scalez(0), scalez(1),$
         top=253,bottom=bottom,invalid=0
 oplot,[0,0,1,1,0],zran([0,1,1,0,0]), thick=3, color=255,/noclip
  ;
 if nox then begin
   xti=''
   xtickn=replicate(' ',30)
 end else begin
   if keyword_set(xtit) then xti=xtit else xti='Time (h)'
   xtickn=''
 end
  ;
 plot,/noerase,[min(x),max(x)],[min(y),max(y)],/nodata,$        ;vykresli graf
      ylog=logy,xst=4,yst=4,pos=posi
 TVImage,sp,posi, zran(0), zran(1),$
         top=253,bottom=bottom,invalid=0
  !y.ticklen=-0.005/(pos(2)-pos(0))
  !x.ticklen=-0.005/(pos(3)-pos(1))
 plot,/noerase,[min(x),max(x)],[min(y),max(y)],tit=tit,/nodata,$
      xtit=xti,xtickn=xtickn,ytit=ytit,ystyle=1,ylog=logy,$
      pos=posi,ythick=3, xthick=3,xstyle=1,color=255,$
      yticks=yticks,ytickv=ytickv,ytickname=ytickname
 if logy then $
   for i=0,n_elements(mticks)-1 do begin
     plots,min(x)+[(min(x)-max(x))*0.006,0.],[mticks[i],mticks[i]],color=255,thick=3
     plots,max(x)+[(max(x)-min(x))*0.006,0.],[mticks[i],mticks[i]],color=255,thick=3
   end
end


pro color_palette, bw, rp,gp,bp, colors
;+
; color_palette
;   defines  a color palette and several particular colors
; Nov, 29, 1999, O.Santolik, written
; Jun, 06, 2000, O.Santolik, bw scale corrected
;-
 if bw then begin
   rp=[255,round(findgen(63)/62.*250.),255-round(findgen(190)/189.*255.),85,0]
   gp=rp & bp=rp
 endif else begin
  rp1=[255,255,250,246,242,238,234,230,226,222,218,214,210,206,202,198,194,190,$
       186,182,178,174,170,165,161,157,153,149,145,141,137,133,129,125,121,117,$
       113,109,105,101, 97, 93, 89, 85, 80, 76, 72, 68, 64, 60, 56, 52, 48, 44,$
    40, 36, 32, 28, 24, 20, 16, 12,  8,  4,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0]
  rp2=[  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  4,  8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64,$
    68, 72, 76, 80, 85, 89, 93, 97,101,105,109,113,117,121,125,129,133,137,$
       141,145,149,153,157,161,165,170,174,178,182,186,190,194,198,202,206,210]
  rp3=[214,218,222,226,230,234,238,242,246,250,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,128,  0]
  rp=[rp1,rp2,rp3]
  gp1=[255,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  7, 15, 23, 31, 39, 47, 55,$
    63, 71, 79, 87, 95,103,111,119,127,135,143,151,159,167,175,183,191,199]
  gp2=[207,215,223,231,239,247,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255]
  gp3=[255,255,255,255,255,255,255,255,255,255,255,250,246,242,238,234,230,226,$
       222,218,214,210,206,202,198,194,190,186,182,178,174,170,165,161,157,153,$
       149,145,141,137,133,129,125,121,117,113,109,105,101, 97, 93, 89, 85, 80,$
    76, 72, 68, 64, 60, 56, 52, 48, 44, 40, 36, 32, 28, 24, 20, 16, 12,  8,$
     4,  0,128,  0]
  gp=[gp1,gp2,gp3]
  bp1=[255,  0,  4,  8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64,$
    68, 72, 76, 80, 85, 89, 93, 97,101,105,109,113,117,121,125,129,133,137,$
       141,145,149,153,157,161,165,170,174,178,182,186,190,194,198,202,206,210,$
       214,218,222,226,230,234,238,242,246,250,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255]
  bp2=[255,255,255,255,255,255,255,246,238,230,222,213,205,197,189,180,172,164,$
       156,148,139,131,123,115,106, 98, 90, 82, 74, 65, 57, 49, 41, 32, 24, 16,$
     8,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0]
  bp3=[  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,128,  0]
  bp=[bp1,bp2,bp3]
 endelse
 colors={white:0,black:255,grey:254,red:1,blue:64,$
     green:125,yellow:190,cyan:96,violet:33,orange:220}
end


;+
; NAME:
;       TVImage
;
; PURPOSE:
;      Universal imaging of an array
;
; CATEGORY:
;       Image.
;
; CALLING SEQUENCE:
;     TVImage, M , Pos , MMin, MMax
;
; INPUTS:
;    M....input array of single precision elements
;    Pos..position (normal coord.)
;    MMin..min value of M, elements below will be set to (bottom) color
;    MMax..max value of M, elements above will be set to (top) color
;          NaN elements (invalid data) will be set to (invalid) color.
;
; KEYWORD PARAMETERS:
;
;     MaxXPix,MaxYPix.. input: maximum number of pixels of a scalable device
;                    in each coordinate
;                (used only if the dimension of the input matrix is larger)
;                    default value: 5000
;    Inverse.. inversion of color indices
;
;    Top..     maximum color index, default !d.n_colors-1
;    Bottom..  minimum color index , default 0
;    Invalid.. color index for NaN elements, default (bottom)
;
; PROCEDURE:
;       Scales an image and plots it, device independent.
;
; MODIFICATION HISTORY:
;      written   O.S. V/94
;      added the keywords MaxDt, Top & Inverse O.S. XII/94
;      keywords bottom and invalid added,
;      change of behaviour for elements>mmax, O.Santolik XI/99
;      MaxDt changed to MaxXPix and MaxYPix, O.S. I/2000
;      any dimensions of the input array accepted, O.S. II/2000
;-
;
pro TVImage, M, Pos, MMin, MMax, MaxXPix=MaxXPix, MaxYPix=MaxYPix,$
         Inverse=inve, Top=top,Bottom=bottom, Invalid=invalid
 ;on_error,2
 ;__keywords
 if NOT n_elements(top) then top=!d.n_colors-1
 if NOT n_elements(bottom) then bottom=0
 if NOT n_elements(invalid) then invalid=bottom
 inv =  keyword_set(inve)
 if top LE bottom then begin
    exch=bottom & bottom=top & top=exch
    inv= NOT inv
  end
 ;__scale to byte array
 i_nan=where(finite(M,/nan),n_nan)
 b=bottom+bytscl(M, Min=MMin, Max=MMax*(1.00433),top=top-bottom,/nan)
          ;od verze 5.2, drive 1.00005
 if n_nan gt 0 then b(i_nan)=invalid
 if keyword_set(inve) then b=!d.n_colors-1-b
 ;__getting a matrix at least 2x2
   sz = size(b)
   if sz(0) eq 0 then b=replicate(b,2,2) $
   else if sz(0) eq 1 then b=rebin(b,sz(1),2,/sample) $
   else if sz(0) gt 2 then begin
        b=reform(b) & sz = size(b)
        b=rebin(b,sz(1),sz(2),/sample)
    endif
   sz = size(b)
   if sz(1) le 1 then b=rebin(b,2,sz(2),/sample) & sz = size(b)
   if sz(2) le 1 then b=rebin(b,sz(1),2,/sample)
 ;__scale to the destination rectangle
 if ((!D.flags and 1L) EQ 0)  then begin
  ; ...non-scalable device (screen)
     px= round([Pos(0),Pos(2)] * !d.x_vsize)
     py= round([Pos(1),Pos(3)] * !d.y_vsize)
     sx = float(px(1)-px(0));
     sy = float(py(1)-py(0));
     if (sx GT 0) AND (sy GT 0) then $
     TV, poly_2d(b, $
        [[0,0],[sz(1)/sx,0]], [[0,sz(2)/sy], $
        [0,0]], 0, sx, sy), px(0), py(0)
   endif else begin
   ; ...scalable device (post script printer)
     if N_elements(maxxpix) EQ 0 then maxxpix=5000. $
     else if maxxpix gt 2. then maxxpix=float(maxxpix) else maxxpix=2.
     if N_elements(maxypix) EQ 0 then maxypix=5000. $
     else if maxypix gt 2. then maxypix=float(maxypix) else maxypix=2.
     ;print,'MaxDt,Sz1,Sz2',maxxpix,maxypix,sz(1),SZ(2)
     if sz(1) gt maxxpix then begin
       b= poly_2d(b,[[0,0],[sz(1)/MAXxpix,0]], $
         [[0,1.],[0,0]],0,maxxpix,sz(2) )
       sz(1)=maxxpix
    end
     if sz(2) gt maxypix then begin
       b= poly_2d(b,[[0,0],[1.,0]], $
           [[0,sz(2)/maxypix],[0,0]],0,sz(1),maxypix )
       sz(2)=maxypix
       end
     XSize=Pos(2)-Pos(0)
     YSize=Pos(3)-Pos(1)
     if (!d.x_vsize*xsize gt 1.)and(!d.y_vsize*ysize gt 1.) then $
      TV, b, $
    Pos(0), Pos(1), XSize=XSize,$
    YSize=YSize, /norm
   endelse
end




;----------------------------------------------------------------aux procs
;----------------------------------------------------------------aux procs
;----------------------------------------------------------------aux procs
;----------------------------------------------------------------aux procs
;----------------------------------------------------------------aux procs
;----------------------------------------------------------------aux procs







;------------------------------------------- lemi tests -----------------------------

pro subtract, name,sub=sub
  ; subtracts
  ;
  tmfile=name
  rdwf, tmfile,wf,fs
  ; ___ plot data
  if n_elements(wf) gt 0 then begin
      wf.d(2,*)=(wf.d(0,*)-wf.d(1,*))/sqrt(2)
      plot_bdat,wf,fs,'SUB_'+tmfile,/logf,sub=sub
      plot_bsp,wf,fs,'SUBLOGSP_'+tmfile,/logf,sub=sub
      plot_bsp,wf,fs,'SUBLINSP_'+tmfile,sub=sub
  end
end




;------------------------------------------- b plots -----------------------------


pro plotbwf,wf,fs,file,donotclose=donotclose,sub=sub,logf=logf
;
; plots B 3-channel waveforms+spectra from dat
; sub plots also subintervals
;
  if not keyword_set(donotclose) then plot_init, strmid(file,0,strpos(file,'.',/reverse_s))+'.ps'
    nt=n_elements(wf.d)/7
    wfm=reform(wf.d(0:2,*),3,nt)
  next_page
  plot_wf_mult,3,wfm,fs, 0.1,0.65, tit='ELMAVAN Waveforms from '+file+' processed '+systime()
  plot_spect_mult,3,wfm,fs,  0.7,0.995, tit='Power spectra',logf=logf
;
;---- subintervals
if keyword_set(sub) then begin
  subl=5000; length in samples
  for i=0l,nt/subl do begin
     next_page
     plot_wf_mult,3,wfm,fs,  0.1,0.65, i*subl,(i+1)*subl-1, $
      tit='ELMAVAN Waveforms from '+file+' processed '+systime()+'  part '+$
           strtrim(string(i+1),2)+'/'+strtrim(string(nt/subl+1),2)
     if (i+1)*subl-1 le nt then $
       plot_spect_mult,3,wfm,fs,  0.7,0.995, i*subl,(i+1)*subl-1, tit='Power spectra',logf=logf
  end
end
  if not keyword_set(donotclose) then plot_close
end


pro plotbwfsp,wf,fs,file,donotclose=donotclose,sub=sub,logf=logf
;
; plots B 3-channel spectra from dat
; sub plots also subintervals
;
  if not keyword_set(donotclose) then plot_init, strmid(file,0,strpos(file,'.',/reverse_s))+'.ps'
    nt=n_elements(wf.d)/7
    wfm=reform(wf.d(0:2,*),3,nt)
  next_page
  plot_spect_mult,3,wfm, fs, 0.1,0.995,/wfparam,logf=logf, $
      tit='Power spectra of ELMAVAN waveforms from '+file+' processed '+systime()
;
;---- subintervals
if keyword_set(sub) then begin
  subl=5000; length in samples
  for i=0l,nt/subl do begin
     next_page
     if (i+1)*subl-1 le nt then $
       plot_spect_mult,3,wfm, fs, 0.1,0.995, i*subl,(i+1)*subl-1,/wfparam,logf=logf,$
               tit='Power spectra of ELMAVAN waveforms from '+file+' processed '+systime()+'  part '+$
           strtrim(string(i+1),2)+'/'+strtrim(string(nt/subl+1),2)
  end
end
  if not keyword_set(donotclose) then plot_close
end


