;--------------
; processing of handyscope files
;   compile: blesk.pro, rdh.pro
;   run: process, '*.bin' or process,'filename.bin'
; 21/02/23 interferometry removed 
; 21/02/25 improved parsing of the gps sentence, station to small figures and spectrogram, spectrogram from 10 kHz
; 21/02/26 fixed bug in latitude parsing, fixed LS processing, added "Temp:" readout as t2
; 21/02/27 getlines removed, temperature writeout and spectra titles improved, 'UT' added 
;-------------

pro process,files,plotLSsevan=plotLSsevan,nant=nant,azim=azim,$
  sm_plot=sm_plot, sm_fcentr_exclude_Hz=sm_fcentr_exclude_Hz, $
  sm_fbandw_exclude_Hz=sm_fbandw_exclude_Hz,$
  sm_interval_sec=sm_interval_sec,sm_thres_stdev=sm_thres_stdev,sm_maxfig=sm_maxfig
; files... file name(s) of the handyscope bin files
;          with full path and extension, wildcards allowed
;           (examples:'*.bin','/mnt/hpdelo/20130812_144419.bin')
;
  f=file_search(files,count=nf)
  print,'---'+strtrim(string(nf),2)+' files found for '+files
  for i=0,nf-1 do begin
    print,string(i+1,form='(i4)')+' '+f[i]+' ---'
    rdh,f[i],plotLSsevan=plotLSsevan,nant=_nant,azim=_azim,$
     sm_plot=sm_plot, sm_fcentr_exclude_Hz=sm_fcentr_exclude_Hz, $
     sm_fbandw_exclude_Hz=sm_fbandw_exclude_Hz,$
     sm_interval_sec=sm_interval_sec,sm_thres_stdev=sm_thres_stdev,sm_maxfig=sm_maxfig
  end
end


pro rdh,file,plotLSsevan=plotLSsevan,nant=_nant,azim=_azim,$
  sm_plot=sm_plot, sm_fcentr_exclude_Hz=sm_fcentr_exclude_Hz, $
  sm_fbandw_exclude_Hz=sm_fbandw_exclude_Hz,$
  sm_interval_sec=sm_interval_sec,sm_thres_stdev=sm_thres_stdev,sm_maxfig=sm_maxfig
  ;
  ; reads one handyscope  file and an associated auxilliary txt file
  ; outputs *_sp#.ps and *_overview#.ps
  ; copies them to sp#.ps and overview#.ps
  ; adds text information to overview.txt
  ; adds aux information necessary for generation of overview#.ps to overview#.idl
  ;
  ; sm_plot... keyword, if set then plot small waveforms
    ;    sm_interval_sec=1e-3 ; duration in seconds
    ;    sm_thres_stdev=5;  threshold for max abs dev in units of std
    ;    sm_maxfig=200 ; maximum number of small waveforms to be plotted
    ;    sm_fcentr_exclude_Hz [fc1,fc2,...fcn] excluded lines in integrated small waveforms
    ;    sm_fbandw_exclude_Hz [bw1,bw2,...bwn] bandwidths of excluded lines 
    ;     example RMC at Vestale 216 kHz+-10 kHz
    ;     sm_fcentr_exclude_Hz=[216e3], sm_fbandw_exclude_Hz=[20e3]
  ; plotLSsevan..keyword, if set then LS sevan plots are added
  ; nant ... number of antennas 
  ;      ... overrides the value found (or not found) in the lightning data description file
  ; azim ... [1..nant] azimuth (from the North toward East) from antenna to a 
  ;          positive downward vertical current causing maximum positive B  ....
  ;          azimuth of the ant plane, to the right when looking along dB/dt>0
  ;      ... overrides the values found (or not found) in the lightning data description file
  ;   
  ;   VE Rustrel 2 antennas heliport till 2013/08/30, without antialiasing filters (+7dB)
  ;       nant=2 & azim=[[49.,319.]   
  ;   VE Rustrel 2 antennas 29.4m interferometry since 2013/08/30 till 2013/09/26, w/o aa filters (+7dB)
  ;       nant=2 & azim=[65.,335] & sepazim= 335 degrees; interferometry separation azimuth 1->2
  ;   VE Rustrel 2 antennas 89m interferometry since 2013/09/26 till 2014/10/15, 
  ;       nant=2 & azim=[245., 335.] & sepazim= 343 degrees; interferometry separation azimuth 1->2
  ;   VE Rustrel 4 antennas interferometry since 2014/10/15
  ;       nant=4 & azim=[235., 344, 333., 237.]
  ;       
  ;   ER Ersa since 2015/10/09 1 antenna with BLESKA    
  ;  
  ;   LS Lomnicky stit since 2016/08/06
  ;       nant=4 & azim=[145.,55.] 
  ;       
  ;   TW Ter Wish (LOFAR) since 2018/09/24 till 2020/08/13
  ;       nant=2 & azim=[150.,66.]
  ;   TW Ter Wish (LOFAR) since 2020/08/13
  ;       nant=2 & azim=[330.,246.]
  ;
  ;   ML Milesovka since 2017/06/21: 50.5552972N, 13.9310831E, 809m n.m.,  
  ;      nant=2 & azim=[40.,310.]  ; questionable 
  ;
  ;   DL Dlouhá Louka since 2020/08/14:   50.6484261N, 13.6556014E, 885m n.m., 
  ;      nant=2 & azim=[298,208]
  ;
  ;   KO Kopisty since 2020/09/08:    50.5440539N, 13.6231053E, 248m n.m., 
  ;       nant=2 & azim=[100.,190.]
  ;
  ;   KR Bohosudov since 2021/02/09:    50.6809864N, 13.8731911E, 258m n.m., 
  ;      nant=2 & azim=[184., 274.]
  ; 
      durus=52 ; duration of trigger detail in us
      fhp=2.99e3 ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
           ; lower cutoff frequency for spectra
  ;
  ;  CALIBRATION SINCE 2013/09/26,  by 7 dB higher before that date
  ;           0=23dB,1=29dB,2=35dB,4=41dB,8=47dB,  9=30dB loGain;
      gaindB=[23.,29.,35.,0.,41.,0.,0.,0.,47.,30.]
      AntSurface = 0.228 ; = 0.435 * 0.525 m^2
  ;
  ; _ auxilliary txt file
   rdth,file,trigtimes,waitt, $; trig time string, wait time in seconds
              nsam, $  ; number of samples in each channel
              postsam, $; number of post-trigger samples in each channel
              terr,$ ; time accuracy
              t1,t2,gain,poe,fs,$; sampling fr in Hz
              bits,$
              range,$ ;fltarr[2] oscilloscope range in each channel +-range V
              trigl0,trigl1,$ ;fltarr[2] trig level in each channel, fraction of -range..range
              timeout, $;0 if triggered
              station,stationcode,$ ; name, 2 letter code
              pos,$ ;position from the GP GGA sentence
              ngps,$ ;number of satellites from the GP GGA sentence
              nant,$ ; number of antennas
              azim ; [1..nant] azimuth (from the North toward East) from antenna to a
              ;          positive downward vertical current causing maximum positive B  ....
              ;          azimuth of the ant plane, to the right when looking along dB/dt>0
  if keyword_set(_nant) then nant=_nant       
  if keyword_set(_azim) then azim=_azim
  ;
  ;_trigtime -> first sample
  str_tns, trigtimes,trigtime
  tns_str, trigtime-double(nsam-postsam)/fs*1d9, times
  ;
  ;_data reading
  nsamr=(file_info(file)).size /nant/2l
  ar=uintarr(nsamr,nant)
  openr,1,file &  readu,1,ar &  close,1
  wf= fltarr(nsam, nant) ; nsam-nsamr filled by zeros
  wf[0:(nsamr<nsam)-1,*]=float(ar[0:(nsamr<nsam)-1,*])-2.^(bits-1) & ar=0
  cal=fltarr(nant)
  for i=0,nant-1 do begin
   ;_calibration  A/D units -> nT/s;   nT/s= cal*(A/D +offs)
    cal[i]=range[i]/2d0^(bits-1)/AntSurface/1d1^(gaindB[gain]/20.)*1d9
  end
  ;
  ; sevan
  if keyword_set(plotLSsevan) then begin
    ; plotting Sevan at LS, setting nant=2
    for i=2,3 do plot_one_sevan, reform(wf[*,i]),i,file,fs,fhp,$
                    times,nsam,postsam,bits,trigl0[i],trigl1[i],timeout
    nant=2
  endif
  ;
  ;_antenna cycle
  intwf=fltarr(nsam,nant)
  for i=0,nant-1 do begin
   ; _plot each
    plot_one_antenna, reform(wf[*,i]),  cal[i], intwfi, sm_plot=sm_plot,$
          sm_fcentr_exclude_Hz=sm_fcentr_exclude_Hz, $
          sm_fbandw_exclude_Hz=sm_fbandw_exclude_Hz,$
          sm_interval_sec=sm_interval_sec,sm_thres_stdev=sm_thres_stdev,sm_maxfig=sm_maxfig,$
          i, azim[i], file, fs, fhp,times, durus,$ ; ant no, azimuth, file,...
              trigtimes,waitt,nsam,postsam,terr,$ ; configuration...
              t1,t2,gaindB[gain],poe,bits,range[i],trigl0[i],trigl1[i],timeout,$
              station,stationcode,pos,ngps
    intwf[*,i]=intwfi
  end
 ;
  textoverview,file,stationcode,wf,intwf,cal,times,waitt,terr,timeout
end


pro rdth,file,trigtimes,waitt,nsam,postsam,terr,$
              t1,t2,gain,poe,fs,bits,range,trigl0,trigl1,timeout,$
              station,stationcode,pos,ngps,nant,azim
  ;
  ; reads auxilliary txt file
  ;
  li=''
  times='2013-08-09T00:00:00.000.000.000'
  waitt='0ns' & nsam=33553920l & postsam=nsam
  terr='0ns' &  t1='0 C'&  t2='0 C' & gain=4b & poe='0V'
  fs=2e8 & bits=12 & range=[4.,4.,4.,4.]& trigl0=[0.,0.,0.,0.] & trigl1=[0.,0.,0.,0.]
  timeout=0
  station='' & stationcode='  ' & pos='' & ngps='' 
  nant=2 & azim=[0,0]
  channo=0 ; switching between channels
  ;
 txtfile=strmid(file,0,strlen(file)-4)+'.txt'
 if file_test(txtfile) then begin
  openr,1,txtfile
  while (not eof(1)) do begin
    readf,1,li
    if strmid(li,0,9) eq 'Station: ' then $
      station=strmid(li,9,strlen(li)-9)
    if strmid(li,0,14) eq 'Station code: ' then $
      stationcode=strmid(li,14,2)
    if strmid(li,0,4) eq 'GPS:' then begin
       gga=strsplit(li,',',/ex)
       if n_elements(gga) ge 10 then begin
         reads,strmid(gga[2],2,7),mnlat & reads,strmid(gga[4],3,7),mnlon
         pos=strmid(gga[2],0,2)+'.'+strtrim(string(round(mnlat/60.*1e6)),2)+gga[3]+', '+ $
             strmid(gga[4],0,3)+'.'+strtrim(string(round(mnlon/60.*1e6)),2)+gga[5]+', '+$
             gga[9]+'m'        
         ngps=gga[7]
       endif
    endif  
    if strmid(li,0,25) eq 'Number of SLAVIA sensors:' then begin
      reads,strmid(li,25,strlen(li)-25),nant
      azim=intarr(nant)
    endif  
    aux=0l
    for i=0,nant-1 do $
      if strmid(li,0,32) eq 'Loop azimuth to +CG of SLAVIA '+strtrim(string(i+1),2)+':' then begin
         reads,strmid(li,32,strlen(li)-32),aux
         azim[i]=aux
      endif   
    if strmid(li,0,15) eq 'Est time error ' then $
        terr=strmid(li,15,strlen(li)-15)
    if strmid(li,0,7) eq 'Core 0:' then $
        t2=strtrim(strmid(li,7,strlen(li)-7),2)
    if strmid(li,0,5) eq 'Temp:' then $
        t2=strtrim(strmid(li,5,strlen(li)-5),2)
    if strmid(li,0,6) eq 'temp2:' then $
        t2=strtrim(strmid(li,6,strlen(li)-6),2)
    ii=strpos(t2,'°C') & if ii ne -1 then t2=strmid(t2,0,ii)
    ii=strpos(t2,'° C') & if ii ne -1 then t2=strmid(t2,0,ii)
    ii=strpos(t2,'°') & if ii ne -1 then t2=strmid(t2,0,ii)
    if strmid(li,0,6) eq 'temp1:' then $
        t1=strtrim(strmid(li,6,strlen(li)-6),2)
    if strmid(li,0,12) eq 'gain set to '  then $
        reads,strmid(li,12,strlen(li)-12),gain
    if strmid(li,0,15) eq 'Preamp voltage:'  then $
        poe=strmid(li,15,strlen(li)-15)
    if strmid(li,0,12) eq 'PoE voltage:'  then $
        poe=strmid(li,12,strlen(li)-12)
    ;
    if strmid(li,0,11) eq 'Resolution:' then $
        reads,strmid(li,11,strlen(li)-11),bits
    if strmid(li,0,14) eq 'Sampling freq:'   then $
        reads,strmid(li,14,strlen(li)-14),fs
;
    if strmid(li,0,8) eq 'Channel '   then $
        reads,strmid(li,8,1),channo
    aux=0.
    if strmid(li,0,7) eq ' Range:' then begin
        reads,strmid(li, 7,strlen(li)-7),aux
        if channo gt 0 then range[channo-1]=aux
    end
    if strmid(li,0,15) eq '  Trig 0 level:' then begin
        reads,strmid(li,15,strlen(li)-15),aux
        if channo gt 0 then trigl0[channo-1]=aux
     end
    if strmid(li,0,15) eq '  Trig 1 level:' then begin
        reads,strmid(li,15,strlen(li)-15),aux
        if channo gt 0 then trigl1[channo-1]=aux
    end
    ;
    if strmid(li,0,12) eq 'Trigger at: ' then $
        trigtimes=strmid(li,12,23)+'.'+strmid(li,35,3)+'.000'
    if strmid(li,0,11) eq 'Waited for ' then $
        waitt=strmid(li,11,strlen(li)-11)
    if strmid(li,0,13) eq 'Postsamples: '  then $
        reads,strmid(li,13,strlen(li)-13),postsam
    if strmid(li,0,15) eq 'Valid samples: '  then $
        reads,strmid(li,15,strlen(li)-15),nsam
   if strmid(li,0,18) eq 'SW trigger timeout' then timeout=1
  end
  close,1
 end else print,'  ***'+txtfile+' NOT FOUND ***'
end



																		  
pro plot_one_sevan,wfi, i, file, fs, fhp, times, nsam, postsam, bits,  trigl0i, trigl1i, timeout
; plotting Sevan data from LS 
;
    psfile=strmid(file,0,strlen(file)-4)+'_sp'+strtrim(string(i),2)+'.ps'
    plot_init,psfile,yfrac=1
    !p.charsize=1. & !p.multi=0 &  !p.thick=2 & !x.thick=3 & !y.thick=3				 
     trigs=['Trigger','Timeout']
     miny=min(wfi,max=maxy) & yr=[miny,maxy]+[-1.,1.]*(maxy-miny)*0.1
     plot_wf,wfi,times+' UT',fs=fs,fhp=fhp,pos=[0.09,0.06,0.995,0.93],/noerase,$
                yran=yr,ytit=' Raw data [TM]',tit='SEVAN CH'+strtrim(string(i-1),2)
               ;'!d '+trigs[timeout]+' at '+$
               ;strtrim(string(round((trigl0i*2d0-1d0)*2d0^(bits-1))),2)+' and '+$
               ;strtrim(string(round((trigl1i*2d0-1d0)*2d0^(bits-1))),2)+' units, '+$
               ;strtrim(string(round(double(nsam-postsam)/fs*1d3)),2)+' ms!n'
        xyouts,/norm,al=1.,0.606,0.086,col=255,$
        '!c!mm!3 = '+strtrim(string(total(double(wfi))/double(nsam),form='(f30.2)'),2)+$
        ', !ms!3 = '+strtrim(string(stddev(wfi,/double),form='(f30.2)'),2)+' units  '
        xyouts,/norm,al=1.,0.836,0.064,col=255,$
         strtrim(string(min(wfi),form='(i30)'),2)+' to '+$
         strtrim(string(max(wfi),form='(i30)'),2)+' units  '		  		 
     plot_close
end



pro plot_one_antenna, wfi, cali, intwfi, $  ; data, calibration in, integrated out
          i, azimi, file, fs, fhp,times,durus,$ ; ant no, azimuth, file,...
              trigtimes,waitt,nsam,postsam,terr,$ ; configuration to report...
              t1,t2,gdB,poe,bits,rangei,trigl0i,trigl1i,timeout,$
              station,stationcode,pos,ngps,$
              sm_plot=sm_plot,sm_fcentr_exclude_Hz=sm_fcentr_exclude_Hz, $
              sm_fbandw_exclude_Hz=sm_fbandw_exclude_Hz, $
              sm_interval_sec=sm_interval_sec,sm_thres_stdev=sm_thres_stdev,sm_maxfig=sm_maxfig
              
    ;
    ; ps plot and overview spectrogram for one antenna
    ; outputs integrated waveform
    ; 
    ; if sm is set: small plots 
    ;    sm_interval_sec=1e-3 ; duration in seconds
    ;    sm_thres_stdev=5;  threshold for max abs dev in units of std
    ;    sm_maxfig=200 ; maximum number of small waveforms to be plotted
    ;    sm_fcentr_exclude_Hz [fc1,fc2,...fcn] excluded lines in integrated small waveforms
    ;    sm_fbandw_exclude_Hz [bw1,bw2,...bwn] bandwidths of excluded lines 
    ;
    psfile=strmid(file,0,strlen(file)-4)+'_sp'+strtrim(string(i),2)+'.ps'
    plot_init,psfile,yfrac=1
    !p.charsize=1. & !p.multi=0 &  !p.thick=2 & !x.thick=3 & !y.thick=3

    ;
    ;_title and global log spectrum
    trigs=['TRIGGER','SURVEY']
    plot_spect,/logf,fs=fs,wfi,pos=[0.09,0.6,0.995,0.93],$
               0,nsam-1, ytit='PSD [dB from ('+strtrim(string(round(cali/1d3)),2)+$
                               ' !mm!3T/s)!u2!n / '+strtrim(string(round(fs/2d6)),2)+'MHz]!c'+$
              ' !d 0dB...Noise with sigma at 1 dig. unit, Bandwidth = '+strtrim(string(round(fs/2d6)),2)+' MHz!n',$
               tit=stationcode+ strtrim(string(i+1),1)+' (loop azimuth '+strtrim(string(round(azimi)),2)+'!uo!n)  '+$
                    strtrim(string(long(nsam/fs*1e3)),2)+' ms from '+times+' UT after '+waitt+' '+trigs[timeout]+'            !c!d'+$
               ' !d'+file+', U='+strmid(poe,1,4)+$
               'V, T='+t2+', G='+strtrim(string(round(gdB)),2)+'dB, R='+$
               strtrim(string(round(rangei*1e3)),2)+'mV, N='+strtrim(string(nsam),2)+$
               ' at '+strtrim(string(fs/1e6,form='(f30.2)'),2)+'MHz, bits='+strtrim(string(bits),2)+$
				', GPS='+pos+', nGPS='+ngps+'     !n',col=255,freq=freq,spec=spec , xran=[fhp,fs/2.]
				smoo=smooth(spec,50) & smoo[0:315]=smoo[315] & smoo[748:898]=smoo[748]
				save,freq,spec,smoo,file='spec'+ strtrim(string(i+1),2)+'.sav'
    ;
    ;_integrated global waveform
    plot_wf,wfi*cali,times+' UT',fs=fs,fhp=fhp,out_wf=intwfi,$
               /int,out_yran=yri,ytit='Integrated data [nT]!c!d Highpass filtered at '+$
               strtrim(string(round(fhp/1d3)),2)+' kHz',$
               pos=[0.09,0.293,0.536,0.52],/noerase,/notimeaxis
       xyouts,/norm,al=1.,0.536,0.519,col=255,$
        '!c!ms!3 = '+strtrim(string(stddev(intwfi,/double),form='(f30.2)'),2)+' nT  '
       xyouts,/norm,al=1.,0.536,0.297,col=255,$
        strtrim(string(min(intwfi),form='(f30.2)'),2)+' to '+$
        strtrim(string(max(intwfi),form='(f30.2)'),2)+' nT  '
    ;
    ;_raw global waveform
     trigs=['Trigger','Timeout']
     miny=min(wfi,max=maxy) & yr=[miny,maxy]+[-1.,1.]*(maxy-miny)*0.1
     plot_wf,wfi,times+' UT',fs=fs,fhp=fhp,pos=[0.09,0.06,0.536,0.287],/noerase,$
                yran=yr,ytit=' Raw data ['+strtrim(string(round(cali/1d3)),2)+' !mm!3T/s]!c'+$
               '!d '+trigs[timeout]+' at '+$
               strtrim(string(round((trigl0i*2d0-1d0)*2d0^(bits-1))),2)+' and '+$
               strtrim(string(round((trigl1i*2d0-1d0)*2d0^(bits-1))),2)+' units, '+$
               strtrim(string(round(double(nsam-postsam)/fs*1d3)),2)+' ms!n'
        xyouts,/norm,al=1.,0.536,0.286,col=255,$
        '!c!mm!3 = '+strtrim(string(total(double(wfi))/double(nsam),form='(f30.2)'),2)+$
        ', !ms!3 = '+strtrim(string(stddev(wfi,/double),form='(f30.2)'),2)+' units  '
        xyouts,/norm,al=1.,0.536,0.064,col=255,$
         strtrim(string(min(wfi),form='(i30)'),2)+' to '+$
         strtrim(string(max(wfi),form='(i30)'),2)+' units  '
     ;
     ;_details around trigger
     plot_wf,reform(intwfi),trigtimes+' UT',/cent,$
                nsam-postsam-round(1d-6*durus*fs),nsam-postsam+round(1d-6*durus*fs), fs=fs,fhp=fhp,$
               yran=yri,pos=[0.54,0.293,0.995,0.52],/noerase,/notimeaxis,/noyaxis
     plot_wf,reform(wfi),trigtimes+' UT',/cent,yran=yr,$
                nsam-postsam-round(1d-6*durus*fs),nsam-postsam+round(1d-6*durus*fs),$
                fs=fs,fhp=fhp,pos=[0.54,0.06,0.995,0.287],/noerase,/noyaxis
     plot_close
     ;
     ; overview spectrogram (from blesk.pro)
     file_copy,psfile,'sp'+strtrim(string(i),2)+'.ps',/overwrite
     overviewspect,file, times,tunits=' UT',double(nsam)/fs/3600.,freq,spec, min_fr_Hz=1e4,max_fr_Hz=1e8,$
           source=stationcode+ strtrim(string(i+1),1)+' (loop azimuth '+strtrim(string(round(azimi)),2)+'!uo!n)  ',$
           scaletit='Power spectral density [dB from ('+strtrim(string(round(cali/1d3)),2)+$
           ' !mm!3T/s)!u2!n / '+strtrim(string(round(fs/2d6)),2)+'MHz]',$
           avetim=72ull  *60000000000ull, $  ;72 min, 20 columns per day  
           savfile='overview'+strtrim(string(i),2)+'.idl', $
           psfile='overview'+strtrim(string(i),2)+'.ps',/logy,trig=(timeout eq 0)
     ;
     ;-------------------------------
     ; plot small parts of the waveform, only for trigger if sm_plot is set
     ;
    if keyword_set(sm_plot) and (timeout eq 0) then begin 
      if not keyword_set(sm_interval_sec) then sm_interval_sec=1e-3 
        ; default: parts of 1 ms 
      if n_elements(sm_thres_stdev) eq 0 then sm_thres_stdev=5
        ; threshold for max abs dev in units of std
      if not keyword_set(sm_maxfig) then sm_maxfig=200
        ; maximum number of small waveforms to be plotted
      pal=long(fs*sm_interval_sec) ; number of samples in one part
      ;
      ; exclude frequencies from int waveform
      ; example RMC at Vestale 216 kHz+-10 kHz
      ;  sm_fcentr_exclude_Hz=[216e3], sm_fbandw_exclude_Hz=[20e3]
      nfex= n_elements(sm_fcentr_exclude_Hz)
      if n_elements(sm_fbandw_exclude_Hz) ne nfex then begin
         sm_fbandw_exclude_Hz=replicate(20e3,nfex)
         print,'------------------ inconsistent bandwidths for line exclusion in small plots, used 20 kHz'  
      end   
      sbandw='' & scentr='' & if nfex gt 0 then fexind=lonarr(2,nfex)
      for ifex=0,nfex-1 do begin
        fexind[0,ifex]=(sm_fcentr_exclude_Hz[ifex] -sm_fbandw_exclude_Hz[ifex]/2)*nsam/fs
        fexind[1,ifex]=(sm_fcentr_exclude_Hz[ifex] +sm_fbandw_exclude_Hz[ifex]/2)*nsam/fs   
        if ifex gt 0 then begin
          sbandw=sbandw+', ' & scentr=scentr+', '
        endif
        sbandw=sbandw+strtrim(string(round(sm_fbandw_exclude_Hz[ifex]/1e3)),2)
        scentr=scentr+strtrim(string(round(sm_fcentr_exclude_Hz[ifex]/1e3)),2)
      endfor
      if sbandw eq '' then sexcl='' else sexcl=sbandw+' kHz around '+scentr+' kHz excluded
      ;
      ; plotting
      psfile=strmid(file,0,strlen(file)-4)+'_sm'+strtrim(string(i),2)+'.ps'
      ;
      nsm=long(nsam/pal) ; number of parts 
      dsm=reform(wfi[0:nsm*pal-1],pal,nsm) ; decomposition to parts
      offss=median(dsm); over all parts 
      maxsm=max(abs(dsm-offss),dim=1,indsm); for each part separately
      thrsm=sm_thres_stdev*stdev(dsm) ; std over all parts
      partsm=where(maxsm gt thrsm,nprtsm) ; parts over threshold
      if nprtsm gt 0 then begin
        nsam2=2L^ceil(alog10(nsam)/alog10(2))
        addsam1=(nsam2-nsam)/2L & addsam2=nsam2-nsam-addsam1
        if addsam1 eq 0 then $
            if addsam2 eq 0 then ffti=fft(reform(intwfi),-1) $
                            else ffti=fft([reform(intwfi),fltarr(addsam2)],-1) $
                        else $
            if addsam2 eq 0 then ffti=fft([fltarr(addsam1),reform(intwfi)],-1) $
                            else ffti=fft([fltarr(addsam1),reform(intwfi),fltarr(addsam2)],-1)
        for ifex=0,nfex-1 do begin ;exclude frequencies from int waveform
         ffti[fexind[0,ifex]:fexind[1,ifex]]=complex(0.,0.)
         ffti[nsam-fexind[1,ifex]:nsam-fexind[0,ifex]]=complex(0.,0.)
        endfor
        intwfi_=(float(fft(ffti,1)))[addsam1:nsam+addsam1-1]
      ;
        plot_init,psfile,yfrac=1
        !p.charsize=1. & !P.MULTI = [0, 2, 2, 0, 0] &  !p.thick=2 & !x.thick=3 & !y.thick=3
        for iprtsm =0, (sm_maxfig < (nprtsm-1) ) do begin
          plot_wf,reform(wfi),times+' UT',fs=fs,tunit=1e-3, partsm[iprtsm]*pal, (partsm[iprtsm]+1)*pal-1,$
               ytit='Part '+strtrim(string(iprtsm+1),2)+'/'+$
               strtrim(string(nprtsm),2)+' |wf-'+strtrim(string(round(offss)),2)+'|>'+$
               strtrim(string(round(thrsm)),2)+' units '+$
               '!c !c  Raw data ['+strtrim(string(round(cali/1d3)),2)+' !mm!3T/s]',$
               tit=stationcode+ strtrim(string(i+1),1)+' (loop azimuth '+strtrim(string(round(azimi)),2)+'!uo!n)  '+$
                    strtrim(string(long(nsam/fs*1e3)),2)+' ms after '+waitt+' '+trigs[timeout]
          plot_wf,reform(intwfi_),times+' UT',fs=fs,tunit=1e-3, partsm[iprtsm]*pal, (partsm[iprtsm]+1)*pal-1,$
               ytit=sexcl+'!c!c Integrated data [nT]',tit=file+'   nGPS='+ngps
        end
        if nprtsm ge sm_maxfig then  $
        print,'------------------ more than '+strtrim(string(sm_maxfig),2)+' '+$
         strtrim(string(pal/fs*1e6,form='(f20.1)'),2)+'-us parts to plot ----------------------'
         ;
         plot_close
      endif ; nprtsm gt 0
    endif;  keyword_set(sm_plot)
end




pro textoverview,file,stationcode,wf,intwf,cal,times,waitt,terr,timeout
       ;
    ; text overview
    minraw=[min(wf[*,0]),min(wf[*,1])]
    maxraw=[max(wf[*,0]),max(wf[*,1])]
    sigraw=[stddev(wf[*,0],/double),stddev(wf[*,1],/double)]
    minint=[min(intwf[*,0]),min(intwf[*,1])]
    maxint=[max(intwf[*,0]),max(intwf[*,1])]
    sigint=[stddev(intwf[*,0],/double),stddev(intwf[*,1],/double)]
    txt='overview.txt'
    if file_test(txt) then openu,1,txt,/append else begin
      openw,1,txt
      printf,1,'File            ST UT of the first sample         TR Deadt'+$
               '  Min1  Min2  Max1  Max2  Sigma1  Sigma2 Mn1[nT] Mn2[nT]'+$
               ' Mx1[nT] Mx2[nT] Sg1[nT] Sg2[nT] Cal1 Cal2'
    end
    pos0=(strpos(file,'/',/reverse_s)>strpos(file,'\',/reverse_s))+1
    printf,1,strmid(file,pos0,15)+' '+stationcode+' '+times+$
      string(timeout,waitt,minraw,maxraw,sigraw,minint,maxint,sigint,form='(i2,5i6,8f8.1)')+$
      string(round(cal/1d3),form='(2i5)')
    close,1
end


