; history
;
;  2012-06-19 written
;  2012-06-22 histogram reports replaced by pro histo, pro debug added
;  2012-06-26 reports for separate events changed
;  2012-06-28 bug in reports for separate events fixed
;  2012-06-29 debug,/event and reports for separate events improved,
;  2012-07-15 bug in survey: wrong reporting max and min time and data in a block
;  2012-07-16 bug in bin->asc command coversion, wrongly reported EW and EB offsets (if >0 and SLoffs<0);
;             fixed bug in time strings (tns_str); timing of extremes id now reported
;  2012-07-26 SL, EW and EB offsets added to time (change of TM description for EW and EB);
;             bugs in timing of Eb and EW extremes and event indices fixed
;  2012-07-27 orb no. reported as integer
;  2012-07-29 unix/windows version switch included,
;  2012-08-30 renamed to blesk, problems for empty events fixed, altitude removed
;  2012-09-06 parameters of pro blesk reorganized, web figures sl.ps and ew.ps added
;             problem with a 1-s delayed log file name fixed
;  2012-09-08 asci file ew.txt and web figure ewspg.ps (4-day history spectrogram) added
;  2012-09-17 setting of the output ps and spectrogram range modified
;  2012-09-20 EW spectrogram modified
;  2012-09-21 bug in EW reading procedure for zero data in an event fixed
;  2012-09-29 bug in time axis of ewspg.ps fixed, b1 offset increased, output to blesklog.txt file added
;  2012-10-01 output to blesklog.txt replaced by output to blesk.txt
;  2012-10-06 bug in insertion of data gaps fixed in ewspg.ps
;  2012-11-08 processing of the dalis files added (input 'dal')
;  2012-11-15 description of analog gain parameter modified in write_conf_asc
;  2012-11-19 bug in 'dal' processing fixed
;  2012-12-14 'mex' processing added
;  2013-01-18 *.bin + *.sim extension handling added to con_conf
;  2013-01-25 gain info updated, characteristics of separate SR packets added
;  2013-01-29 bug in 'dal' batch processing fixed
;  2013-02-08 gain info updated, characteristics of separate SL blocks added, SL rank added
;  2013-02-14 SL rank changed to floating point
;  2013-02-20 board temperature added to ascii log files
;  2013-04-10 anomalies and init packet changed
;  2013-04-18 endtime in read RW corrected
;  2013-04-30 ascii output for Martin added
;  2013-05-07 ms replaced by pages in the *.asc control file,
;             high-pass filtering of integrated waveforms,
;             filename removed from blesk.txt
;  2013-05-08 the 1st order IIR filter in high-pass filtering of integrated waveforms
;             replaced by an inverse boxcar average filter
;  2013-05-15 radio EW events added;
;  2013-06-03 radio EW events modified, flag moved from the Taranis header to the event header
;  2013-06-10 radio end time bug fixed
;  2013-06-16 b1s added, oveview EW plots: added start time
;  2013-06-18 b1s calibration updated
;  2013-06-19 b1s calibration updated, packet header modified: Tf added
;  2013-06-20 histogram reporting improved by sfunct (header)
;  2013-08-10 fs dependent spectra (not for 80 MHz)
;  2013-08-12 -3dB on log spectra -> now 0dB corresponds to a gaussian noise with Veff = 1 quantization unit
;             bug in plotting long waveforms fixed (x axis to double)
;  2013-08-14 plot_wf modified to combine plots
;  2013-08-17 dB title in the EW log spectrum modified
;  2013-08-20-25 bug in "Configuration check: number of HFPE/RW intervals" fixed
;                bug in absolute time (usecs right shifted by 4 bits) fixed
;                in read_bin and fill_time, OLD *.dat need to be reprocessed
;  2013-08-31 switch overplotcolor added to plot_wf
;  2013-09-03 switch overplotsym added to plot_wf
;  2013-09-07 selected 2us figures added to plot_ew, averages and x-axis marks added to spectrograms
;  2013-09-12 bug in plot spect (ovl float->double) fixed
;  2013-09-14 averaging interval for overview ->40 min
;  2013-09-16 10MHz counters in event headers linked to the event data (each packet separately):
;             EM4 design change implemented; evtime corrected for processing delay
;  2013-09-18 XO relative timing implemented, real XO calibrated fs implemented,
;             consistency check for time differences between survey blocks added, data consistency error reporting improved
;  2013-09-28 uncorrected time reporting added, bugs in XO relative time fixed
;  2013-10-23 bug in XO relative time error reporting for EB data fixed
;  2013-11-09 --- OK --- in the txt log for init a final pakets,
;             XO relatve time for min max SL a EW, derivative min max added,
;             new header bits for memory and configuration errors implemented
;  2013-11-22 data type bugs in time corrections fixed, 'dau', 'em4' added
;  2014-03-20 hex command output in one single line, updated anomalies and gain steps, SA delay correction, ascii command file
;  2014-04-14 XO time calculated also for time before the last second
;  2014-05-19 first ini packet from a  multi-packet 1224 file processed. error keyword added for batch processing
;  2014-05-29 header analysis added for 1230-1 (radio) apids
;  2014-06-16 EMC options (3bits) added to configuration
;  2014-06-19 description of EMC options modified in asc conf file,  offset subtracted from radio I-Q data
;  2015-03-19 added averages in the SA and SB txt log files, updated tsens and tfpga
;  2015-03-28 dummy XO correction defined from init packet (1224) even if no survey data are received
;  2015-03-31 XO correction handling improved (done only from SA and SR, otherwise not done),
;             sequence of filters reversed in outputs according to the data where the filters are ordered 12...1
;  2015-06-06 bug in XO initialization fixed (initialized only if doesn't exist from before)
;  2015-09-12 b2s source (ersa station) added
;  2015-10-08 RW/PE test pattern options and new sampling frequencies implemented, SDRAM CRC error reported (RW,PE,EW,EB)
;  2016-06-29 10 Mhz counter in the event heder in read_Event changed from 32bit to 31bit
;  2016-11-12 RW and PE comments section in ASC command files updated to new bandwidths (5,10,20,80 kHz)
;  2017-01-11 Med: Ave: Std: added to EW txt output
;  2017-01-13 SL: removed '::' from txt output
;  2017-01-15 SA: added channel max&min for every block
;  2017-02-07 unknown APID bug fixed, default params corrected to EW buff 113, MaxEB 31
;  2017-02-08  RW gain bug fixed to mask the error bits
;  2017-02-10 SA block temperature writeout bug fixed
;  2017-03-15 a bug in smoothing integrated waveform before plotting it fixed
;  2017-03-30 31 bits 10 MHz counter in event descriptors implemented
;  2017-03-31 event data structure related to the firmware version, header checks added,
;             reordering packets in the sequence of packet number added
;  2017-04-01 bug fix: using 10MHz counter when it overflows just after the 1Hz pulse, end of sequence;
;             test pattern checks added, delays adjusted
;  2017-04-03 delay times adjusted: 175 ns for EW, 275 ns for SL, 75 ns for SR
;             (the same for test pattern and data where other delays compensate each other)
;  2017-04-05 preamp gain corrected in 1224 files; Conversion from HEX to asc and bin commands added
;  2017-04-06  bug in conversion from HEX to asc and bin commands fixed
;  2017-04-11  bug in design version recogtition fixed, rel. time correction report added
;  2017-04-13 bug in checking time in the oldest design fixed
;  2017-04-17 reltime cycle now allows more than 2^15 packets (error occurs only for EM with 2x event memory
;             or for more than 7 hours of survey data), output lists of packets limited to packet numbers
;  2017-04-26 keywords for plotting small intervals of EW added
;  2017-05-17 string formating compatible with idl version 8.6 and greater
;  2017-07-24 bug fix: reordering of event packets with only one gap
;  2021-02-25 spectrograms: source and min_fr_Hz added
;  2021-02-27 spectrograms: scaletit and tunits added

function ver
  return,'2021-02-27'
end


;---------------------------------------------------------------
;------------           measurement cycle         --------------
;---------------------------------------------------------------


pro blesk, input, source, filepath, error=error,$
           thrcsm=thrcsm,maxnsm=maxnsm,tintsm=tintsm

;
; IME-HF / BLESKA processing cycle
; EXAMPLE: blesk,'mex'
; parameters
;   input  'usb'... direct usb connection to BLESKA,
; 					a single data aquisition cycle is done.
;					Needs server.exe running, with server.cfg,
;					fw.ihx, and IME-HF*.rbf config files;
; 					Needs client10.exe and its config file jakub.txt
;                   to be present in the active directory  as well as the
;                   command file cmd.asc
;                   data occur in the subdirectory 'data'
;          'bin'... batch processing (from the active directory or from a '/' terminated
;                   path) of  all *.bin files w/o UT information and the correspondng *.log
;                   files from IME-HF / BLESKA
;          'dat'... batch processing of  *.dat (with UT information) files
;          'mex'... batch processing of  *.dat (with UT information) files, with mexic naming convention
;          'dal'... batch processing (from the active directory or from a '/' terminated
;                   path) of  all *.bin files w/o UT information (added from file name)
;                   from IME-HF / DALIS
;          'dau'... direct USB connection to dalis  processing
;                   ; 					a single data aquisition cycle is done.
; 					Needs command file cmd.asc to be present in the active directory
;
;   source 'b1m'... BLESKA-1 receiver (rustrel) with a 1-m magnetic loop antenna
;          'b1s'... BLESKA-1 receiver (rustrel) , gain 3 (1.22x) with the SLAVIA sensor, gain 4 (90x)
;          'b2s'... BLESKA-2 receiver (ersa), gain 1 (1x) with the SLAVIA sensor, gain 2
;          'em3'... TARANIS IME-HF <=EM3 receiver
;          'em4' of undefined... TARANIS IME-HF >=EM4 receiver (including FM)
;
;   filepath        a '/' terminated path for *.bin and *.log, or *.dat files
;   keywords:
;           error  output keyword: set of any error in data structure is found
;           tintsm, default=2e-6 - duration of small intervals in seconds
;           thrcsm, default=10   - threshold for max abs dev in units of std
;           maxnsm, default=450  - maximum number of smallest waveforms to be plotted
;----------------------- HW parameters ------------------
error=0
if not keyword_set(filepath) then filepath=''
if not keyword_set(input) then input='usb'
if not keyword_set(source) then source='fm'
case source of
  'b1m': begin
    off=4.9 ; correction of digitization offset in A/D units
    cal=1.28e6   ; offset corrected A/D units -> nT/s;   nT/s= cal*(A/D +offs)
    int=1
    tit='Integrated waveform (nT)'
    SL_thres=40 ; abs(waveform-offset)threshold for generating a new sl.ps figure
                ; noise is ~5 , should be above it
    fs=8e7      ; sampling frequency in Hz
    fhp=2e3     ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
   end
  'b1s': begin
    off=3.8 ; correction of digitization offset in A/D units
    cal=3.99e4   ; offset corrected A/D units -> nT/s;   nT/s= cal*(A/D +offs)
    int=1
    tit='Integrated waveform (nT)'
    SL_thres=600 ; abs(waveform-offset)threshold for generating a new sl.ps figure
                ; noise is ~10, monte carlo is ~100 , should be above it
    fs=8e7      ; nominal sampling frequency in Hz
    fhp=2e3     ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
   end
  'b2s': begin
    off= -5. ; correction of digitization offset in A/D units
    cal=2.975e4   ; offset corrected A/D units -> nT/s;   nT/s= cal*(A/D +offs)
    int=1
    tit='Integrated waveform (nT)'
    SL_thres=600 ; abs(waveform-offset)threshold for generating a new sl.ps figure
                ; noise is ~10, monte carlo is ~100 , should be above it
    fs=8e7      ; nominal sampling frequency in Hz
    fhp=2e3     ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
   end
  'em3': begin
    off=0 ; correction of digitization offset in A/D units
    cal=1   ; offset corrected A/D units -> result;   result= cal*(A/D +offs)
    int=1
    tit='Integrated waveform'
    SL_thres=40 ; abs(waveform-offset)threshold for generating a new sl.ps figure
                ; noise is ~5 , should be above it
    fs=8e7      ; nominal sampling frequency in Hz
    fhp=2e3     ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
   end
  'em4': begin
    off=0 ; correction of digitization offset in A/D units
    cal=1   ; offset corrected A/D units -> result;   result= cal*(A/D +offs)
    int=1
    tit='Integrated waveform'
    SL_thres=40 ; abs(waveform-offset)threshold for generating a new sl.ps figure
                ; noise is ~5 , should be above it
    fs=8e7      ; nominal sampling frequency in Hz
    fhp=2e3     ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
   end
  'fm' : begin
    off=0 ; correction of digitization offset in A/D units
    cal=1   ; offset corrected A/D units -> result;   result= cal*(A/D +offs)
    int=1
    tit='Integrated waveform'
    SL_thres=40 ; abs(waveform-offset)threshold for generating a new sl.ps figure
                ; noise is ~5 , should be above it
    fs=8e7      ; nominal sampling frequency in Hz
    fhp=2e3     ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
   end
end
;--------------------------------------------------------
;
  print,'-----------------> Blesk ver.'+ver()+'. Source='+source+', input='+input
  case input of
  'usb': begin
      conv_conf,'cmd.asc'
      if !VERSION.OS_FAMILY eq 'Windows' then $
        spawn,'client10 <jakub.txt >out.txt' else $
        spawn,'wine client10 147.231.75.214 <jakub.txt >out.txt'
    read_log,'out.txt', t0
    f=file_search('*.bin',count=nf)
    if nf gt 0 then begin
      logfile=strmid(f[0],0,strpos(f[0],'.',/REVERSE_SEARCH)-5)+'.log'
      file_copy,'out.txt','data/'+logfile
      for i=0,nf-1 do begin
        fill_time,f[i],t0
        file_move,strmid(f[i],0,strpos(f[i],'.',/REVERSE_SEARCH))+'.dat','data'
        file_delete,f[i]
      end
     end
    time=strmid(f[0],strpos(f[0],'.',/REVERSE_SEARCH)-20,15)
    path='data/'+time+'*.dat'
   end
  'dau': begin
      conv_conf,'cmd.asc'
      tns_str,(systime(/julian) -2400000.5d0) *86400000000000ull ,snow
      dalfile=strmid(snow[0],5,2)+'_'+strmid(snow[0],8,2)+'_'+strmid(snow[0],2,2)+'_'+$
              strmid(snow[0],11,2)+'_'+strmid(snow[0],14,2)+'_'+strmid(snow[0],17,2)
      spawn,'ifsim -a cmd.blc -o '+dalfile+'.bin'
      ;spawn,'ifsim -a cmd.blc
      if not file_test(dalfile+'.bin',/zero) then parse,dalfile+'.bin',dal=dal
      path=dal+'*.dal'
   end
  'dat': begin
      path=filepath+'*.dat'
    end
  'mex': begin
      path=filepath+'*.DAT'
    end
  'bin': begin
      path=filepath+'*.bin'
    end
  'dal': begin
      filedir=file_dirname(filepath)
      filebase=file_basename(filepath)
      f=file_search(filepath+'*.bin',count=nf)
      if nf gt 0 then  f=f[sort(f)]
      for  i=0,nf-1 do if not file_test(f[i],/zero) then parse,f[i],dal=dal
      if filebase eq '' then if filedir eq '.' then path='*.dal' $
                                               else path=filedir+path_sep()+'*.dal'$
                        else if filedir eq '.' then path=dal+'*.dal'  $
                                               else path=filedir+path_sep()+dal+'*.dal'
    end
  end
  ;
  curtime='' & inftime=''
  f=file_search(path,count=nf) & if nf gt 0 then f=f[sort(f)]
  for i=0,nf-1 do if not file_test(f[i],/zero) then begin
    case input of
      'mex': begin
               apid=strmid(f[i],strpos(f[i],'.',/REVERSE_SEARCH)-43,4)
               time=strmid(f[i],strpos(f[i],'.',/REVERSE_SEARCH)-31,15)
             end
      else: begin
             apid=strmid(f[i],strpos(f[i],'.',/REVERSE_SEARCH)-4,4)
             time=strmid(f[i],strpos(f[i],'.',/REVERSE_SEARCH)-20,15)
           end
     endcase
    ;
    if (inftime ne '')and((time ne curtime)or(i eq nf-1)) then begin
        openw,1,'blesk.txt'
        if n_elements(infew) eq 0 then infew='       '
        if n_elements(infsl) eq 0 then infsl='       '
        if n_elements(infsr) eq 0 then infsr='                     '
        printf,1,  inftime+ infew+infsl+infsr
        close,1
        inftime=''& infsr='                     '& infsl='       '& infew='       '
    end
    if time ne curtime then begin
     case input of
      'dal': begin ; DALIS file
	     print,'time: ',time
         t0=strmid(time,0,4)+'-'+strmid(time,4,2)+'-'+strmid(time,6,5)+$
                 ':'+strmid(time,11,2)+':'+strmid(time,13,2)+'.000'
         infos='DALIS file'
        end
      'dau': begin ; DALIS file
         t0=strmid(time,0,4)+'-'+strmid(time,4,2)+'-'+strmid(time,6,5)+$
                 ':'+strmid(time,11,2)+':'+strmid(time,13,2)+'.000'
         infos='DALIS file'
        end
      'mex': begin ; MEXIC file
         t0=strmid(time,0,4)+'-'+strmid(time,4,2)+'-'+strmid(time,6,5)+$
                 ':'+strmid(time,11,2)+':'+strmid(time,13,2)+'.000'
         infos='MEXIC file'
        end
       else: begin ; serching logfile for the same time or for a time by 1s sooner
        t0=0ul & infos='- !No log file found! -'
        logfile=strmid(f[i],0,strpos(f[i],'.',/REVERSE_SEARCH)-5)+'.log'
        if file_test(logfile,/read) and not file_test(logfile,/zero) $
         then read_log,logfile, t0,infos $
         else begin
          str_tns,strmid(time,0,4)+'-'+strmid(time,4,2)+'-'+strmid(time,6,5)+$
                  ':'+strmid(time,11,2)+':'+strmid(time,13,2)+'.000.000.000',tns
          tns=ulong64(tns-1000000000ll)
          tns_str,tns,ts
          logfile1=strmid(f[i],0,strpos(f[i],'.',/REVERSE_SEARCH)-20)+$
                  strmid(ts,0,4)+strmid(ts,5,2)+strmid(ts,8,2)+'_'+strmid(ts,11,2)+$
                  strmid(ts,14,2)+strmid(ts,17,2)+'.log'
          if file_test(logfile1,/read) and not file_test(logfile1,/zero) $
            then  read_log,logfile1, t0,infos $
            else print,'-----!! No logfile found: '+logfile+', '+logfile1+' ------'
         endelse
       end
      endcase
      curtime=time
    end
   if ((input eq 'bin') or (input eq 'dal')or (input eq 'dau')) then if (t0 ne 0) then begin
      fill_time,f[i],t0
      file_delete,f[i]
      f[i]=strmid(f[i],0,strpos(f[i],'.',/REVERSE_SEARCH))+'.dat'
   end
   if keyword_set(t0) then begin
    if n_elements(XO) eq 0 then XO={n:0} ; default, no data
    case apid of
    '1224': begin
              read_ini,f[i], conf,rwconf,nrwconf,a, err
              write_rep,f[i],a,err,error
            end
    '1225': begin
              read_Survey,f[i], XO, 1225, d,nd, dind,dtimes,dtimer,dfs,dgains, a,err, conf
              write_rep,f[i],a,err,error
              plot_SA, f[i], nd,d,dind,dtimes,dtimer,dgains,infos,err
            end
    '1226': begin
             read_Survey,f[i], XO, 1226, d,nd, dind,dtimes,dtimer,dfs,dgains, a,err, conf
              write_rep,f[i],a,err,error
              plot_SR, f[i], nd,d,dind,dtimes,dtimer,dgains,infos,err,$
                       int=int,cal=cal,off=off,tit=tit, inftime=inftime,infmax=infsr
            end
    '1227': begin
             read_Survey,f[i], XO, 1227, d,nd, dind,dtimes,dtimer,dfs,dgains, a,err, conf
              write_rep,f[i],a,err,error
              plot_SL, f[i], nd,d,dind,dtimes,dtimer,dfs,dgains,infos,err,$
                       int=int,cal=cal,tit=tit,off=off,SL_thres=SL_thres,infmax=infsl
            end
    '1228': begin
             read_Event,f[i], XO, 1228, ev,nev, evind, evfs,evtimer,evtimes,evtypes,evgains, a, err, conf
             write_rep,f[i],a,err,error
             plot_EB,f[i], nev,ev,evind,evtimes,evtimer,evtypes,evgains,infos,err,fs=evfs
            end
    '1229': begin
             read_Event,f[i], XO, 1229, ev,nev, evind, evfs,evtimer,evtimes,evtypes,evgains, a, err, conf
             write_rep,f[i],a,err,error
             plot_EW,f[i], nev,ev,evind,evtimes,evtimer,evtypes,evgains,infos,err,$
                            int=int,cal=cal,off=off,tit=tit,infmax=infew,fs=evfs,fhp=fhp,$
                            thrcsm=thrcsm,maxnsm=maxnsm,tintsm=tintsm
            end
    '1230': begin
              read_rw,1230,f[i],rwi,rwq, nrw,rwind, nrwi,rwtnss,rwtnse, rwbw, rwrate, $
                  amp,fre, rws, rwgain,rwagain,rwdgain,$
                  rwconf,nrwconf,a,err
              write_rep,f[i],a,err,error
              wav_RW, f[i], rwi,rwq, nrw,rwind, nrwi,rwtnss, rwrate
              plot_RW,f[i],rwi,rwq, nrw,rwind, nrwi,rwtnss,rwtnse, rwbw, rwrate, $
                  amp,fre, rws, rwgain,rwagain,rwdgain,$
                  rwconf,nrwconf
            end
    '1231': begin
              read_rw,1231,f[i],rwi,rwq, nrw,rwind, nrwi,rwtnss,rwtnse, rwbw, rwrate, $
                  amp,fre, rws, rwgain,rwagain,rwdgain,$
                  rwconf,nrwconf,a,err
              write_rep,f[i],a,err,error
              wav_RW, f[i], rwi,rwq, nrw,rwind, nrwi,rwtnss, rwrate
              plot_RW,f[i],rwi,rwq, nrw,rwind, nrwi,rwtnss,rwtnse, rwbw, rwrate, $
                  amp,fre, rws, rwgain,rwagain,rwdgain,$
                  rwconf,nrwconf
            end
    '1232': begin
              read_fin,f[i],  conf,rwconf,nrwconf,a,err
              write_rep,f[i],a,err,error
            end
     else: print,'-------- UNKNOWN APID:'+apid+' --------'
    end
   end else error=1
  end
  if nf eq 0 then begin
    print,'NO DATA FOUND <-----------------'
    error=1
  end
end


;------- parsing DALIS file

pro parse,file,dal=dal
 ;date time and apid coding: 11_07_12_15_43_35.bin -> 20120830_033103_1229.dal
 dot=strpos(file,'.',/REVERSE_SEARCH)
 dal='20'+strmid(file,dot-11,2)+strmid(file,dot-17,2)+strmid(file,dot-14,2)+$
          strmid(file,dot-9,3)+strmid(file,dot-5,2)+strmid(file,dot-2,2)+'_'
 ;dal=strmid(file,dot-15,15)+'_'
 print,'--- Parsing DALIS file '+file
 read_bin,file,npacks,packs,tns
 print,'---   '+string(npacks)+' packets.'
 h=histogram(packs.id,omin=m,rever=r)
 for i=0,n_elements(h)-1 do if h[i] gt 0 then begin
   dalfile=dal+string(m+i,form='(i4.4)')+'.dal'
   openw,1,dalfile,/SWAP_IF_LITTLE_ENDIAN
   writeu,1,packs[r[r[i]:r[i+1]-1]]
   close,1
   print,'      ->'+dalfile+string(h[i])+' packets.'
 end
 print,'--- Parsing done.
end



;-----------------------cal -------------------------

function calib_waveform, wf, cal=cal
;   Keyword cal
;     for a magnetic loop, radius R:
;         dB/dt [T/s] =U[V]/(pi*R[m]^2)= 1.97V/1966/(pi*R[m]^2)*Value
;         -> cal = 3.19e-4/R[m]^2  (default R=0.5 ->  cal=1.28e-3)
;     for an electric dipole, length L between geometrical
;     midpoints of the two monopoles:
;         E[V/m] = U[V]/L[m] = 1.97V/1966/L[m]*Value
;         -> cal = 1.00e-3/L[m]
;
;   calibration
    if not keyword_set(cal) then cal=1; 1.28e-3
    return, wf*cal
end



;---------------------------------------------------------------
;------------    Survey filterbank data SA    ------------------
;------------   Survey  regular waveforms SR  ------------------
;------------   Survey selected waveforms SL  ------------------
;---------------------------------------------------------------




pro plot_SA, file, nd,d,dind,dtimes,dtimer,dgains,infos,err
  ;
  ; multipage ps overview of all SA events in the file
  ; output goes to a ps file with the same name
  ; y axis autoscaled
  ;
  if nd gt 0 then begin
   plot_init,strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.ps'
   for j =0,nd-1 do begin
     np=dind[j+1]-dind[j]
     x=findgen(np)*0.00738
     ymax=1>max(d[*,dind[j]:dind[j+1]-1],min=ymin)
     if err then errs='(!!)' else errs=''
     plot,[0,1],/nodata,yran=[ymin-(ymax-ymin)/15. ,ymax],xran=[x[0],x[np-1]],xst=1,yst=19,$
               xtit='seconds from '+dtimes(j)+$
                ' (XO:'+strtrim(string(dtimer[j],form='(f30.9)'),2)+')',$
               ytit='SA block '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nd),2)+$
               '!c !c Filter bank data (TMU)',$
               tit='!d'+infos+$
               ' !dFile '+file+errs+' Gain='+dgains[j]+'!n',color=255
      for k =0,11 do oplot,x,d[k,dind[j]:dind[j+1]-1],color=253-20*(11-k),thick=3
      ytxt=ymin-(ymax-ymin)/16.
      xyouts,x[np/80],ytxt,'Channels:',color=255,charsize=0.7
      for k =0,11 do  xyouts,x[np/15*(k+2)],ytxt,string(k+1),color=253-20*k,charsize=0.7
   end
   plot_close
  end
end


pro sr_arrows, apfs
    for i=0,17 do begin
      if i gt 0 then oplot,1360d0*[i,i]/apfs,!y.crange,color=254,line=0
      arrow,color=254,1360d0*(i+1)/apfs,(!y.crange[0]*29+!y.crange[1])/30.,$
                      1360d0*i/apfs,(!y.crange[0]*59+!y.crange[1])/60.,/data,hsize=-0.4,/solid
    end
end


pro plot_SR, file, nd,d,dind,dtimes,dtimer, dgains,infos,err,$
            der=der, int=int,cal=cal,off=off,tit=tit,inftime=inftime,infmax=infmax
  ;
  ; multipage ps overview of all SR packets in the file
  ; output goes to a ps file with the same name
  ; y axis autoscaled
  ; keywords: ytit - marks y axis by it
  ;           int - integrates waveform, 'int' is added to
  ;                 the file name
  ;           der - derivative of the waveform, 'der' is added to
  ;                 the file name
  ;           cal - calibration coefficient, see function calib_waveform
  ;           off - offset in TMU
  ;
  if nd gt 0 then begin
   plot_init,strmid(file,0,strlen(file)-4)+'.ps'
   if err then errs='(!!)' else errs=''
   if not keyword_set(off) then off=0
   if not keyword_set(cal) then cal=1
   if not keyword_set(tit) then tit='Waveform (TMU)'
   ;
    apfs=1360d0*18/15.0109200d0
   dmax=0 & jmax=-1
   for j =0,nd/18-1 do begin
     plot_wf, d(dind(j*18):dind((j+1)*18)-1), dtimes(j*18)+$
                ' (XO:'+strtrim(string(dtimer[j*18],form='(f30.9)'),2)+')',$
               tunits=1,fs=apfs,fhp=0,ytit='SR block '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nd/18),2)+$
               '!c !c Waveform (TMU)',noisescale=20,$
               tit='!d'+infos+$
               ' !dFile '+file+errs+' Gain='+dgains[j]+'!n'
     sr_arrows,apfs
     plot_spect,d(dind(j*18):dind((j+1)*18)-1)+off,$
               ytit='Power spectrum (dB)'
     plot_wf,calib_waveform(d(dind(j*18):dind((j+1)*18)-1)+off,cal=cal*apfs/8e7),dtimes(j*18),$
               der=der,int=int,tunits=1,fs=apfs,ytit=tit,noisescale=20*cal*apfs/8e7
     sr_arrows,apfs
     dmax=max([dmax,abs(d[dind(j*18):dind((j+1)*18)-1]+off)],imax)
     if (imax ne 0) then jmax=j
   end
   for j =0,nd-1 do begin
    plot_wf,d(dind(j):dind(j+1)-1),dtimes(j)+$
                ' (XO:'+strtrim(string(dtimer[j],form='(f30.9)'),2)+')',$
               tunits=1e-6,ytit='SR packet '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nd),2)+$
               '!c !c Waveform (TMU)',noisescale=20,$
               tit='!d'+infos+$
               ' !dFile '+file+errs+' Gain='+dgains[j]+'!n'
    plot_spect,d(dind(j):dind(j+1)-1)+off,$
               ytit='Power spectrum (dB)'
    plot_wf,calib_waveform( d [dind(j):dind(j+1)-1]+off,cal=cal), dtimes(j),$
               der=der,int=int,$
               tunits=1e-6,ytit=tit,noisescale=20*cal

   end
   plot_close
       ;
    ;txt file--------------------
    avewf=off+(moment(d,sdev=sdewf))(0)
 ;   inftime=strmid(file,0,strlen(file)-9)+' '+strmid(dtimes(0),11,8)+'-'+$
 ;             strmid(dtimes(nd-1),11,8)
   inftime=strmid(dtimes(0),11,8)+'-'+strmid(dtimes(nd-1),11,8)
    ;infmax=string(dmax,jmax,avewf,sdewf,form='(i5,i3,f7.2,f7.2)')
    infmax=string(dmax,avewf,sdewf,form='(i5,f7.2,f7.2)')
   ;
  end
end


pro plot_SL, file, nd,d,dind,dtimes,dtimer,dfs,dgains,infos,err,$
            der=der, int=int,cal=cal,off=off, tit=tit,SL_thres=SL_thres,infmax=infmax
  ;
  ; multipage ps overview of all SL blocks in the file
  ; output goes to a ps file with the same name
  ; y axis autoscaled
  ; keywords: ytit - marks y axis by it
  ;           int - integrates waveform, 'int' is added to
  ;                 the file name
  ;           der - derivative of the waveform, 'der' is added to
  ;                 the file name
  ;           cal - calibration coefficient, see function calib_waveform
  ;           off - offset in TMU
  ;
  if nd gt 0 then begin
   plot_init,strmid(file,0,strlen(file)-4)+'.ps'
   if err then errs='(!!)' else errs=''
   ; whole blocks first
   if not keyword_set(cal) then cal=1
   if not keyword_set(off) then off=0
   if not keyword_set(tit) then tit='Waveform (TMU)'
   if not keyword_set(SL_thres) then SL_thes=10
   ;
   ;
   dmax=0 & jmax=-1
   for j =0,nd-1 do begin
    plot_wf,d(dind(j):dind(j+1)-1),dtimes(j)+$
                ' (XO:'+strtrim(string(dtimer[j],form='(f30.9)'),2)+')',fs=dfs[j],$
               tunits=1e-6,ytit='SL block '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nd),2)+$
               '!c !c Waveform (TMU)',noisescale=20,$
               tit='!d'+infos+$
               ' !dFile '+file+errs+' Gain='+dgains[j]+' BLESK-'+ver()+'!n'
    plot_spect,d(dind(j):dind(j+1)-1)+off,fs=dfs[j],$
               ytit='Power spectrum (dB)'
    plot_wf,calib_waveform(d(dind(j):dind(j+1)-1)+off,cal=cal),fs=dfs[j],dtimes(j),$
               der=der,int=int,tunits=1e-6,ytit=tit,noisescale=20*cal
    dmax=max([dmax,abs(d[dind(j):dind(j+1)-1]+off)],imax)
    if (imax gt 0) then jmax=j
   end
  ; small parts
   pal=1360L ; parts of one packet = 1360 samples
   for j =0,nd-1 do $
    for i=0,(dind(j+1)-dind(j))/pal-1 do begin
       plot_wf,d(dind(j):dind(j+1)-1),dtimes(j)+$
                ' (XO:'+strtrim(string(dtimer[j],form='(f30.9)'),2)+')',fs=dfs[j],$
               i*pal,(i+1)*pal-1,$
               noisescale=20,tunits=1e-6,$
               ytit='SL block '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nd),2)+$
              '  Packet '+strtrim(string(i+1),2)+'/'+$
               strtrim(string((dind(j+1)-dind(j))/pal),2)+$
               '!c !c  Waveform (TMU)',$
               tit='!d'+infos+$
               ' !dFile '+file+errs+' Gain='+dgains[j]+' BLESK-'+ver()+'!n'
       plot_spect,d(dind(j):dind(j+1)-1)+off,fs=dfs[j],$
               i*pal,$
               min([(i+1)*pal-1,dind(j+1)-dind(j)-1]),$
               ytit='Power spectrum (dB)'
       plot_wf,calib_waveform(d(dind(j):dind(j+1)-1)+off,cal=cal),dtimes(j),fs=dfs[j],$
               i*pal,(i+1)*pal-1,$
               der=der,int=int,noisescale=20*cal,tunits=1e-6,$
               ytit=tit
    end
   plot_close
    ;
    ;txt file--------------------
    ;infmax=string(dmax,jmax,form='(i5,i3)')
    infmax=string(dmax,form='(i5)')
   ;
   ; sl ps file
   if (dmax gt SL_thres) then begin
    slpsfile=strmid(file,0,strlen(file)-4)+'_sl.ps'
    plot_init,slpsfile,yfrac=0.6
    !p.multi=0 &  !p.thick=2 & !x.thick=3 & !y.thick=3
    plot_wf,calib_waveform(d(dind(jmax):dind(jmax+1)-1)+off,cal=cal),$
               strmid(dtimes(jmax),0,27),$
               der=der,int=int,tunits=1e-6,ytit='SL block '+$
               strtrim(string(jmax+1),2)+'/'+strtrim(string(nd),2)+$
               '!c !c '+tit,noisescale=20*cal,$
               tit=strtrim(string((dind(jmax+1)-dind(jmax))/80l),2)+$
               ' !mm!3s from '+strmid(dtimes(jmax),0,27)+'!c!d'+infos+$
               ' !dFile '+file+errs+' Gain='+dgains[jmax]+$
               ' BLESK-'+ver()+'!n', pos=[0.1,0.13,0.99,0.88]
    plot_close
    file_copy,slpsfile,'sl.ps',/overwrite
   end
  end
end




pro read_Survey,file, XO, apid, d,nd, dind,dtimes,dtimer,dfs,dgains, a,err, confi
  ; common survey data reading routine
  ;
  ;
  ;
; extracts SA filterbank or SR or SL waveform data (output nd - number of continuous data intervals,
;                        output dind - array of nd starting
;                              indices of data intervals & last index
;                        output d -  SA filterbank data, array[12,dind(nd)-1])
;                                    SR,SL waveform data, array[dind(nd)-1])
;  dtime of the beginning of each continuous data interval
;   is stored in and array of strings
;   dgains  array of strings
;
  err=0 & nd=0
  read_bin,file,npacks,packs,tns
  ;
  case apid of
   1225: begin
     pgpk=18ul; packets per block
     dpk=18ul; packets per continuous interval
     pkv=1356ul; values per packet,
     dats='SA' ;
     dtdat=615ul*113ul*120ul;  time (100ns) between packets
    check_headers,file,npacks, packs,apid, tns,  t100nsc,  a, err,XO=XO; XO is only read from SA or SR
    end
   1226: begin
     pgpk=18ul; packets per block
     dpk=1ul; packets per continuous interval
     pkv=1360ul; values per packet,
     dats='SR' ;
     dtdat=8339400ul;  time (100ns) between packets
    check_headers,file,npacks, packs,apid, tns,  t100nsc,  a, err,XO=XO; XO is only read from SA or SR
    end
   1227: begin
     pgpk=18ul; packets per block
     dpk=18ul; packets per continuous interval
     pkv=1360ul; values per packet,
     dats='SL' ;
     dtdat=170ul;  time (100ns) between packets
    check_headers,file,npacks, packs,apid, tns,  t100nsc,  a, err
   end
  end
  ;
  id=where(packs.id eq apid, ndp)
  nbl=ndp/pgpk & expacks=npacks-nbl*pgpk
  ;
  a=[a,' ','SURVEY DATA READER','-- File '+file+', '+strtrim(string(nbl),2)+$
    ' '+dats+' data blocks found. Extra packets: '+strtrim(string(expacks),2)+'.']
  if expacks gt 0 then begin
    err=1; expacks should be 0
    a=[a,'   ERROR -  Extra packets']
  end
  ;
  if nbl gt 0 then begin
    id=id(0:nbl*pgpk-1)
    get_sud, packs[id], id, dtdat, tns[id], C10Mhz, rank,  blno, err,a, sl=(apid eq 1227)
    reltime, XO, packs[id].hzcnt,C10Mhz, rt,corr,a=a
    ;
    case apid of
    1225: begin
     get_filterbank,packs[id], nbl*pgpk, d
     dpfs=corr/(12d-6 * 615) ; fs for each packet: 615 12us samples into one SA record
     tns_str,tns[id]-long64(7381900ll/corr),dptime
     dptimer=rt-7.3819d-3/corr      ;7.3729ms  delay
     a=[a,'--','   MIU or XO time of the 1st SA sample = header (MIU or corrected XO) time - 7.3819ms delay']
     end
    1226: begin
     ; Experiment with time triggered EW and SR (IDEE tests): SR has the same data with exactly 9 samples (112.5 ns)
     ; less delay than EW ( => SR pattern delay would be 5 samples) .
     ; Pattern (when jakub sequence unusually started later):  SR pattern delay should be 6 samples (75 ns)
     ; With MIU SR 1st sample of the first packet has a delay 300 ns (24samp) after 10MHz counter =0 with jakub the delay is 100 ns (8samp)
     ; EW experiment: 1s gps pulses seen with jakub in EW (forked at the input do Virtex, as if they go from the output of
     ; fpga-xp .......delay of 180 ns .... approx the same as the pattern->EWdata delay, 175ns
     ; it means that the analog propagation delay + digital delay up to the pattern inserting point is compensated by the
     ; delay of 1Hz clock propagation inside Virtex
     ; definition:
     ;             SR delay 75  ns   6  samples  (= pattern-> SRdata)
     ;             EW delay 175 ns   14 samples  (= pattern->EWdata)
     ;             SL delay 275 ns   22 samples  (= pattern->SLdata)
     ;
     ; experiment with pattern: SL needs to have by 8 samples (by 100ns) more delay than EW
     ; SL pattern delay is 22 samples (275 ns)
     ;
     get_waveforms,packs[id], nbl*pgpk, d
     dpfs=corr*8d7 ; fs for each packet
     tns_str,tns[id]-long64(75d0/corr),dptime
     dptimer=rt-75d-9/corr
     a=[a,'--','   MIU or XO time of the 1st SR sample = header (MIU or corrected XO) time - 75 ns delay']
     end
    1227: begin
       get_waveforms,packs[id], nbl*pgpk, d
     dpfs=corr*8d7 ; fs for each packet
     tns_str,tns[id]+long64(( (0ll < long64(confi.Time_Off_SL*4000ll) )-275ll )/corr),dptime
     dptimer=rt+( (0d0 < confi.Time_Off_SL*4d-6) -275d-9 )/corr
     a=[a,'--','   MIU or XO time of the 1st SL sample = header (MIU or corrected XO) time +offset(if negative) - 275 ns delay']
     end
    end

    ;
    a=[a,'   Samplig frequency correction factor in all '+dats+' packets: '+strtrim(string(min(corr), form='(f30.9)'),2)+$
             ' - '+strtrim(string(max(corr), form='(f30.9)'),2)]
    a=[a,'   Min 1st sample MIU time among all '+dats+' packets: '+min(dptime)]
    a=[a,'   Max 1st sample MIU time among all '+dats+' packets: '+max(dptime)]
    a=[a,'   Min 1st sample XO time among all '+dats+' packets: '+strtrim(string(min(dptimer), form='(f30.9)'),2)]
    a=[a,'   Max 1st sample XO time among all '+dats+' packets: '+strtrim(string(max(dptimer), form='(f30.9)'),2)]
    a=[a,'   TM values from '+strtrim(string(min(d)),2)+$
      ' to '+strtrim(string(max(d)),2)+' median:'+strtrim(string(median(d)),2)]
    ;
    ;
    nd=nbl*pgpk/dpk   ; continuous data blocks
    dindp=lindgen(nd+1)*dpk  ; indices of packets starting cont. data blocks
    dtimes=dptime[dindp[0:nd-1]]
    dtimer=dptimer[dindp[0:nd-1]]
    dfs=dpfs[dindp[0:nd-1]]
    gains=['1x','4x','13x','_13x','1/2x(ant4)','4/2x(ant4)','13/2x(ant4)','_13/2x(ant4)','1/2x(ant3)','4/2x(ant3)','13/2x(ant3)','_13/2x(ant3)','1/15x','4/15x','13/15x','_13/15x']
    dgains=gains[packs[id[dindp[0:nd-1]]].gains_bw/16b]
    ;
  mode=(packs[id[0]].iinf mod 16) &  tst=((mode ge 10)and(mode le 13)) ;____if TEST pattern
 ;
  case apid of
   1225: begin
     d=reform(d,12,pgpk*nbl*pkv/12)    ;*0.83252; 512/615
     dind=dindp*pkv/12
     if tst then  a=[a,'--','TEST PATTERN MODE DETECTED -> test pattern analysis of blocks of SA data packets:'] $
            else  a=[a,'--','Characteristics in blocks of SA data packets:']
     for i =0L, nbl-1 do begin
        a=[a,'',' - Block number: '+strtrim(string(i+1),2)+' starting at '+dtimes(i)+$
                     ' (XO: '+ strtrim(string(dtimer[i], form='(f30.9)'),2)+$
                    ' - '+strtrim(string(dtimer[i]+(pkv*pgpk/12)/dfs[i], form='(f30.9)'),2)+$
                    ', t='+strtrim(string(packs[id[i*pgpk]].tb*1.608-273.15, form='(f30.1)'),2)+' degC)']
        str_tns,dtimes(i),t0
        if tst then begin ; test pattern mode
        for j = 0,10 do begin ;__check differences between channels
          devj=0+d[j+1,i*pkv*pgpk/12:(i+1)*pkv*pgpk/12-1]-d[j,i*pkv*pgpk/12:(i+1)*pkv*pgpk/12-1]
          iok=where( (devj eq 341), nok)
          perc=100.*nok/(pkv*pgpk/12)
          if perc lt 70. then begin
            a=[a,'   Channel '+string(12-j-1,form='(i2.2)')+' - Channel '+string(12-j,form='(i2.2)')+$
                 ': '+strtrim(string(perc, form='(f30.0)'),2)+'% of differences = 341. ERROR - less than 70%']
               err=1
          end else a=[a,'   Channel '+string(12-j-1,form='(i2.2)')+' - Channel '+string(12-j,form='(i2.2)')+$
                      ': '+strtrim(string(perc, form='(f30.0)'),2)+'% of differences = 341  - OK - More than 70%']
        end
        devj=0+d[0,i*pkv*pgpk/12+1:(i+1)*pkv*pgpk/12-1]-d[0,i*pkv*pgpk/12:(i+1)*pkv*pgpk/12-2] ;__check time sequence of channel 1
        iok=where( (devj eq 615), nok)
          perc=100.*nok/(pkv*pgpk/12)
          if perc lt 60. then begin
            a=[a,'   Channel *12*' +$
                 ': '+strtrim(string(perc, form='(f30.0)'),2)+'% of steps = 615. ERROR - less than 60%']
               err=1
          end else a=[a,'   Channel *12*'+$
                      ': '+strtrim(string(perc, form='(f30.0)'),2)+'% of steps = 615  - OK - More than 60%']
        end
         for j = 0,11 do begin
          tmmax=max(d[j,i*pkv*pgpk/12:(i+1)*pkv*pgpk/12-1])
          tmmin=min(d[j,i*pkv*pgpk/12:(i+1)*pkv*pgpk/12-1])
          tmmed=median(d[j,i*pkv*pgpk/12:(i+1)*pkv*pgpk/12-1])
          a=[a,'   Channel '+string(12-j,form='(i2.2)')+$
              ': TM values - Median:'+strtrim(string(tmmed),2)+$
              ' Max:'+string(tmmax, form='(i5)')+$
              ' Min:'+string(tmmin, form='(i5)')]
         end
      end
    end
   1226: begin
    d=reform(d,pkv*nbl*pgpk)
    dind=dindp*pkv
    if tst then  a=[a,'--','TEST PATTERN MODE DETECTED -> test pattern analysis of separate SR data packets:'] $
           else  a=[a,'--','Characteristics of separate SR data packets:']
    for i =0L, nbl*pgpk-1 do begin
      a=[a,'',' - Packet number: '+strtrim(string(i+1),2)+' starting at '+dptime[i]+$
                    ' (XO: '+ strtrim(string(dptimer[i], form='(f30.9)'),2)+$
                    ' - '+strtrim(string(dptimer[i]+pkv/dfs[i], form='(f30.9)'),2)+$
                    ', t='+strtrim(string( packs[id[i]].tb*1.608-273.15, form='(f30.1)'),2)+' degC)']
       ;if  (dptimer[i] gt 13.998)and (dptimer[i] lt 14.4) then begin &  sr= d[i*pkv:(i+1)*pkv-1] & tsr=dptimer[i] & save,tsr,sr & end
      str_tns,dptime[i],t0
      if tst then begin
        devj=d[i*pkv+1:(i+1)*pkv-1]-d[i*pkv:(i+1)*pkv-2] ;__check time sequence
        inok=where( (devj ne 1)and(devj ne -4095), nnok)
        if nnok eq 0 then a=[a,'   Monotonous time sequence, no irregularities - OK -'] $
        else begin
          a=[a,'   Found '+strtrim(string(nnok),2)+' time sequence irregularities:']
          trel=dblarr(nnok) & tval=intarr(nnok)
          for ii=0,nnok-1 do begin
            tns_str,t0+ulong64(1d9*inok[ii]/dpfs[i]),tmit
            trel[ii]=dptimer[i]+(inok[ii]+1)/dpfs[i]
            tval[ii]=d[i*pkv+1+inok[ii]]
            a=[a,'   Value of '+strtrim(string(tval[ii]),2)+$
                 ' succeeding '+strtrim(string(d[i*pkv+inok[ii]]),2)+$
                 ' at offset of '+strtrim(string(1d6*(inok[ii]+1)/dpfs[i],form='(f230.3)'),2)+' usec -> '+tmit+$
                 ' (XO: '+ strtrim(string(trel[ii], form='(f30.9)'),2)+')']
          end
          ADCcorr=0d-9
          if (nnok eq 2) then $
            if (abs(trel[0]-round(trel[0])+ADCcorr) lt 13d-9)and(abs(trel[0]-trel[1]) lt 13d-9)and(tval[0] eq 0)and(tval[1] eq 0) $
            then a=[a,'   Expected reset of the test sequence by two zeros at entire second - OK -'] $
            else begin
              err=1 & a=[a,'   ERROR: unexpected irregularities (sequence should be reset by two zeros at entire second).']
            end  $
          else begin
              err=1 & a=[a,'   ERROR: too many irregularities (sequence should be reset only at entire seconds).']
          end
        endelse
      end  ; __not  pattern: Extremes
        tmmax=max(d[i*pkv:(i+1)*pkv-1],min=tmmin)
        tmave=total(d[i*pkv:(i+1)*pkv-1])/float(pkv)
        tmstd=sqrt(total( (d[i*pkv:(i+1)*pkv-1]-tmave)^2 )/float(pkv-1) )
        a=[a,'   TM values: Max = '+string(tmmax, form='(i5)')+'  Min = '+string(tmmin, form='(i5)')+$
                         '  Ave = '+string(tmave, form='(f8.2)')+'  Std = '+string(tmstd, form='(f8.2)')]
     end
    end
   1227: begin
    d=reform(d,pkv*nbl*pgpk)
    dind=dindp*pkv
    if tst then  a=[a,'--','TEST PATTERN MODE DETECTED -> test pattern analysis of separate SL data blocks:'] $
           else  a=[a,'--','Extremes and waveform characteristics in blocks of SL data packets:']
    for i =0L, nbl-1 do begin
      if rank[dindp[i]] eq 0 then rerr=' --- ERROR: ZERO RANK ---' else rerr=''
      a=[a,'',' - Block number: '+strtrim(string(i+1),2)+' starting at '+dtimes(i)+$
                    ' (XO: '+ strtrim(string(dtimer[i], form='(f30.9)'),2)+$
                    ' - '+strtrim(string(dtimer[i]+(pkv*pgpk)/dfs[i], form='(f30.9)'),2)+$
                    ', t='+strtrim(string(packs[id[i*pgpk]].tb*1.608-273.15, form='(f30.1)'),2)+' degC)'+$
                    '  Rank: '+ string(rank[dindp[i]], form='(I5.5)') + rerr]
      str_tns,dtimes(i),t0
      if tst then begin ; test pattern mode
        devj=d[i*pkv*pgpk+1:(i+1)*pkv*pgpk-1]-d[i*pkv*pgpk:(i+1)*pkv*pgpk-2] ;__check time sequence
        inok=where( (devj ne 1)and(devj ne -4095), nnok)
        if nnok eq 0 then a=[a,'   Monotonous time sequence, no irregularities - OK -'] $
        else begin
          a=[a,'   Found '+strtrim(string(nnok),2)+' time sequence irregularities:']
          trel=dblarr(nnok) & tval=intarr(nnok)
          for ii=0,nnok-1 do begin
            tns_str,t0+ulong64(1d9*inok[ii]/dfs[i]),tmit
            trel[ii]=dtimer[i]+(inok[ii]+1)/dfs[i]
            tval[ii]=d[i*pkv*pgpk+1+inok[ii]]
            a=[a,'   Value of '+strtrim(string(tval[ii]),2)+$
                 ' succeeding '+strtrim(string(d[i*pkv*pgpk+inok[ii]]),2)+$
                 ' at offset of '+strtrim(string(1d6*(inok[ii]+1)/dfs[i],form='(f230.3)'),2)+' usec -> '+tmit+$
                 ' (XO: '+ strtrim(string(trel[ii], form='(f30.9)'),2)+')']
          end
          ADCcorr=0d-9
          if (nnok eq 2) then $
            if (abs(trel[0]-round(trel[0])+ADCcorr) lt 13d-9)and(abs(trel[0]-trel[1]) lt 13d-9)and(tval[0] eq 0)and(tval[1] eq 0) $
                 then a=[a,'   Expected reset of the test sequence by two zeros at entire second - OK -'] $
            else begin
              err=1 & a=[a,'   ERROR: unexpected irregularities (sequence should be reset by two zeros at entire second).']
            end  $
          else begin
              err=1 & a=[a,'   ERROR: too many irregularities (sequence should be reset only at entire seconds).']
          end
        endelse
      end  ; __not  pattern: Extremes
        tmmax=max(d[i*pkv*pgpk:(i+1)*pkv*pgpk-1],min=tmmin)
        tmave=total(d[i*pkv*pgpk:(i+1)*pkv*pgpk-1])/float(pkv*pgpk)
        tmstd=sqrt(total( (d[i*pkv*pgpk:(i+1)*pkv*pgpk-1]-tmave)^2 )/float(pkv*pgpk-1) )
        a=[a,'   TM values: Max = '+string(tmmax, form='(i5)')+'  Min = '+string(tmmin, form='(i5)')+$
                         '  Ave = '+string(tmave, form='(f8.2)')+'  Std = '+string(tmstd, form='(f8.2)')]
        tmv=min(d[i*pkv*pgpk:(i+1)*pkv*pgpk-1],tmi)
        tns_str,t0+ulong64(1d9*tmi/dfs[i]),tmit
        a=[a,'   Min TM value:'+strtrim(string(tmv),2)+' at an offset of '+$
             strtrim(string(1d6*tmi/dfs[i],form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(dtimer[i]+tmi/dfs[i], form='(f30.9)'),2)+')']
        tmv=max(d[i*pkv*pgpk:(i+1)*pkv*pgpk-1],tmi)
        tns_str,t0+ulong64(1d9*tmi/dfs[i]),tmit
        a=[a,'   Max TM value:'+strtrim(string(tmv),2)+' at an offset of '+$
             strtrim(string(1d6*tmi/dfs[i],form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(dtimer[i]+tmi/dfs[i], form='(f30.9)'),2)+')']
            ;
        savgolFilter = SAVGOL(4, 4, 1, 2)*corr[dindp[i]]*8d1 ; use 16,16,1,4 for more smoothing
        mx=max(min=mn, /nan,$
               CONVOL(d[i*pkv*pgpk:(i+1)*pkv*pgpk-1], savgolFilter, /EDGE_TRUNCATE,/nan),$
               tmx,sub=tmn)
        tns_str,t0+ulong64(1d9*tmx/dfs[i]),tmit
        a=[a,'   Max dTM/dt [1/us]:'+strtrim(string(mx, form='(e10.2)'),2)+' at an offset of '+$
             strtrim(string(1d6*tmx/dfs[i],form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(dtimer[i]+tmx/dfs[i], form='(f30.9)'),2)+')']
        tns_str,t0+ulong64(1d9*tmn/dfs[i]),tmit
        a=[a,'   Min dTM/dt [1/us]:'+strtrim(string(mn, form='(e10.2)'),2)+' at an offset of '+$
             strtrim(string(1d6*tmn/dfs[i],form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(dtimer[i]+tmn/dfs[i], form='(f30.9)'),2)+')']
      end
    end
  end
;
    errstr=['-- OK --','-- ERROR --']
    print,'   --> ', strtrim(string(nbl),2),' '+dats+' data blocks read.   '+errstr[err]
    ;
  end
 end



pro get_sud, packs, id, dtdat, tns,  C10Mhz, rank, blno, err,a,sl=sl
;
; extracts the survey data descriptor information
;   last 4 bytes of the 10 Mhz counter ... C10Mhz,
;   last 16 bits of rank for SL, 0 otherwise ...rank
;   block number...blno,
; from
;   packet records...packs
; reports anomalies in
;   a anomaly diagnistic string array...a
;   error flag...err
;
    C10Mhz=packs.dat(0)*16777216ul +packs.dat(1)*65536ul+packs.dat(2)*256ul+packs.dat(3)
    ranke=packs.dat(4)/32
    rank=2^ranke * ( (packs.dat(4)-ranke*32u)*256u+packs.dat(5))
    blno=packs.dat(6)*256u+packs.dat(7)
    np=n_elements(packs)
    ;
    a=[a,'--','Data structure of blocks of survey data packets:']
    for i =0L, np/18-1 do begin
     a=[a,'',' - Block number: '+strtrim(string(i+1),2)+' from the sequence of 18 packets.']
     a=[a,'   Block numbers in the data:']
       histo,blno[i*18:(i+1)*18-1], id[i*18:(i+1)*18-1], a, vald,fred
       if (n_elements(vald) ne 1) or (fred[0] ne 18) then begin
         err=1
         a=[a,'   !!! Inconsistent block numbers in the data.']
       end else a[n_elements(a)-1]=a[n_elements(a)-1]+'  -- OK --'
     ;
     a=[a,'   Ranks for SL, 0 for SR and SA:']
       histo,rank[i*18:(i+1)*18-1], id[i*18:(i+1)*18-1],  a, vald,fred
       if (n_elements(vald) ne 1) or (fred[0] ne 18) then begin
         err=1
         a=[a,'   !!! Inconsistent rank values in the data.']
       end else a[n_elements(a)-1]=a[n_elements(a)-1]+'  -- OK --'
     ;
     a=[a,'   10 MHz counter from '+strtrim(string(min(C10MHz[i*18:(i+1)*18-1])/1d7, form='(f30.7)'),2)+' to '+$
                                    strtrim(string(max(C10MHz[i*18:(i+1)*18-1])/1d7, form='(f30.7)'),2)]
     dC10MHz=long(C10MHz[i*18+1:(i+1)*18-1]-C10MHz[i*18:(i+1)*18-2])
     a=[a,'      Difference between neigboring packets: ']
       histo,dC10MHz, id[i*18:(i+1)*18-2],  a, vald,fred,point=7
       if (n_elements(vald) ne 1) or (fred[0] ne 17) or (vald[0] ne dtdat) then begin
         err=1
         a=[a,'   !!! Inconsistent values of the 10 MHz counter in the data.']
       end  else a[n_elements(a)-1]=a[n_elements(a)-1]+'  -- OK --'
     if (i+1)*18 lt np then begin
          nextbl=long(C10MHz[(i+1)*18]-C10MHz[(i+1)*18-1])
          a=[a,'      Difference to the following block: '+$
             strtrim(string(double(nextbl)/1d7, form='(f30.7)' ),2)]
       if keyword_set(sl) then $
          if (nextbl le 2l*150109200l-170l*18l)and(nextbl ge 170l*18l)$
          and( packs[(i+1)*18].hzcnt-packs[(i+1)*18-1].hzcnt le 31) $
           then a[n_elements(a)-1]=a[n_elements(a)-1]+'  -- OK for SL --' $
           else begin
            err=1
            a=[a,'   !!! Inconsistent values of the 10 MHz or 1 Hz counter in the data.']
           end $
       else $
          if (nextbl eq dtdat)$
          and( packs[(i+1)*18].hzcnt-packs[(i+1)*18-1].hzcnt le 1) $
           then a[n_elements(a)-1]=a[n_elements(a)-1]+'  -- OK for SA or SR --' $
           else begin
            err=1
            a=[a,'   !!! Inconsistent values of the 10 MHz or 1 Hz counter in the data.']
           end
     end
     ;
     tns_str,min(tns[i*18:(i+1)*18-1]),mintime
     tns_str,max(tns[i*18:(i+1)*18-1]),maxtime
     ;
     a=[a,'   MIU time from '+mintime+' to '+maxtime]
     a=[a,'      Difference between neigboring packets: ']
       histo,(tns[i*18+1:(i+1)*18-1]-tns[i*18:(i+1)*18-2])/1000, $
            id[i*18:(i+1)*18-2],  a, vald,fred,point=6
     if (i+1)*18 lt np then $
          a=[a,'      Difference to the following block: '+$
              strtrim(string(long64(tns[(i+1)*18]-tns[(i+1)*18-1])/1d9,form='(f30.6)' ),2)]
     ;

    end
end






;---------------------------------------------------------------
;---------  Event data: filterbank EB, waveforms EW  -----------
;---------------------------------------------------------------





pro plot_EB, file, nev,ev,evind,evtimes,evtimer,evtypes,evgains,infos,err,fs=fs
  ;
  ; multipage ps overview of all EB events in the file
  ; output goes to a ps file with the same name
  ; y axis autoscaled
  ;
  if nev gt 0 then begin
   plot_init,strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.ps'
   for j =0,nev-1 do begin
     np=evind[j+1]-evind[j]
     if np gt 0 then begin
       x=findgen(np)/fs[j]*1e3
       ymax=1>max(ev[*,evind[j]:evind[j+1]-1],min=ymin)
       if err then errs='(!!)' else errs=''
       plot,[0,1],/nodata,yran=[ymin-(ymax-ymin)/15. ,ymax],xran=[x[0],x[np-1]],xst=1,yst=19,$
               xtit='milliseconds from '+evtimes(j)+' (XO:'+strtrim(string(evtimer[j],form='(f30.9)'),2)+')',$
               ytit='EB event '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nev),2)+$
               '!c !c Filter bank data (TMU)',$
               tit='!d'+infos+$
               '!c!dFile '+file+errs+' Gain='+evgains[j]+' Type='+evtypes[j]+' BLESK-'+ver()+'!n',color=255
        for k =0,11 do oplot,x,ev[k,evind[j]:evind[j+1]-1],color=253-20*(11-k),thick=3
        ytxt=ymin-(ymax-ymin)/16.
        xyouts,x[np/80],ytxt,'Channels:',color=255,charsize=0.7
        for k =0,11 do  xyouts,x[np/15*(k+2)],ytxt,string(k+1),color=253-20*k,charsize=0.7
      end
   end
   plot_close
  end
end



pro plot_EW, file, nev,ev,evind,evtimes,evtimer,evtypes,evgains,infos,err,$
            der=der, int=int,cal=cal,off=off,tit=tit, infmax=infmax,fs=fs,fhp=fhp,$
            thrcsm=thrcsm,maxnsm=maxnsm,tintsm=tintsm

  ;
  ; multipage ps overview of all EW events in the file
  ; output goes to a ps file with the same name
  ; y axis autoscaled
  ; keywords:
  ;           int - integrates waveform, 'int' is added to
  ;                 the file name
  ;           der - derivative of the waveform, 'der' is added to
  ;                 the file name
  ;           cal - calibration coefficient, see function calib_waveform
  ;           off - offset in TMU
  ;           tintsm, default=2e-6 - duration of small intervals in seconds
  ;           thrcsm, default=10   - threshold for max abs dev in units of std
  ;           maxnsm, default=450  - maximum number of smallest waveforms to be plotted

  ;
  if nev gt 0 then begin
   plot_init,strmid(file,0,strlen(file)-4)+'.ps'
   if err then errs='(!!)' else errs=''
   ; whole events first
   if not keyword_set(cal) then cal=1
   if not keyword_set(off) then off=0
   if not keyword_set(tit) then tit='Waveform (TMU)'

   maxp=8000000L ; maximum continuously plotted waveform (100 ms= 8000000 samples)
   maxwf=fltarr(nev)
   for j =0,nev-1 do begin
    nsub=ceil((evind(j+1)-evind(j))/double(maxp))
    isub= (evind(j+1)-evind(j))/nsub
    for i = 0L, nsub-1 do begin
       if nsub gt 1 then plotno=' Plot '+strtrim(string(i+1),2)+'/'+$
               strtrim(string(nsub),2) else plotno=''
       plot_wf,ev(evind(j):evind(j+1)-1),evtimes(j)+$
                ' (XO:'+strtrim(string(evtimer[j],form='(f30.9)'),2)+')',fs=fs[j],fhp=fhp,$
               i*isub,(i+1)*isub-1,$
               ytit='EW event '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nev),2)+plotno+$
               '!c !c Waveform (TMU)',noisescale=20,$
               tit='!d'+infos+$
               ' !d File '+file+errs+' Gain='+evgains[j]+' Type='+evtypes[j]+$
               ' BLESK-'+ver()+'!n'
       plot_spect,ev(evind(j):evind(j+1)-1)+off,i*isub,(i+1)*isub-1,$
               ytit='Power spectrum (dB)',fs=fs[j]
       plot_wf,calib_waveform(ev(evind(j):evind(j+1)-1)+off,cal=cal),evtimes(j)+$
                ' (XO:'+strtrim(string(evtimer[j],form='(f30.9)'),2)+')',fs=fs[j],fhp=fhp,$
               i*isub,(i+1)*isub-1,$
               der=der,int=int,ytit=tit,noisescale=20*cal
    end
    if evind(j) ne evind(j+1) then maxwf[j]=max(abs(ev[evind(j):evind(j+1)-1]+off))
   end
  ; small parts
   for j =0,nev-1 do begin
    pal=long(fs[j]*1e-3) ; parts of 1ms
    for i=0,(evind(j+1)-evind(j))/pal do if evind(j) lt evind(j+1) then begin
       plot_wf,ev(evind(j):evind(j+1)-1),evtimes(j)+$
                ' (XO:'+strtrim(string(evtimer[j],form='(f30.9)'),2)+')',fs=fs[j],fhp=fhp,$
               i*pal,(i+1)*pal-1,noisescale=20,$
               ytit='EW event '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nev),2)+$
              ' Part '+strtrim(string(i+1),2)+'/'+$
               strtrim(string((evind(j+1)-evind(j))/pal+1),2)+$
               '!c !c Waveform (TMU)',$
               tit='!d'+infos+$
               ' !dFile '+file+errs+' Gain='+evgains[j]+' Type='+evtypes[j]+$
               ' BLESK-'+ver()+'!n'
       plot_spect,ev(evind(j):evind(j+1)-1)+off,$
               i*pal,fs=fs[j],$
               min([(i+1)*pal-1,evind(j+1)-evind(j)-1]),$
               ytit='Power spectrum (dB)'
       plot_wf,calib_waveform(ev(evind(j):evind(j+1)-1)+off,cal=cal),evtimes(j)+$
                ' (XO:'+strtrim(string(evtimer[j],form='(f30.9)'),2)+')',fs=fs[j],fhp=fhp,$
               i*pal,(i+1)*pal-1,$
               der=der,int=int,noisescale=20*cal,$
               ytit=tit
    end
   end
  ;
  ; smallest waveforms
    if not keyword_set(thrcsm) then thrcsm=10l ;threshold for max abs dev in units of std
    if not keyword_set(maxnsm) then maxnsm=450; maximum number of smallest waveforms to be plotted
    if not keyword_set(tintsm) then tintsm=2e-6 ; duration of small intervals in seconds
    ;
    for j =0,nev-1 do begin
     pal=long(fs[j]*tintsm) ; parts of tintsm seconds
     nsm=long((evind(j+1)-evind(j))/pal)
     dsm=reform(ev[evind(j):evind(j)+nsm*pal-1],pal,nsm)
     offss=median(dsm) & maxsm=max(abs(dsm-offss),dim=1,indsm)
     thrsm=thrcsm*stdev(dsm) &   partsm=where(maxsm gt thrsm,nprtsm)
     for iprtsm =0, (maxnsm < nprtsm-1) do begin
       plot_wf,ev(evind(j):evind(j+1)-1),evtimes(j)+$
                ' (XO:'+strtrim(string(evtimer[j],form='(f30.9)'),2)+')',fs=fs[j],fhp=fhp,tunit=1e-6,$
               partsm[iprtsm]*pal,(partsm[iprtsm]+1)*pal-1,noisescale=20,$
               ytit='EW event '+strtrim(string(j+1),2)+'/'+strtrim(string(nev),2)+$
              ' part '+strtrim(string(iprtsm+1),2)+'/'+$
               strtrim(string(nprtsm),2)+' |wf-'+strtrim(string(round(offss)),2)+'|>'+$
               strtrim(string(round(thrsm)),2)+'TMU '+'!c !c Waveform (TMU)'
     end
     if nprtsm ge maxnsm then  $
        print,'------------------ more than '+strtrim(string(maxnsm),2)+' '+$
         strtrim(string(pal/fs*1e6,form='(f20.1)'),2)+'-us parts to plot ----------------------'
   ;
   end ; ev loop
   plot_close
   ;
   ;   ew ps file----------------------
   ewpsfile=strmid(file,0,strlen(file)-4)+'_ew.ps'
   plot_init,ewpsfile,yfrac=0.6
   !p.multi=0 &  !p.thick=2 & !x.thick=3 & !y.thick=3
   j=(sort(maxwf))[nev-1]
   if evind(j) ne evind(j+1) then $
        plot_spect,/logf,fs=fs[j],ev[evind(j):evind(j+1)-1]+off,$
               0,evind(j+1)-evind(j)-1, ytit='EW event '+$
               strtrim(string(j+1),2)+'/'+strtrim(string(nev),2)+$
              '!c !c Power spectrum (dB)',$
               tit=strtrim(string((evind(j+1)-evind(j))/long(fs[j]/1e3)),2)+$
               ' ms from '+strmid(evtimes(j),0,27)+'!c!d'+infos+$
               ' !dFile '+file+errs+' Gain='+evgains[j]+' Type='+evtypes[j]+$
               ' BLESK-'+ver()+'!n', pos=[0.1,0.13,0.99,0.88],$
               col=255, xran=[fhp,8e7/2.],freq=freq,spec=spec
   plot_close
   file_copy,ewpsfile,'ew.ps',/overwrite
    ;
    ;txt file--------------------
    ;infmax=string(maxwf[j],j,form='(i5,i2)')
    infmax=string(maxwf[j],form='(i5)')
    ;
    ; spectrogram-----------------

    overviewspect,file, evtimes[j],(evind(j+1)-evind(j))/2.88e11,freq,spec, $
           savfile='bleska_spectrogram.idl', psfile='ewspg.ps',/logy
  end
end




pro overviewspect,file, ts0,tunits=tunits,durh,freq,spec,min_fr_Hz=min_fr_Hz,max_fr_Hz=max_fr_Hz, scaletit=scaletit,$
  savfile=savfile, psfile=psfile,logy=logy,trig=trig,maxtim=maxtim,avetim=avetim,source=source
;
;__accumulated overview figure-----------------
   if not keyword_set(maxtim) then maxtim=4ull*86400000000000ull ; 4 days
   if not keyword_set(avetim) then avetim=40ull  *60000000000ull ; 40 min = 144 time intervals
   if keyword_set(ts0) then str_tns,ts0+'.000.000',tns0 ; newly added time
  ;
  ;           min_fr_Hz  - if set then force the min of the frequency axis
   ;
   ;__ 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 2 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
       if n_elements(mark) ne n_elements(tim) then mark=bytarr( n_elements(tim))
   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(spec) 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])
      markn=bytarr(ntim)
      spgramn[0:ntim-2,*]=spgram[outt:ntim+outt-2,*]
      markn[0:ntim-2]=mark[outt:ntim+outt-2]
      if keyword_set(ts0) then begin
        spgramn[ntim-1,*]=spec
        markn[ntim-1]=keyword_set(trig)
      end  else begin
        spgramn[ntim-1,*]=spgram[ntim+outt-1,*]
        markn[ntim-1]=mark[ntim+outt-1]
      end
      spgram=spgramn
      mark=markn
    ;__ 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
      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(spec)) & spgramn=spgram & tim=tns0 & frq=freq & mark=keyword_set(trig)
   end
; _saving data
   save,file=savfile,tim,frq,spgram,mark
   !p.charsize=0.65
   if not keyword_set(ts0) then begin
     tns0=tim[ntim-1]
     tns_str,tns0,ts0 & durh=0
   end
; _plotting
   tns_str,tim[0],tss
   if not keyword_set(scaletit) then scaletit='Power (dB)'
   if not keyword_set(tunits) then tunits=''
   spgfile=strmid(file,0,strlen(file)-4)+'_'+psfile
   plot_init,spgfile,yfrac=1.
   !p.multi=0 &  !p.thick=2 & !x.thick=2 & !y.thick=2
   if not keyword_set(source) then source=''
   plot_spectrogram,spgramn,-float(long64(tns0[0]-tim)/1000000000ll)/3600.,$
               frq,logy=logy,min_fr_Hz=min_fr_Hz,max_fr_Hz=max_fr_Hz, $
         xtit='Time (hours) from '+strmid(ts0,0,19)+tunits,ztit=scaletit,$
         ytit='Frequency (Hz)',$
         tit=source+' Spectrogram from '+strmid(tss,0,19)+' to '+strmid(ts0,0,19)+tunits,$
         avex=float(long64(avetim)/1000000000ll)/3600., markx=mark
   plot_close
   file_copy,spgfile,psfile,/overwrite
end




pro read_Event,file, XO, apid, ev,nev, evind, evfs,evtimer,evtimes,evtypes,evgains, a, err, confi
;description:
  ; common event reading routine
;   extracts EB filterbank or EW waveform data
;parameters:
;   input:
;       XO  calibration structure
;       file, apid
;       confi - instrument configuration
;   output:
;       nev  - no of events
;       ev    - filterbank data, array[0:11,0:evind(nev)-1])
;               waveform data,   array[0:evind(nev)-1])
;               data from the ith event are: ev[evind[i]:evind[i+1]-1]
;               for empty (patological) events  evind[i]=evind[i+1]
;       evind - array [0,nev] of  [starting
;               indices of events, last index+1]
;
;   evtimes strarr(nev)  first sample of each event from the absolute time
;   evtypes, strarr(nev)
;   evgains  strarr(nev)
;   evfs   dblarr(nev) corrected sampling frequency
;   evtimer dblarr(nev) first sample of each event from the  XO counter, relative
;
  err=0 & nev=0
  read_bin,file,npacks,packs,tns
  check_headers,file,npacks,packs,apid,tns,  t100nsc,  a, err
 if npacks gt 0 then begin

  datver=0 ; final, FM2 with empty parity bit and variable 31bit c10Mz in the event data descriptor
  if packs[0].ver lt 221 then datver=1  ; EM4 or FM, with variable 32bit c10Mz in the event data descriptor
  if packs[0].ver lt 169 then datver=2  ; preEM4, constant 32 bit c10MHz
  datverinfos=(['   Final (FM2) data structure with variable 31bit C10MHz in the event data descriptor',$
  '   Pre-final (EM4-FM1) data structure with variable 32bit C10MHz in the event data descriptor',$
  '   Old (pre-EM4) data structure constant 32bit C10MHz in the event data descriptor']) [datver]
  datvers=(['   LS 31 bits of C10MHz from descriptor (period 215s) ',$
            '   LS 32 bits of C10MHz from descriptor (period 430s) ',$
            '   LS 32 bits of C10MHz from descriptor (one per event) ']) [datver]
  ;
  case apid of
   1228: begin
     pgpk=3ul; packets per page
     pkv=1356ul; values per packet,
     pgv=3468ul; values per page,
               ; 600 values (900 bytes) in the last packet are filled by 0xFF
     dats='EB' ;
     dtdat=113ul*120ul;  time (100ns) between packets
    end
   1229: begin
     pgpk=204ul; packets per page
     pkv=1360ul; values per packet,
     pgv=pgpk*pkv ; values per page, all are used
     dats='EW' ;
     dtdat=170ul;  time (100ns) between packets
    end
  end
  ;
  iev=where(packs.id eq apid, nevp)
  npag=nevp/pgpk & expacks=npacks-npag*pgpk
  a=[a,'','EVENT DATA READER','-- File '+file+', '+strtrim(string(npag),2)+$
    ' '+dats+' pages found. Extra packets: '+strtrim(string(expacks),2)+'.']
  if expacks gt 0 then begin
    err=1; expacks should be 0
    a=[a,'   ERROR -  Extra packets']
  end
  ;
  a=[a,'   Design version: '+ strtrim(string(packs[0].ver, form='(I5)'),2)+' => Data interpreted as:'+datverinfos]
  ;
  if npag gt 0 then begin
    iev=iev(0:npag*pgpk-1)
    get_evd, packs[iev], iev, C10Mhz_, pagno_,evno,evpagno,miuev_,radio_,err,a, datvers
    c1mhz_=256ul*packs[iev].mhzcnt256+packs[iev].mhzcnt
    c1hz_=packs[iev].hzcnt
    gain_=packs[iev].gains_bw/16b
    tns_str,tns[iev],evptimeh_
    reltime, XO, C1hz_,C10Mhz_, rt,corr   ; only corr needed
    case apid of
      1228: begin
         a=[a,'   Absolute time of the 1st EB sample = header abs. time + offset - one page (3.468 ms) - 1.6 us delay']
        tns_str,tns[iev]+long64( (confi.Time_Off_EB*3468000ll-3469600ll)/corr ) ,evptime_
        get_filterbank,packs[iev], npag*pgpk, ev
       end
      1229: begin
         a=[a,'   Absolute time of the 1st EW sample = header abs. time + offset - one page (3.468 ms) - 175 ns delay']
         tns_str,tns[iev]+long64( (confi.Time_Off_EW*3468000ll-3468175ll)/corr ) ,evptime_
         get_waveforms,packs[iev], npag*pgpk, ev
       end
    end
    a=[a,'   Min '+dats+' event (1st sample) abs. time found among all events: '+min(evptime_)]
    a=[a,'   Min '+dats+' header abs. time found among all events            : '+min(evptimeh_)]
    a=[a,'   Max '+dats+' event (1st sample) abs. time found among all events: '+max(evptime_)]
    a=[a,'   Max '+dats+' header abs. time found among all events            : '+max(evptimeh_)]
    a=[a,'   TM values from '+strtrim(string(min(ev)),2)+$
      ' to '+strtrim(string(max(ev)),2)+' median:'+strtrim(string(median(ev)),2)]
    ;
    a=[a,'--','Data structure of separate events (selected and re-ordered from the packet sequence):']
    evo=intarr(pkv*pgpk,npag)
    C10Mhz=ulonarr(pgpk,npag)
    C1Mhz=ulonarr(pgpk,npag)
    C1hz=ulonarr(pgpk,npag)
    pagno=uintarr(pgpk,npag)
    typeev=intarr(pgpk, npag)
    gain=intarr(pgpk, npag)
    evptime=strarr(pgpk, npag)
    ;
    ipag=0l  & evind=0l
    nev=max(evno)
    evtimer=dblarr(nev) & evfs=dblarr(nev)
    evtimes=strarr(nev) & evtypes=strarr(nev) & evgains=strarr(nev)
    for l=1, nev do begin
      evpag=0
      ievp=where(evno eq l,nevp)
      a=[a,'',' Event no. '+strtrim(string(l),2)+' total number of packets: '+$
                strtrim(string(nevp),2)]
      if nevp gt 0 then begin
        for m=min(evpagno[ievp]), max(evpagno[ievp]) do begin
          ievpp=where(evpagno[ievp] eq m,nevpp)
          if nevpp gt 0 then begin
            a=[a,'    Page no. '+strtrim(string(m),2)+':  '+strtrim(string(nevpp),2)+$
                        ' packets ']
            if nevpp eq pgpk then $
              a[n_elements(a)-1]=a[n_elements(a)-1]+' -- OK; '$
            else begin
              err=1
              a[n_elements(a)-1]=a[n_elements(a)-1]+$
              ' -- NOT OK -- (should be '+strtrim(string(pgpk),2)+' packets); '
            end
            if nevpp ge pgpk then begin
              nevppr=nevpp/pgpk
              iorig=ievp[ievpp[0:nevppr*pgpk-1] ]
              ipord=ipag+lindgen(nevppr)
            ;
              evo[*,ipord]=ev[*,iorig]
              C10MHz[*,ipord]=C10Mhz_[iorig]
              C1MHz[*,ipord]=C1Mhz_[iorig]
              C1Hz[*,ipord]=C1hz_[iorig]
              gain[*,ipord]=gain_[iorig]
              pagno[*,ipord]=pagno_[iorig]
              typeev[*,ipord]=miuev_[iorig]+2*radio_[iorig]
              evptime[*,ipord]=evptime_[iorig]
            ;
              ipag=ipag+nevppr
              evpag=evpag+nevppr
        ;
              a[n_elements(a)-1]=a[n_elements(a)-1]+' memory page(s):'
              histo,pagno[*,ipord], iev[iorig], a, vald,fred
              if n_elements(vald) eq 1 then a[n_elements(a)-1]=a[n_elements(a)-1]+'  --OK-- ' $
              else begin
                a=[a,'      -- NOT OK -- (not constant for all packets within event page)']
                err=1
              end
            end
          end
        end
        evind=[evind,ipag*pgpk]
        evpk=evpag*pgpk
        a=[a,'   Event no. '+strtrim(string(l),2)+' pages: '+$
                strtrim(string(evpag),2)+'  ordered packets: '+strtrim(string(evpk),2)]
        ;
        if evpk gt 0 then begin
         a=[a,'   TM values from '+strtrim(string(min(evo[0:pgv-1,evind[l-1]/pgpk:(evind[l]-1)/pgpk])),2)+$
             ' to '+strtrim(string(max(evo[0:pgv-1,evind[l-1]/pgpk:(evind[l]-1)/pgpk])),2)+$
             ' median: '+strtrim(string(median(evo[0:pgv-1,evind[l-1]/pgpk:(evind[l]-1)/pgpk])),2)]
         i1=where(c1mhz[evind[l-1]:evind[l]-1] eq c1mhz[evind[l-1]],n1)
         i1=where(c1hz[evind[l-1]:evind[l]-1] eq c1hz[evind[l-1]],n3)
         if (n1 eq evpk)and(n3 eq evpk)  then $
             a=[a,'   1Hz and 1MHz counters -- OK -- (constant for all packets of this event):', $
                  '          C1Hz='+strtrim(string(c1hz[evind[l-1]]),2)+$
                  ' C1MHz='+strtrim(string( c1mhz[evind[l-1]]/1d6, form='(f30.6)'),2)] $
         else begin
             a=[a,'   1Hz and 1MHz counters -- NOT OK -- (not constant for all packets of this event). ']
             err=1
             C1mhze=lonarr(evpag) & C1hze=lonarr(evpag)
             adds1m='' & adds1=''
             for kk=0,evpag-1 do begin
              i1=where(C1mhz[evind[l-1]+kk*pgpk:evind[l-1]+(kk+1)*pgpk-1] $
                       eq C1mhz[evind[l-1]+kk*pgpk],n1)
              if n1 eq pgpk then C1mhze[kk]=C1mhz[evind[l-1]+kk*pgpk] else begin
                 C1mhze[kk]=-333
                adds1m='( -0.000333 means that differrent values are found within packets of a given page)'
              end
              i1=where(C1hz[evind[l-1]+kk*pgpk:evind[l-1]+(kk+1)*pgpk-1] $
                       eq C1hz[evind[l-1]+kk*pgpk],n1)
              if n1 eq pgpk then C1hze[kk]=C1hz[evind[l-1]+kk*pgpk] else begin
                 C1hze[kk]=-333
                adds1='( -333 means that differrent values are found within packets of a given page)'
              end
             end
             a=[a,'     1-MHz counter: '+adds1m]
             histo,C1mhze, lindgen(evpag), a, vald,fred, point=6,pack='pages'
             a=[a,'     1-Hz counter: '+adds1]
             histo,C1hze, lindgen(evpag), a, vald,fred,pack='pages'
         end
         ;
         if evpk gt 1 then begin
           a=[a,'   Differences of'+datvers+$
                ' in neigboring packets of this event: ']
           dC10MHz=ulong(C10mhz[evind[l-1]+1:evind[l-1]+evpk-1] - C10mhz[evind[l-1]:evind[l-1]+evpk-2])
           histo,dC10mhz, lindgen(evpk-1), a, vald,fred, point=7
                  ok10m=0
           if datver eq 2 $  ; pre em4 design, const c10MHz
            then if (n_elements(fred) eq 1) and (fred[0] eq evpk-1) and (vald[0] eq 0ul)  $
                 then a=[a,'      -- OK -- (constant for all packets of this event):', $
                    '     '+datvers+': '+strtrim(string(C10MHz[evind[l-1]]/1d7, form='(f30.7)'),2)] $
                 else begin
                   err=1 & a=[a,'      -- ERROR -- (not constant for all packets of this event']
                 end  $
            else begin ; new design with changing c10Mhz
             case apid of
                1228: begin
                       if evpag gt 0 then $
                         if (total(dC10MHz[lindgen(evpag)*pgpk] eq 13560l) eq evpag) and $
                            (total(dC10MHz[lindgen(evpag)*pgpk+1] eq 13560l) eq evpag) $
                         then ok10m=1
                       if evpag gt 1 then $
                         if (total(dC10MHz[lindgen(evpag-1)*pgpk+2] eq 7560l) ne evpag-1) then ok10m=0
                      end
                1229: if total(dC10MHz eq 170l) eq (evpk-1)  then ok10m=1
             end
             if ok10m then $
               a=[a,'      -- OK -- correct differences for all packets of this event', $
               '      '+datvers+' in the first packet of this event: '+$
                strtrim(string(C10MHz[evind[l-1]]/1d7, form='(f30.7)'),2)] $
             else begin
               a=[a,'      -- ERROR -- incorrect differences for packets of this event. ']
               err=1
             end
            end
         end
        ;
         types=['Internal','External','_Radio_','Radio']
         i1=where(typeev[evind[l-1]:evind[l]-1] eq typeev[evind[l-1]],n1)
         if (n1 eq evpag*pgpk) then $
           a=[a,'   Event type -- OK -- (constant for all packets of this event): '+$
                     types[typeev[evind[l-1]]] ] $
         else begin
           a=[a,'   Event type -- NOT OK -- (not constant for all packets of this event). First packet: '+$
                   types[typeev[evind[l-1]]] ] & err=1
         end
        ;
         gains=['1x','4x','13x','_13x','1/2x(ant4)','4/2x(ant4)','13/2x(ant4)','_13/2x(ant4)','1/2x(ant3)','4/2x(ant3)','13/2x(ant3)','_13/2x(ant3)','1/15x','4/15x','13/15x','_13/15x']

         i1=where( gain[evind[l-1]:evind[l]-1] eq gain[evind[l-1]],n1)
         if (n1 eq evpag*pgpk) then $
           a=[a,'   Gain -- OK -- (constant for all packets of this event): '+$
                 gains[gain[evind[l-1]] ] ]$
         else begin
           a=[a,'    Gain -- NOT OK -- (not constant for all packets of this event). First packet: '$
               +gains[gain[evind[l-1]] ] ]& err=1
         end
      ;
         i1=where(evptime[evind[l-1]:evind[l]-1] eq evptime[evind[l-1]],n1)
         if (n1 eq evpag*pgpk) then begin
           a=[a,'   Header abs. time  -- OK -- (constant for all packets of this event): '+$
                     evptimeh_[ievp[0]],   $
                '   Event 1st sample abs. time : '+ evptime[evind[l-1]] +')' ]
         end else begin
           a=[a,'   Header abs. time -- NOT OK -- (not constant for all packets of this event). First packet: '+$
                   evptimeh_[ievp[0]],   $
                '   Event 1st sample abs. time: '+evptime[evind[l-1]]  +')' ] & err=1
         end
        ;
         reltime, XO, C1hz[evind[l-1]], C10mhz[evind[l-1]], rt,corr, bits31=(datver eq 0),a=a
         case apid of
          1228: begin
              if datver eq 2 then begin
                 evtimer[l-1]=rt[0]+(3468d-6*(confi.Time_Off_EB-1)-1600d-9)/corr[0]
                 a=[a,'   XO time of the 1st EB sample = data descriptor XO time + offset - one page (3.468 ms) - 1.6 us delay']
              end  else begin
                 evtimer[l-1]=rt[0]-1600d-9/corr[0]
                 a=[a,'   XO time of the 1st EB sample = data descriptor XO time - 1.6 us delay']
              end
              evfs[l-1]=corr[0]/12d-6
            end
          1229: begin
              if datver eq 2 then begin
                evtimer[l-1]=rt[0]+(3468d-6*(confi.Time_Off_EW-1)-175d-9)/corr[0]
                a=[a,'   XO time of the 1st EW sample = data descriptor XO time + offset - one page (3.468 ms) - 175 ns delay']
              end else begin
                evtimer[l-1]=rt[0]-175d-9/corr[0]
                a=[a,'   XO time of the 1st EW sample = data descriptor XO time - 175 ns delay']
              end
              evfs[l-1]=corr[0]*8d7
            end
         end
         evtimes[l-1]=evptime[evind[l-1]]
         evtypes[l-1]=types[typeev[evind[l-1]] ]
         evgains[l-1]=gains[gain [evind[l-1]] ]
;
         a=[a,'   Event 1st sample XO time (reconstructed seconds from XP start): '+string(evtimer[l-1],form='(f25.9)') ]
         a=[a,'   Event data descriptor XO time (reconstructed seconds from XP start)    : '+string(rt[0],form='(f25.9)') ]
         a=[a,'   Sampling frequency correction factor for this event: '+string(corr[0],form='(f15.9)') ]
        end else begin ; if evpk gt 0
          a=[a,'   No valid pages in this event  -- NOT OK -- ']
          err=1
        end
      end else begin ;  if nevp gt 0
          evind=[evind,ipag*pgpk] ; add the same as before, empty event
          a=[a,'   No packets in this event  -- NOT OK -- ']
          err=1
        end
    end ; for l= 1..nev
     ;
  mode=(packs[iev[0]].iinf mod 16) &  tst=((mode ge 10)and(mode le 13))
  if tst then  a=[a,'--','TEST PATTERN MODE DETECTED -> test pattern analysis:'] $
         else  a=[a,'--','Extremes in separate events:']
         ;
  case apid of
   1228: begin
     ev=reform(evo[0:pgv-1,*],12,pgv/12*npag)
     evind=evind/pgpk *(pgv/12)
     for i =0L, nev-1 do if evind[i] lt evind[i+1] then begin
      str_tns,evtimes(i),t0
      tns_str,t0+ulong64((evind[i+1]-evind[i])/evfs[i]*1d9),evends
      a=[a,'',' - Event number: '+strtrim(string(i+1),2)+' from '+evtimes(i)+' to '+evends+$
                    ' (XO: '+ strtrim(string(evtimer[i], form='(f30.9)'),2)+$
                    ' - '+strtrim(string(evtimer[i]+(evind[i+1]-evind[i])/evfs[i], form='(f30.9)'),2)+')']
       if tst then begin
        for j = 0,10 do begin ;__check differences between channels
          devj=ev[j+1,evind[i]:evind[i+1]-1]-ev[j,evind[i]:evind[i+1]-1]
          inok=where( (devj ne 341)and(devj ne -3755), nnok)
          if nnok gt 0 then begin
            tns_str,t0+ulong64(inok[0]/evfs[i]*1d9),tmit
            a=[a,'   Channel '+string(12-j-1,form='(i2.2)')+' - Channel '+string(12-j,form='(i2.2)')+$
                 ': '+strtrim(string(nnok),2)+' ERROR - wrong difference(s), example :'+$
                  strtrim(string(devj[inok[0]]),2)+' at an offset of '+$
                  strtrim(string(inok[0]/evfs[i]*1d6,form='(f230.3)'),2)+' usec -> '+tmit+$
                   ' (XO: '+ strtrim(string(evtimer[i]+inok[0]/evfs[i], form='(f30.9)'),2)+')']
               err=1
          end else a=[a,'   Channel '+string(12-j-1,form='(i2.2)')+' - Channel '+string(12-j,form='(i2.2)')+$
                      ': Differences - OK -']
        end
        devj=ev[0,evind[i]+1:evind[i+1]-1]-ev[0,evind[i]:evind[i+1]-2] ;__check time sequence of channel 1
        inok=where( (devj ne 1)and(devj ne -4095), nnok)
        if nnok eq 0 then a=[a,'   Channel *12* Monotonous time sequence, no irregularities - OK -'] $
        else begin
          a=[a,'   Channel *12* Found '+strtrim(string(nnok),2)+' time sequence irregularities:']
          for ii=0,nnok-1 do begin
            tns_str,t0+ulong64(inok[ii]/evfs[i]*1d9),tmit
            trel=evtimer[i]+(inok[ii]+1)/evfs[i]
            tval=ev[0,evind[i]+1+inok[ii]]
            a=[a,'   Value of '+strtrim(string(tval),2)+$
                 ' succeeding '+strtrim(string(ev[0,evind[i]+inok[ii]]),2)+$
                 ' at offset of '+strtrim(string((inok[ii]+1)/evfs[i]*1d6,form='(f230.3)'),2)+' usec -> '+tmit+$
                 ' (XO: '+ strtrim(string(trel, form='(f30.9)'),2)+')']
          end
          if (nnok gt 1)or (abs(trel-round(trel)) gt 12e-6) or (tval ne 1) then begin
            err=1
            a=[a,'   ERROR:  unexpected irregularities (sequence should be reset only at entire seconds).']
          end  else a=[a,'   Expected reset of the test sequence at entire second - OK -']
         endelse
       end ; __not  pattern: Extremes
        savgolFilter = SAVGOL(4, 4, 1, 2)*evfs[i] ; use 16,16,1,4 for more smoothing
        for j = 0,11 do begin
          tmv=min(ev[j,evind[i]:evind[i+1]-1],tmi)
          tns_str,t0+ulong64(tmi/evfs[i]*1d9),tmit
          a=[a,'   Channel '+string(12-j,form='(i2.2)')+': Min TM value:'+strtrim(string(tmv),2)+' at an offset of '+$
             strtrim(string(tmi/evfs[i]*1d6,form='(f230.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmi/evfs[i], form='(f30.6)'),2)+')']
          tmv=max(ev[j,evind[i]:evind[i+1]-1],tmi)
          tns_str,t0+ulong64(tmi/evfs[i]*1d9),tmit
          a=[a,'               Max TM value:'+strtrim(string(tmv),2)+' at an offset of '+$
               strtrim(string(tmi/evfs[i]*1d6,form='(f230.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmi/evfs[i], form='(f30.6)'),2)+')']
        ;
          savgolFilter = SAVGOL(4, 4, 1, 2)*evfs[i]/1e6 ; use 16,16,1,4 for more smoothing
          mx=max(min=mn, /nan,$
               CONVOL(reform(ev[j,evind[i]:evind[i+1]-1]), savgolFilter, /EDGE_TRUNCATE,/nan),$
               tmx,sub=tmn)
          tns_str,t0+ulong64(tmx/evfs[i]*1d9),tmit
          a=[a,'               Max dTM/dt [1/us]:'+strtrim(string(mx, form='(e30.2)'),2)+' at an offset of '+$
             strtrim(string(tmx/evfs[i]*1d6,form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmx/evfs[i], form='(f30.6)'),2)+')']
          tns_str,t0+ulong64(tmn/evfs[i]*1d9),tmit
          a=[a,'               Min dTM/dt [1/us]:'+strtrim(string(mn, form='(e30.2)'),2)+' at an offset of '+$
             strtrim(string(tmn/evfs[i]*1d6,form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmn/evfs[i], form='(f30.6)'),2)+')']
          a=[a,'               Median TM value:'+strtrim(string(median(ev[j,evind[i]:evind[i+1]-1])),2)]
        end
     end else a=[a,'',' - Event number: '+strtrim(string(i+1),2)+' --EMPTY-- ']
    end
   1229: begin
    ev=reform(evo[0:pgv-1,*],pgv*npag)
    evind=evind/pgpk* pgv
    for i =0L, nev-1 do if evind[i] lt evind[i+1] then begin
      str_tns,evtimes(i),t0
      tns_str,t0+ulong64((evind[i+1]-evind[i])/evfs[i]*1d9),evends
      a=[a,'',' - Event number: '+strtrim(string(i+1),2)+' from '+evtimes(i)+' to '+evends+$
                    ' (XO: '+ strtrim(string(evtimer[i], form='(f30.9)'),2)+$
                    ' - '+strtrim(string(evtimer[i]+(evind[i+1]-evind[i])/evfs[i], form='(f30.9)'),2)+')']
       if tst then begin
        devj=ev[evind[i]+1:evind[i+1]-1]-ev[evind[i]:evind[i+1]-2] ;__check time sequence
        inok=where( (devj ne 1)and(devj ne -4095), nnok)
        if nnok eq 0 then a=[a,'   Monotonous time sequence, no irregularities - OK -'] $
        else begin
          a=[a,'   Found '+strtrim(string(nnok),2)+' time sequence irregularities:']
          trel=dblarr(nnok) & tval=intarr(nnok)
          for ii=0,nnok-1 do begin
            tns_str,t0+ulong64(inok[ii]/evfs[i]*1d9),tmit
            trel[ii]=evtimer[i]+(inok[ii]+1)/evfs[i]
            tval[ii]=ev[evind[i]+1+inok[ii]]
            a=[a,'   Value of '+strtrim(string(tval[ii]),2)+$
                 ' succeeding '+strtrim(string(ev[evind[i]+inok[ii]]),2)+$
                 ' at offset of '+strtrim(string((inok[ii]+1)/evfs[i]*1d6,form='(f230.3)'),2)+' usec -> '+tmit+$
                 ' (XO: '+ strtrim(string(trel[ii], form='(f30.9)'),2)+')']
          end
          ADCcorr=0d-9
          if (nnok eq 2) then $
            if (abs(trel[0]-round(trel[0])+ADCcorr) lt 13d-9)and(abs(trel[0]-trel[1]) lt 13d-9)and(tval[0] eq 0)and(tval[1] eq 0) $
            ; 6 sampling periods previously substracted from evtimer for the ADC delay - 4 need to be added here
            then a=[a,'   Expected reset of the test sequence by two zeros at entire second - OK -'] $
            else begin
              err=1 & a=[a,'   ERROR: unexpected irregularities (sequence should be reset by two zeros at entire second).']
            end  $
          else begin
              err=1 & a=[a,'   ERROR: too many irregularities (sequence should be reset only at entire seconds).']
          end
        endelse
       end ; __not  pattern: Extremes
        tmv=min(ev[evind[i]:evind[i+1]-1],tmi)
        tns_str,t0+ulong64(tmi/evfs[i]*1d9),tmit
        a=[a,'   Min TM value:'+strtrim(string(tmv),2)+' at an offset of '+$
             strtrim(string(tmi/evfs[i]*1d6,form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmi/evfs[i], form='(f30.9)'),2)+')']
        tmv=max(ev[evind[i]:evind[i+1]-1],tmi)
        tns_str,t0+ulong64(tmi/evfs[i]*1d9),tmit
        a=[a,'   Max TM value:'+strtrim(string(tmv),2)+' at an offset of '+$
             strtrim(string(tmi/evfs[i]*1d6,form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmi/evfs[i], form='(f30.9)'),2)+')']
        ;
        savgolFilter = SAVGOL(4, 4, 1, 2)*evfs[i]/1e6 ; use 16,16,1,4 for more smoothing
        mx=max(min=mn, /nan,$
               CONVOL(ev[evind[i]:evind[i+1]-1], savgolFilter, /EDGE_TRUNCATE,/nan),$
               tmx,sub=tmn)
        tns_str,t0+ulong64(tmx/evfs[i]*1d9),tmit
        a=[a,'   Max dTM/dt [1/us]:'+strtrim(string(mx, form='(e10.2)'),2)+' at an offset of '+$
             strtrim(string(tmx/evfs[i]*1d6,form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmx/evfs[i], form='(f30.9)'),2)+')']
        tns_str,t0+ulong64(tmn/evfs[i]*1d9),tmit
        a=[a,'   Min dTM/dt [1/us]:'+strtrim(string(mn, form='(e10.2)'),2)+' at an offset of '+$
             strtrim(string(tmn/evfs[i]*1d6,form='(f30.3)'),2)+' usec -> '+tmit+$
              ' (XO: '+ strtrim(string(evtimer[i]+tmn/evfs[i], form='(f30.9)'),2)+')']
        evmed=median(ev[evind[i]:evind[i+1]-1])
        evstd=stddev(ev[evind[i]:evind[i+1]-1],/double)
        evave=total(ev[evind[i]:evind[i+1]-1],/double)/double(evind[i+1]-evind[i])
        a=[a,'   Med: '+strtrim(string(evmed),2)+$
             '   Ave: '+strtrim(string(evave),2)+$
             '   Std: '+strtrim(string(evstd),2)]
        ;
      end else a=[a,'',' - Event number: '+strtrim(string(i+1),2)+' --EMPTY-- ']
    end
  end
;
    errstr=['-- OK --','-- ERROR --']
    print,'   --> ', strtrim(string(nev),2),' '+dats+' events read.   '+errstr[err]
    ;
  end
 end
 end



pro get_evd, packs, id, C10Mhz, pagno,evno,evpagno,miuev,radio, err,a, datvers
;
; extracts the event data descriptor information
;   last 4 bytes... 31bits (EM4 and newer) 32bits (before EM4) of the 10 Mhz counter ... C10Mhz,
;   time triggered (radi) EW event flag...radio
;   page number in the onboard memory...pagno,
;   event number...evno,
;   page number within the event...evpagno,
;   MIU event flag...miuev
; from
;   packet records...packs
; reports anomalies in
;   anomaly diagnostic string array...a
;   error flag...err
;
    C10Mhz=packs.dat(0)*16777216ul +packs.dat(1)*65536ul+packs.dat(2)*256ul+packs.dat(3)
    radio=packs.dat(4)/128b
    pagno=(packs.dat(4) mod 16b )*1024u+packs.dat(5)*4u+packs.dat(6)/64u
    evno=fix((packs.dat(6) mod 64b)/2b)
    evpagno=fix((packs.dat(6) mod 2b)*128b+packs.dat(7)/2b)
    miuev=fix(packs.dat(7) mod 2b)
    np=n_elements(packs)
    ;
     a=[a,'--','Data structure of all event data packets in the sequence read from file:']
     a=[a,'   Event numbers from '+strtrim(string(min(evno)),2)+' to '+strtrim(string(max(evno)),2)]
       histo,evno, id,  a, vald,fred
     a=[a,datvers+' from '+strtrim(string(min(C10MHz)/1d7, form='(f30.7)'),2)+' to '+$
                                    strtrim(string(max(C10MHz)/1d7, form='(f30.7)'),2)]
     dC10MHz=long(C10MHz[1:np-1]-C10MHz[0:np-2])
     a=[a,'      Difference between the neighboring packets: ']
       histo,dC10MHz, id, a, vald,fred,point=7
     a=[a,'   Memory page number from '+strtrim(string(min(pagno)),2)+' to '+strtrim(string(max(pagno)),2)]
     dpagno=pagno[1:np-1]-pagno[0:np-2]
     a=[a,'      Difference between the neighboring packets: ']
       histo,dpagno, id, a, vald,fred
     a=[a,'   Event page number from '+strtrim(string(min(evpagno)),2)+' to '+strtrim(string(max(evpagno)),2)]
     dpagno=evpagno[1:np-1]-evpagno[0:np-2]
     h=histogram(dpagno,rev=i) & j=where(h gt 0,nj)
     a=[a,'      Difference between the neighboring packets: ']
       histo,dpagno, id,  a, vald,fred
     a=[a,'   MIU event flag ']
       histo,miuev, id,  a, vald,fred
     a=[a,'   Time triggered event flag ']
       histo,radio, id,  a, vald,fred

end







;---------------------------------------------------------------
;------------      Radio Waveform data     RW     --------------
;---------------------------------------------------------------

pro wav_RW, file, rwi,rwq, nrw,rwind, nrwi,rwtnss, rwrate
  ;
  ; writes i-q  2-channel wav files
  ;
  for i=0,nrw-1 do begin
    d=intarr(2,nrwi[i])
    d(0,*)=rwi[rwind[i]:rwind[i+1]-1]
    d(1,*)=rwq[rwind[i]:rwind[i+1]-1]
    wavfile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'_'+$
         strmid(rwtnss[i],11,2)+strmid(rwtnss[i],14,2)+strmid(rwtnss[i],17,2)+'.wav'
      WRITE_WAV, wavfile, d, rwrate[i]
  end
end



pro plot_RW, file,rwi,rwq, nrw,rwind, nrwi,rwtnss,rwtnse, rwbw, rwrate, $
                  amp,fre, rws, rwgain,rwagain,rwdgain,$
                  rwconf,nrwconf
  ;
  ;
  ; multipage ps overview of all RW intervals
  ;
  psfile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.ps'
  plot_init,psfile
  !p.multi=[0,4,4]
  for i=0,nrw-1 do begin
      if rwconf[i].Radio_bandw ne rwbw[i] then bwerr=' (Config NOT OK)' else bwerr=' (Config OK)'
      if rwconf[i].Radio_gain ne rwgain[i] then gerr=' (Config NOT OK)' else gerr=' (Config OK)'
      ;
      i_val=rwi[rwind[i]:rwind[i+1]-1]
      q_val=rwq[rwind[i]:rwind[i+1]-1]
      n_val=rwind[i+1]-rwind[i]
      i_stp=i_val[1:n_val-1]-i_val[0:n_val-2]
      q_stp=q_val[1:n_val-1]-q_val[0:n_val-2]
      if (i_val[0] eq 0)and(q_val[0] eq 0) and $
         ((where(i_stp ne 1))[0] eq -1)and((where(q_stp ne 1))[0] eq -1) $
         then test= '-- 16bit TEST PATTERN --' else test='-- DATA --'
      iamp=amp[rwind[i]:rwind[i+1]-1]
      ifre=fre[rwind[i]:rwind[i+1]-1]
      info=['  RW interval no. '+string(i+1,form='(i2.2)')$
           ,'     I-Q data points: '+strtrim(string(nrwi[i]),2)+' packets: '+strtrim(string(nrwi[i]/512),2) $
           ,'     From (UT) '+rwtnss[i] $
           ,'     To   (UT) '+rwtnse[i] $
           ,'     Center frequency (MHz): '+string(rwconf[i].Radio_fc_MHz,form='(f10.6)')+' MHz' $
           ,'     Bandwidth: '+string(rwrate[i]/1000.,form='(f6.1)')+' kHz'+bwerr $
           ,'     Gain (a/d/r): '+strtrim(string(0+rwagain[i]),2)+$
                    ' / '+strtrim(string(0+rwdgain[i]),2)+' / '+strtrim(string(0+rwgain[i]),2)+gerr $
           ,'     '+test $
           ,'     I Values  (TMU, 1st / min / max / median): '$
           ,'                  '+strtrim(string(i_val[0]),2)+' / '+strtrim(string(min(i_val)),2)+' / '+$
                                 strtrim(string(max([i_val])),2)+' / ' +strtrim(string(median(i_val)),2) $
           ,'     Q Values  (TMU, 1st / min / max / median): '$
           ,'                  '+strtrim(string(q_val[0]),2)+' / '+strtrim(string(min(q_val)),2)+' / '+$
                                 strtrim(string(max([q_val])),2)+' / ' +strtrim(string(median(q_val)),2) $
           ,'     I-Q Amp (TMU, 1st / min / max / median): '$
           ,'                  '+strtrim(string(iamp[0]),2)+' / '+strtrim(string(min(iamp)),2)+' / '$
                            +strtrim(string(max(iamp)),2)+' / ' +strtrim(string(median(iamp)),2) $
           ,'     I-Q Fre (kHz, 1st / min / max / median): '$
           ,'                  '+strtrim(string(ifre[0]),2)+' / '+strtrim(string(min(ifre)),2)+' / '$
                            +strtrim(string(max(ifre)),2)+' / ' +strtrim(string(median(ifre)),2)]
      plot,[0,1],xst=4,ystyle=4,/nodata,color=255
      xyouts,replicate(0,n_elements(info)),1-0.07*indgen(n_elements(info)),info,color=255,charsize=0.45
      ;
      plot_wf, iamp ,rwtnss(i),fs=rwrate[i],$
               ytit=' I-Q Amplitude',noisescale=20,$
               tit='File '+file+'    RW interval '+$
               strtrim(string(i+1),2)+'/'+strtrim(string(nrw),2)
      ;
      plot_wf, ifre,rwtnss(i),fs=rwrate[i],$
               ytit=' I-Q Frequency (kHz)',$
               strtrim(string(i+1),2)+'/'+strtrim(string(nrw),2)
      plot_io,-rwrate[i]/1000./2.+findgen(512)*rwrate[i]/1000./511.,color=255,$  ;  f= -bw*255/256,...0,...+bw
               total(rws[*,rwind[i]/512:rwind[i+1]/512-1],2)/nrwi[i], xtit='Frequency (kHz)',ytit='Power'
   end
   plot_close
end


pro read_RW, apid,file,rwi,rwq, nrw,rwind, nrwi,rwtnss,rwtnse, rwbw, rwrate, $
                  amp,fre, rws, rwgain,rwagain,rwdgain,$
                  rwconf,nrwconf,a,err
;
; extracts RW waveforms (output nrw - number of intervals,
;                        output rwind - array of nrw
;                           starting indices of intervals & last index
;                        output rw - waveform, array[ewind(new-1)])
;  from an array of records of packets (input packs)
;  time of the beginning of each interval
;   is stored in and array of strings (output rwtns)
;   based on time strings of packets (input tns)
;
   rates=[4882.8125d0, 9765.625d0, 19531.25d0 ,78125d0]
;
  read_bin,file,npacks,packs,tns
  check_headers,file,npacks,packs,apid,tns,  t100nsc,  a, err
  ;
  irw=where((packs.id eq 1230)or(packs.id eq 1231) , nrwp)
  if nrwp gt 1 then begin
    ;
    ;__16 bit I-Q from packets
    rwi=intarr(512,nrwp)
    rwq=intarr(512,nrwp)
    rwi(*,*)= 256*packs(irw).dat(lindgen(512)*4)+packs(irw).dat(lindgen(512)*4+1)
    rwq(*,*)= 256*packs(irw).dat(lindgen(512)*4+2)+packs(irw).dat(lindgen(512)*4+3)
    ;
    ;__baseband signal and power spectrum from I-Q
    rwb=intarr(1024,nrwp)
    rws=fltarr(512,nrwp)
    for j=0,nrwp-1 do begin
     ;rrr=poly_fit(findgen(512),rwi(*,j),0,yfit=rwitrend)
     ;rrr=poly_fit(findgen(512),rwq(*,j),0,yfit=rwqtrend)
     rwitrend=total(rwi(*,j))/512.
     rwqtrend=total(rwq(*,j))/512.
     ft=fft(hanning(512)*complex(rwi(*,j)-rwitrend,$
                                 rwq(*,j)-rwqtrend)) ; forward transform I-Q
     fto=rotate([ft(257:511),ft(0:256)],2) ; re-ordered  f= -bw*255/256,...0,...+bw
     rws[*,j]=float(fto)^2+imaginary(fto)^2
     rwb[*,j]=round(fft([complex(0,0),fto,conj(fto[510-indgen(511)])],/inverse)) ;inverse
    end
    ;
    ;__continuous intervals : have the same time in the header
    rwindp=where(tns[irw[1:nrwp-1]] ne tns[irw[0:nrwp-2]],nrw) +1l
    if nrw gt 0 then rwind=[0,rwindp*512l,nrwp*512l] else  rwind=[0,nrwp*512l]
    nrw=nrw+1
    ;
    ;__ bandwidth and gain for each interval
    rwbw=packs(irw(rwind[0:nrw-1]/512l)).gains_bw mod 4b
    rwrate=rates[rwbw]
    rwgain=packs(irw(rwind[0:nrw-1]/512l)).rw_gain mod 16b
    rwagain=packs(irw(rwind[0:nrw-1]/512l)).gains_bw/16b
    rwdgain= (packs(irw(rwind[0:nrw-1]/512l)).gains_bw mod 16b) / 4b
    ;
    ;__ instantaneous amplitude and frequency from analytical I-Q signal
    amp=sqrt(rwi^2.+rwq^2.)
    phas=atan(rwq,rwi)
    dphas=-phas[1:512L*nrwp-1]+phas[0:512L*nrwp-2]
    ineg=where(dphas lt -(!pi),nneg) & if nneg gt 0 then dphas(ineg)=dphas(ineg)+2*!pi
    ipos=where(dphas gt  (!pi),npos) & if npos gt 0 then dphas(ipos)=dphas(ipos)-2*!pi
    fre=[dphas,0.]/2./!pi
    for i=0,nrw-1 do begin
      fre[rwind[i]:rwind[i+1]-2]= fre[rwind[i]:rwind[i+1]-2]*rwrate[i]/1000. ; kHz
      fre[rwind[i+1]-1] =fre[rwind[i+1]-2]
    end
    ;
    ;__start & stop time
    rwtns=tns(irw(rwind[0:nrw-1]/512l))
    tns_str,rwtns,rwtnss ; ___start time
    nrwi=rwind[1:nrw]-rwind[0:nrw-1]
    tns_str,rwtns+long64(1000000000ll*nrwi/long64(rates(rwbw))),rwtnse ;___stoptime
    ;
    print, string(nrw,form='(i2)'),'RW intervals read.'
  endif else begin
    nrw=0
    print,'-----------ERROR: Less than two RW/HFPE packets found --------------'
  end
end

pro fine, file, fr,ave
  ; fr=[-165,-155] frequency range in Hz
  ; ave=16 averaging
  read_rw,1230,file,rwi,rwq, nrw,rwind, nrwi,rwtnss,rwtnse, rwbw, rwrate, $
                  amp,fre, rws, rwgain,rwagain,rwdgain,$
                  rwconf,nrwconf,a,err
 niq=4
 for iq=0,niq-1 do begin
  for i = 0, nrw-1 do begin
    nw= rwind[i+1]-rwind[i] & sh=0l
    nw=nw/niq & sh=nw/niq*iq
    ft=fft(hanning(nw)*complex(rwi[rwind[i]+sh:rwind[i]+sh+nw-1],rwq[rwind[i]+sh:rwind[i]+sh+nw-1]) )
    fto=rotate([ft[nw/2+1:nw-1],ft[0:nw/2]],2)
    sw=float(fto)^2+imaginary(fto)^2
    f=findgen(nw)*float(rwrate[i])/nw - float(rwrate[i])/2.
    af=rebin(f,nw/ave)
    asw=rebin(sw,nw/ave)
    j=where((af gt fr[0])and(af lt fr[1]),nj)
    if nj gt 0 then $
      if (i eq 0)and(iq eq 0) then plot_io,af[j],asw[j],xtit='Hz' $
                else oplot,af[j],asw[j],line=i
  endfor
 endfor
end


;---------------------------------------------------------------
;--------------------  INI and FIN packets  --------------------
;---------------------------------------------------------------


pro read_ini, file,  conf,rwconf,nrwconf,a,err
;
; extracts config from ini packet file
;
 err=0
 read_bin,file,npacks,packs,tns
 check_headers,file,npacks,packs, 1224, tns, t100nsc, a, err
 err_=1
 if (npacks ge 1) and (packs[0].id eq 1224 )  then begin
   if npacks eq 1 then err_=0
   confb=packs[0].dat[0:18]
   rwconfb=packs[0].dat[20+indgen(6)#replicate(1,20)+replicate(1,6)#(8*indgen(20))]
   conf_from_blc,conf,rwconf, confb,rwconfb,nrw=nrw
   conffile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.asc'
   write_conf_asc,conffile,conf,rwconf
   print,'-----------------> INI packet <-------------------'
   a=[a,'','Configuration from the INI packet written to '+conffile]
   bitstr=['OK','ERROR']
     cerr=packs[0].dat[19]/128b & if cerr eq 1 then err=1
   a=[a,'','Configuration check: mode validity (other than 0, 14 or 15): '+bitstr[cerr] ]
     cerr=(packs[0].dat[19]*2b)/128b & if cerr eq 1 then err=1
   a=[a,   'Configuration check: completeness of CONFHFA (19 bytes)    : '+bitstr[cerr] ]
     cerr=(packs[0].dat[19]*4b)/128b & if cerr eq 1 then err=1
   a=[a,   'Configuration check: completeness of HFPE/RW (19+120 bytes): '+bitstr[cerr] ]
   a=[a,   'Configuration check: number of received CONFHFA bytes      : '+strtrim(string(fix(packs[0].dat[19]*8b)/8b),2) ]
   a=[a,   'Configuration check: number of HFPE/RW intervals (max. 20) : '+strtrim(string(fix(packs[0].dat[20+159])),2) ]
 ;  memory test
   memory_test, packs,a,err
 end
 if err_ eq 1 then begin
   err=1
   print,'---> Error reading INI packet'
   a=[a,'','Error reading configuration from the INI packet:',$
          '  number of packets='+strtrim(string(npacks),2)+', APID='+strtrim(string(packs[0].id),2)]
 end
end


pro read_fin, file,  conf,rwconf,nrwconf,a,err
;
; extracts config from ini packet file
;
 err=0
 read_bin,file,npacks,packs,tns
 check_headers,file,npacks,packs,1232,tns, t100nsc, a, err
 if (npacks eq 1) and (packs[0].id eq 1232 )  then begin
   confb=packs[0].dat[0:18]
   rwconfb=packs[0].dat[20+indgen(6)#replicate(1,20)+replicate(1,6)#(8*indgen(20))]
   conf_from_blc,conf,rwconf, confb,rwconfb,nrw=nrw
   conffile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.asc'
   write_conf_asc,conffile,conf,rwconf
   print,'-----------------> FIN packet <-------------------'
   a=[a,'','Configuration from the FIN packet written to '+conffile]
   bitstr=['OK','ERROR']
     cerr=packs[0].dat[19]/128b & if cerr eq 1 then err=1
   a=[a,'','Configuration check: mode validity (other than 0, 14 or 15): '+bitstr[cerr] ]
     cerr=(packs[0].dat[19]*2b)/128b & if cerr eq 1 then err=1
   a=[a,   'Configuration check: completeness of CONFHFA (19 bytes)    : '+bitstr[cerr] ]
     cerr=(packs[0].dat[19]*4b)/128b & if cerr eq 1 then err=1
   a=[a,   'Configuration check: completeness of HFPE/RW (19+120 bytes): '+bitstr[cerr] ]
   a=[a,   'Configuration check: number of received CONFHFA bytes      : '+strtrim(string(fix(packs[0].dat[19]*8b)/8b),2) ]
   a=[a,   'Configuration check: number of HFPE/RW intervals (max. 20) : '+strtrim(string(fix(packs[0].dat[20+159])),2) ]
 ;
 ; memory test
   memory_test, packs,a,err, /fin
 ;
 ; alert statistics
   a=[a,'','Parameters of the first 16 alerts:']
   for i= 0, 15 do begin
     t100ns=strtrim(string( (packs[0].dat[1088+i*60+3]*4294967296ull + packs[0].dat[1088+i*60+4]*16777216ull + $
            packs[0].dat[1088+i*60+5]*65536ull + packs[0].dat[1088+i*60+6]*256ull + $
            packs[0].dat[1088+i*60+7])/1d7, form='(f30.7)'),2)
     ai='' & mi=''
     for j=0,11 do ai=ai+string(packs[0].dat[1088+i*60+8+j*2]*256u + packs[0].dat[1088+i*60+9+j*2],form='(i5)')
     for j=0,11 do mi=mi+string(packs[0].dat[1088+i*60+32+j*2]*256u + packs[0].dat[1088+i*60+33+j*2],form='(i5)')
     R=strtrim(string(packs[0].dat[1088+i*60+56]*16777216ull + packs[0].dat[1088+i*60+57]*65536ull + $
      packs[0].dat[1088+i*60+58]*256ull + packs[0].dat[1088+i*60+59]),2)
     a=[a,' - Alert '+string(i+1,form='(i2.2)')+' 10-MHz-counter = '+t100ns+' Rank = '+R,$
          '            ai = '+ai,'            mi = '+mi]
   endfor
 end else begin
   err=1
   print,'---> Error reading FIN packet'
   a=[a,'','Error reading configuration from the FIN packet:',$
          '  number of packets='+strtrim(string(npacks),2)+', APID='+strtrim(string(packs[0].id),2)]
 end
end


pro memory_test, packs,a, err, fin=fin
 ;
 ; reads memory test data from FIN and INIT packets
 ;
 if keyword_set(fin) then nt=44 else nt=92
 ;
 ; stack errors
   ers=[ 'OK', 'ERROR']
   s10err=packs[0].dat[183]
   s20err=packs[0].dat[181]
   s11err=packs[0].dat[182]
   s21err=packs[0].dat[180]
   if (s10err ne 0) or (s20err ne 0) or (s11err ne 0) or (s21err ne 0) then err=1
   if (s10err ne 0) then s10errs='ERROR-'+string(s10err,form='(z2.2)') else s10errs='OK'
   if (s11err ne 0) then s11errs='ERROR-'+string(s11err,form='(z2.2)') else s11errs='OK'
   if (s20err ne 0) then s20errs='ERROR-'+string(s20err,form='(z2.2)') else s20errs='OK'
   if (s21err ne 0) then s21errs='ERROR-'+string(s21err,form='(z2.2)') else s21errs='OK'

   a=[a,'','Stack1 memory test: writing 0: '+ s10errs+'    writing 1: '+ s11errs,$
           'Stack2 memory test: writing 0: '+ s20errs+'    writing 1: '+ s21errs]
 ;
 ; SDRAM tests
  a=[a,'','SDRAM memory test: '] & sderr=0
  for i= 0,nt do $
   ; test vector: "00000000FFFFFFFF55555555AAAAAAAA ".
   if (packs[0].dat[192+i*20] ne '00'xb) or (packs[0].dat[193+i*20] ne '00'xb) or $
      (packs[0].dat[194+i*20] ne '00'xb) or (packs[0].dat[195+i*20] ne '00'xb) or $
      (packs[0].dat[196+i*20] ne 'FF'xb) or (packs[0].dat[197+i*20] ne 'FF'xb) or $
      (packs[0].dat[198+i*20] ne 'FF'xb) or (packs[0].dat[199+i*20] ne 'FF'xb) or $
      (packs[0].dat[200+i*20] ne '55'xb) or (packs[0].dat[201+i*20] ne '55'xb) or $
      (packs[0].dat[202+i*20] ne '55'xb) or (packs[0].dat[203+i*20] ne '55'xb) or $
      (packs[0].dat[204+i*20] ne 'AA'xb) or (packs[0].dat[205+i*20] ne 'AA'xb) or $
      (packs[0].dat[206+i*20] ne 'AA'xb) or (packs[0].dat[207+i*20] ne 'AA'xb) then $
   begin
     sderr=1
     a=[a,'   ERROR at address '+string(packs[0].dat[188+i*20],form='(Z2.2)')+$
         string(packs[0].dat[189+i*20],form='(Z2.2)')+string(packs[0].dat[190+i*20],form='(Z2.2)')+$
         string(packs[0].dat[191+i*20],form='(Z2.2)')+': Read '+  $
         string(packs[0].dat[192+i*20],form='(Z2.2)')+string(packs[0].dat[193+i*20],form='(Z2.2)')+$
         string(packs[0].dat[194+i*20],form='(Z2.2)')+string(packs[0].dat[195+i*20],form='(Z2.2)')+$
         string(packs[0].dat[196+i*20],form='(Z2.2)')+string(packs[0].dat[197+i*20],form='(Z2.2)')+$
         string(packs[0].dat[198+i*20],form='(Z2.2)')+string(packs[0].dat[199+i*20],form='(Z2.2)')+$
         string(packs[0].dat[200+i*20],form='(Z2.2)')+string(packs[0].dat[201+i*20],form='(Z2.2)')+$
         string(packs[0].dat[202+i*20],form='(Z2.2)')+string(packs[0].dat[203+i*20],form='(Z2.2)')+$
         string(packs[0].dat[204+i*20],form='(Z2.2)')+string(packs[0].dat[205+i*20],form='(Z2.2)')+$
         string(packs[0].dat[206+i*20],form='(Z2.2)')+string(packs[0].dat[207+i*20],form='(Z2.2)')+$
         '  (should be 00000000FFFFFFFF55555555AAAAAAAA)']
   end
  if sderr eq 0 then a[n_elements(a)-1]=a[n_elements(a)-1]+'OK' else err=1
  for i=0,ceil(nt/10.) do begin
    if i eq 0 then adr='Addr ' else adr='     '
    for j=0,10< (nt-10*i) do $
      adr=adr+string(packs[0].dat[188+(i*10+j)*20],form='(Z3.2)')+$
         string(packs[0].dat[189+(i*10+j)*20],form='(Z2.2)')+$
         string(packs[0].dat[190+(i*10+j)*20],form='(Z2.2)')+$
         string(packs[0].dat[191+(i*10+j)*20],form='(Z2.2)')
    a=[a,adr]
  end
end

;---------------------------------------------------------------
;---------              TC configuration             -----------
;---------------------------------------------------------------


pro conv_conf,file
 ;
 ; .blc, .bin <-> .asc
 ;
 ;
 case strmid(file,strpos(file,'.',/REVERSE_SEARCH),strlen(file)) of
 '.bin': begin
    ofile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.asc'
    read_conf_blc,file, conf,rwconf,nrw=nrw
    write_conf_asc,ofile, conf,rwconf,nrw=nrw
    print,'-- Converted '+file+'  ->  '+ofile+string(nrw)+ ' radio windows.'
    file_sim=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.sim'
    write_conf_sim,file_sim, conf,rwconf
    print,'-- Converted '+file+'  ->  '+file_sim
    file_hex=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.hex'
    write_conf_hex,file_hex, conf,rwconf,nrw=nrw+string(nrw)+ ' radio windows.'
    print,'-- Converted '+file+'  ->  '+file_hex
  end
 '.blc': begin
    ofile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.asc'
    read_conf_blc,file, conf,rwconf,nrw=nrw
    write_conf_asc,ofile, conf,rwconf,nrw=nrw
    print,'-- Converted '+file+'  ->  '+ofile+string(nrw)+ ' radio windows.'
    file_sim=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.sim'
    write_conf_sim,file_sim, conf,rwconf
    print,'-- Converted '+file+'  ->  '+file_sim
    file_hex=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.hex'
    write_conf_hex,file_hex, conf,rwconf,nrw=nrw
    print,'-- Converted '+file+'  ->  '+file_hex+string(nrw)+ ' radio windows.'
  end
 '.asc': begin
    ofile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.blc'
    read_conf_asc,file, conf,rwconf,nrw=nrw
    write_conf_blc,ofile, conf,rwconf,nrw=nrw
    print,'-- Converted '+file+'  ->  '+ofile+string(nrw)+ ' radio windows.'
    file_sim=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.sim'
    write_conf_sim,file_sim, conf,rwconf
    print,'-- Converted '+file+'  ->  '+file_sim
    file_hex=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.hex'
    write_conf_hex,file_hex, conf,rwconf,nrw=nrw
    print,'-- Converted '+file+'  ->  '+file_hex+string(nrw)+ ' radio windows.'
  end
 '.hex': begin
    ofile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.blc'
    read_conf_hex,file, conf,rwconf,nrw=nrw
    write_conf_blc,ofile, conf,rwconf,nrw=nrw
    print,'-- Converted '+file+'  ->  '+ofile+string(nrw)+ ' radio windows.'
    file_sim=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.sim'
    write_conf_sim,file_sim, conf,rwconf
    print,'-- Converted '+file+'  ->  '+file_sim
    ofile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.asc'
    write_conf_asc,ofile, conf,rwconf,nrw=nrw
    print,'-- Converted '+file+'  ->  '+ofile+string(nrw)+ ' radio windows.'
  end
  else: print,'-- No conversion for this suffix: '+file
 end
end

pro read_conf_blc,file_blc, conf,rwconf,nrw=nrw
;
; reads a binary file containing conf
;
  confb=bytarr(19) & rwconfb=bytarr(120)
  openr,1,file_blc
   readu,1,confb
   b=0b & c19=1
   for i=0,119 do if (not eof(1)) then begin
     readu,1,b & rwconfb[i]=b & c19=0
   end
  close,1
  if c19 then nrw=0 else nrw=20
  conf_from_blc,conf,rwconf, confb,rwconfb,nrw=nrw
end


pro write_conf_blc,file_blc, conf,rwconf,nrw=nrw
;
; writes a binary file containing conf
;
  conf_into_blc,conf,rwconf, confb,rwconfb
  openw,1,file_blc
  c19=0
  if n_elements(nrw) gt 0 then if nrw eq 0 then c19=1
  if c19 eq 1 $
   then writeu,1,confb $
   else writeu,1,confb,rwconfb
  close,1
end


pro write_conf_sim,file_sim, conf,rwconf
;
; writes an ascii file containing at each line
; a binary-coded byte from conf
;
  conf_into_blc,conf,rwconf, confb,rwconfb
  openw,1,file_sim
   printf,1,confb,rwconfb, format='(B8.8)'
  close,1
end

pro write_conf_hex,file_hex, conf,rwconf,nrw=nrw
;
; writes an ascii file containing a sequence of
; hex-coded bytes from conf
;
  conf_into_blc,conf,rwconf, confb,rwconfb
  openw,1,file_hex
  c19=0
  if n_elements(nrw) gt 0 then if nrw eq 0 then c19=1
  if c19 eq 1 $
   then printf,1,confb, format='(19Z2.2)' $
   else printf,1,confb,rwconfb, format='(139Z2.2)'
  close,1
end


pro read_conf_hex,file_hex, conf,rwconf,nrw=nrw
;
; reads an ascii file containing a sequence of
; hex-coded bytes from conf
;
  confb=bytarr(19) & rwconfb=bytarr(120)
  fi=file_info(file_hex)
  if fi.read then begin
    openr,1,file_hex
    if fi.size ge 38 then $
     if fi.size ge 278 $
       then begin
        readf,1,confb,rwconfb, format='(139Z2.2)'
        nrw=20
       end else begin
         readf,1,confb, format='(19Z2.2)'
         nrw=0
       end
    close,1
  end
  conf_from_blc,conf,rwconf, confb,rwconfb,nrw=nrw
end


pro read_conf_asc,file_asc, conf,rwconf,nrw=nrw
;
; reads an ascii file containing conf
;
  confb=bytarr(19) & rwconfb=bytarr(120)
  conf_from_blc,conf,rwconf, confb,rwconfb,nrw=20  ; get conf structures
  openr,1,file_asc
   readf,1,conf,form='(///33(22x,g9/))'
   i=0
   while (not eof(1)) and (i lt 20) do begin
     rw=rwconf[i]
     readf,1,rw,form='(////4(22x,g9/),22x,g9)'
     rwconf[i]=rw
     i=i+1
   end
  close,1
  if i lt 20 then begin
    rwconf[i].Radio_tstart_s=2047l*2l
    rwconf[i].Radio_tstop_s=2047l*2l
  end
  nrw=i
end



pro write_conf_asc,file_asc, conf,rwconf,nrw=nrw
;
; writes an ascii file containing conf
;
  openw,1,file_asc
   printf,1,'---------- IME-HF / BLESKA configuration file ----------'
   printf,1,'                                               '
   printf,1,'------Parameter------|--Value--|---Comments--->'
   printf,1,'   Instrument mode   |'+string(conf.mode,form='(i9)')+'| 1..BW-RL, 2..BW-AL, 3..BW-AR, 4..RW-W-RL, 5..RW-W-AL, 6..RW-W-AR, 7..PE-W-RL, 8..PE-W-AL, 9..PE-W-AR, 10..TestBW-RL, 11..TestBW-AL, 12..TestRW-W-AR, 13..TestPE-W-RL'
   printf,1,'   Analog gain       |'+string(conf.Analog_Gain,form='(i9)')+'|  0:1x, 1:4x, 2:13x, 3:13x, 4-7: same/2 (ant4), 8-11:same/2 (ant3), 12-15: same/22 (passive)
   printf,1,'   Digital gain      |'+string(conf.Digital_Gain,form='(i9)')+'| 0..1x, 1..2x, 2..4x, 3..1x'
   printf,1,'   Enable int. events|'+string(conf.En_Int_Events,form='(i9)')+'| 0..disable, 1..enable'
   printf,1,'   Enable ext. events|'+string(conf.En_Ext_Events,form='(i9)')+'| 0..disable, 1..enable'
   printf,1,'   Offset SL (4usec) |'+string(conf.Time_Off_SL,form='(i9)')+'| min..-76, max..127; in steps of 4us'
   printf,1,'   Offset EB (pages) |'+string(conf.Time_Off_EB,form='(i9)')+'| min..-128 (depends on EB circ. buff. memory), max..+127 ; page=3.468ms'
   printf,1,'   Offset EW (pages) |'+string(conf.Time_Off_EW,form='(i9)')+'| min..-128 (depends on EW circ. buff. memory), max..+127 ; page=3.468ms'
   printf,1,'   Log2(M_averaged)  |'+string(conf.Lg2_Pnt_Aver,form='(i9)')+'| min..3 (96us), max..31 (7.2h), default..20 (12.6s)'
   printf,1,'   Log2(N_maximum)   |'+string(conf.Lg2_Pnt_Max,form='(i9)')+'| min..1 (24us), max..7 (1.5ms), default..4 (192us)'
   printf,1,' Weight of channel 12|'+string(conf.Weight_C12,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel 11|'+string(conf.Weight_C11,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel 10|'+string(conf.Weight_C10,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  9|'+string(conf.Weight_C9,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  8|'+string(conf.Weight_C8,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  7|'+string(conf.Weight_C7,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  6|'+string(conf.Weight_C6,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  5|'+string(conf.Weight_C5,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  4|'+string(conf.Weight_C4,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  3|'+string(conf.Weight_C3,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  2|'+string(conf.Weight_C2,form='(i9)')+'| min..0, max..15'
   printf,1,' Weight of channel  1|'+string(conf.Weight_C1,form='(i9)')+'| min..0, max..15'
   printf,1,' Event trigger thres.|'+string(conf.Thres_Event,form='(i9)')+'| min..0, max..491520'
   printf,1,'   EMC options       |'+string(conf.EMC_opt,form='(i9)')+'| min..0, max..7, bits(msb->lsb) 2.5V(0..608; 1..625kHz) / 1.2V(0..604; 1..625kHz) / Filterbank-off'
   printf,1,'   Radio test options|'+string(conf.TST_opt,form='(i9)')+'| min..0, max..7, bits(msb->lsb) pattern(11,00:inp12b;01:out16b;10:out8b) / forced-data-ready '
   printf,1,'  Index to EW event 1|'+string(conf.Event_pos_1,form='(i9)')+'| min..0(not recorded), max..31'
   printf,1,'  Index to EW event 2|'+string(conf.Event_pos_2,form='(i9)')+'| min..0(not recorded), max..31'
   printf,1,'  Index to EW event 3|'+string(conf.Event_pos_3,form='(i9)')+'| min..0(not recorded), max..31'
   printf,1,'Event trig. algorithm|'+string(conf.Trig_WF,form='(i9)')+'| min..0(filter_bank), 1(waveform_gain_1x) max..7(waveform_gain_64x)'
   printf,1,'  Number of EB events|'+string(conf.Max_EB_events,form='(i9)')+'| min..0,  max..31, default..31'
   printf,1,'Event duration(pages)|'+string(conf.Event_pages,form='(i9)')+'| min..1, max..117(1EW) or 255(EB only), default..12 ; page=3.468ms'
   printf,1,'EW circ. buff.(pages)|'+string(conf.EW_Circ_buf_pages,form='(i9)')+'| min..0 (EB or RW only) , max..120,  default..113 ; page=3.468ms'
   printf,1,'                                               '
   allwindows=0
   for i=0,19 do begin
     if ((rwconf[i].Radio_tstart_s/2l eq  2047) and ( rwconf[i].Radio_tstop_s/128l eq 31)) then allwindows=1
     if n_elements(nrw) gt 0 then if (i ge nrw) then allwindows=1
     if allwindows eq 0 then begin
       printf,1,'                                               '
       printf,1,'-------- IME-HF / BLESKA Radio interval no. '+string(i+1,form='(i2.2)')+$
                                                     ' ---------'
       printf,1,'                                               '
       printf,1,'------Parameter------|--Value--|---Comments--->'
       printf,1,'Start time - T0 (sec)|'+string(rwconf[i].Radio_tstart_s,form='(i9)')+'| min..0,  max..4092 seconds; 4094: end of sequence of radio intervals'
       printf,1,'Stop  time - T1 (sec)|'+string(rwconf[i].Radio_tstop_s,form='(i9)')+'| min..0,  max..4092 seconds; 4094: Radio triggered EW event at T0'
       printf,1,'  Radio bandwidth    |'+string(rwconf[i].Radio_bandw,form='(i9)')+'| 0..5 kHz, 1..10 kHz, 2..20 kHz, 3..80 kHz'
       printf,1,'  Central freq (MHz) |'+string(rwconf[i].Radio_fc_MHz, form='(f9.6)')+'| min..0, max..40 MHz, step..~38Hz'
       printf,1,'  Radio gain         |'+string(rwconf[i].Radio_gain,form='(i9)')+'| 0..1x, max..7 (128x,5 & 10kHz) ..6(64x,20 kHz) ..5(32x,80 kHz) '
     end
   end
   close,1
end

pro conf_from_blc,conf,rwconf, confb,rwconfb,  nrw=nrw
  ;
  ; transform configuration structure from its binary form
  ;
  conf={mode:0b, Analog_Gain:0b, Digital_Gain:0b, $
      En_Int_Events:0b, En_Ext_Events:0b, $
      Time_Off_SL:0, Time_Off_EB:0, Time_Off_EW:0, $
      Lg2_Pnt_Aver:0b, Lg2_Pnt_Max:0b, $
      Weight_C12:0b, Weight_C11:0b, Weight_C10:0b, Weight_C9:0b, $
      Weight_C8:0b, Weight_C7:0b, Weight_C6:0b, Weight_C5:0b, $
      Weight_C4:0b, Weight_C3:0b, Weight_C2:0b, Weight_C1:0b, $
      Thres_Event:0l, EMC_opt:0b, TST_opt:0b,$
      Event_pos_1:0b, Event_pos_2:0b, Event_pos_3:0b, $
      Trig_WF:0b, Max_EB_events:0b, Event_pages:0b, EW_Circ_buf_pages:0b }
  rwconf=replicate({Radio_tstart_s:0u, Radio_tstop_s:0u, $
                    Radio_bandw:0b, $
                    Radio_fc_MHz:0d0, Radio_gain:0b},20)
;
  conf.mode=confb[0] mod 16b
  conf.Analog_gain=(confb[1] /16b) mod 16b
  conf.Digital_Gain=(confb[1] /4b) mod 4b
  conf.En_Int_Events=(confb[1] /2b) mod 2b
  conf.En_Ext_Events=confb[1] mod 2b
  if confb[2] ge 128 then conf.Time_Off_SL=(-(confb[2] xor 255)-1 )    else conf.Time_Off_SL=confb[2]
  if confb[3] ge 128 then conf.Time_Off_EB=(-(confb[3] xor 255)-1 ) else conf.Time_Off_EB=confb[3]
  if confb[4] ge 128 then conf.Time_Off_EW=(-(confb[4] xor 255)-1 ) else conf.Time_Off_EW=confb[4]
  conf.Lg2_Pnt_Aver=confb[5] / 8b
  conf.Lg2_Pnt_Max=confb[5] mod 8b
  conf.Weight_C12=confb[6] / 16b
  conf.Weight_C11=confb[6] mod 16b
  conf.Weight_C10=confb[7] / 16b
  conf.Weight_C9=confb[7] mod 16b
  conf.Weight_C8=confb[8] / 16b
  conf.Weight_C7=confb[8] mod 16b
  conf.Weight_C6=confb[9] / 16b
  conf.Weight_C5=confb[9] mod 16b
  conf.Weight_C4=confb[10] / 16b
  conf.Weight_C3=confb[10] mod 16b
  conf.Weight_C2=confb[11] / 16b
  conf.Weight_C1=confb[11] mod 16b
  conf.Thres_Event=(confb[12] mod 16b)*2l^(confb[12]/16b)
  conf.Event_pos_1=confb[13] mod 32b
  conf.EMC_opt=confb[13] / 32b
  conf.TST_opt=confb[14] / 32b
  conf.Event_pos_2=confb[14] mod 32b
  conf.Event_pos_3=confb[15] mod 32b
  conf.Trig_WF=confb[16] /32b
  conf.Max_EB_events=confb[16] mod 32b
  conf.Event_pages=confb[17]
  conf.EW_Circ_buf_pages=confb[18]
  crw=1 & if n_elements(nrw) gt 0 then if nrw eq 0 then crw=0
 if crw then begin
  nrw=20
  for i=0,19 do begin
    rwconf[i].Radio_tstart_s=2u*(rwconfb[6*i]*8u+$
                                 rwconfb[6*i+1]/32u)
    rwconf[i].Radio_tstop_s=2u*((rwconfb[6*i+1] mod 32u)*64u+$
                                rwconfb[6*i+2]/4u)

    rwconf[i].Radio_bandw=rwconfb[6*i+2] mod 4b
    rwconf[i].Radio_fc_MHz=80d0/2d0^21*(rwconfb[6*i+3]*4096ul + $
                  rwconfb[6*i+4] * 16ul + rwconfb[6*i+5] / 16ul)
    rwconf[i].Radio_gain=rwconfb[6*i+5] mod 16ul
    if (rwconf[i].Radio_tstart_s/2l eq  2047) and ( rwconf[i].Radio_tstop_s/128l eq 31) then nrw=i
  end
 end
end


pro conf_into_blc,conf,rwconf, confb,rwconfb
  ;
  ; transform configuration structure into a binary form
  ;
  confb=bytarr(19)
  rwconfb=bytarr(120)
  ;
  confb[0]=conf.mode mod 16b
  confb[1]=(conf.Analog_gain mod 16b)*16b +$
           (conf.Digital_gain mod 4b)*4b +$
           (conf.En_Int_Events mod 2b)*2b +$
           conf.En_Ext_Events mod 2b
  confb[2]=byte(round(conf.Time_Off_SL))
  confb[3]=byte(round(conf.Time_Off_EB))
  confb[4]=byte(round(conf.Time_Off_EW))
  confb[5]=(conf.Lg2_Pnt_Aver mod 32b)*8b + $
            conf.Lg2_Pnt_Max mod 8b
  confb[6]=(conf.Weight_C12 mod 16b)*16b +conf.Weight_C11 mod 16b
  confb[7]=(conf.Weight_C10 mod 16b)*16b +conf.Weight_C9 mod 16b
  confb[8]=(conf.Weight_C8 mod 16b)*16b +conf.Weight_C7 mod 16b
  confb[9]=(conf.Weight_C6 mod 16b)*16b +conf.Weight_C5 mod 16b
  confb[10]=(conf.Weight_C4 mod 16b)*16b +conf.Weight_C3 mod 16b
  confb[11]=(conf.Weight_C2 mod 16b)*16b +conf.Weight_C1 mod 16b
    if conf.Thres_Event le 0 then Thres_Event_e=0b $
         else Thres_Event_e=alog(conf.Thres_Event)/alog(2.)-3
    if Thres_Event_e lt 0. then Thres_Event_e=0b $
              else Thres_Event_e=byte(Thres_Event_e)
    Thres_Event_m = byte( conf.Thres_Event / 2.^Thres_Event_e)
  confb[12]=(Thres_Event_e mod 16b)*16b+Thres_Event_m mod 16b
  confb[13]=(conf.EMC_opt mod 8b)*32b+(conf.Event_pos_1 mod 32b)
  confb[14]=(conf.TST_opt mod 8b)*32b+(conf.Event_pos_2 mod 32b)
  confb[15]=conf.Event_pos_3 mod 32b
  confb[16]=(conf.Trig_WF mod 8b)*32b + conf.Max_EB_events mod 32b
  confb[17]=byte(round(conf.Event_pages) )
  confb[18]=byte(round(conf.EW_Circ_buf_pages) )
  ;
  for i=0,19 do begin
    rwconfb[6*i]=byte(rwconf[i].Radio_tstart_s/16b)
    rwconfb[6*i+1]=(rwconf[i].Radio_tstart_s/2b  mod 8b)*32b+$
                (rwconf[i].Radio_tstop_s/128b)  mod 32b
    rwconfb[6*i+2]=(rwconf[i].Radio_tstop_s/2b  mod 64b)*4b+$
                rwconf[i].Radio_bandw  mod 4b
      Radio_fc=ulong(round(rwconf[i].Radio_fc_MHz*2d0^21/80d0,/l64))
    rwconfb[6*i+3]= byte(Radio_fc/4096ul)
    rwconfb[6*i+4]= byte((Radio_fc/16ul) mod 256ul)
    rwconfb[6*i+5]= byte(Radio_fc mod 16ul)*16b + $
                        rwconf[i].Radio_gain mod 16b
  end
end


;---------------------------------------------------------------
;------------     general data reading routines    -------------
;---------------------------------------------------------------


pro read_bin,file,npacks,packs,tns,oldabstime=oldabstime
;
; reads a bin file (input: file)
; into an array of records (output: packs)
; of TARANIS packets (number od packtes - output: npacks)
; output: tns...header time: modified jul. date in nanoseconds
; ;
; each packet is a record
  pack={ $
;___header 27 byte
  id:0u,no:0u,totno:0u,hzcnt:0u,mhzcnt256:0u,mhzcnt:0b,$
  Year:0u, Month:0b, Day:0b, Hour:0b, Minute:0b, Sec:0b, usec16:0u, usec_ind:0b,$
  Orb:0u, Orbinf:0b,iinf:0b,ianom:0b, ver:0b,$
;___HK params 13 byte
  tenMHzcnt256:0ul,tenMHzcnt:0b,tb:0b, t3:0b ,t4:0b,tf:0b, miu_ev:0b, hfa_al:0b, rw_gain:0b, gains_bw:0b,$
;___data
  dat:bytarr(2048)}
  npacks=0
  packs=pack
;
  inf=file_info(file)
  if inf.read then begin
;
;   number of packets
    npacks=inf.size/2088
   if npacks gt 0 then begin
    packs=replicate(pack,npacks)
;
;   read
    openr,l,file,/get_lun,/SWAP_IF_LITTLE_ENDIAN
    readu,l,packs
    close,l
    free_lun,l
;
;    print file info
    print,'-- File '+file+' read  '+strtrim(string(npacks),2)+' packets found.'
   if  ARG_PRESENT(tns) then begin
    if keyword_set(oldabstime) then $
      tns=into_tns(packs.usec16*256000ul + (packs.usec_ind)*1000ul,packs.Sec,packs.Minute,$
                packs.Hour,packs.Day,packs.Month,packs.Year>1) $
    else $
      tns=into_tns(16000ul*packs.usec16 + 1000ul*(packs.usec_ind/16b),packs.Sec,packs.Minute,$
                packs.Hour,packs.Day,packs.Month,packs.Year>1)
    tns_str,min(tns),mintns &  print,'   Min packet time: '+mintns
    tns_str,max(tns),maxtns &  print,'   Max packet time: '+maxtns
   end
  endif else print,'-- File '+file+' has zero packets.'
 endif else begin
    print,'-- File '+file+' does not exist or is not readable.'
    npacks=0
  endelse
end


pro get_filterbank,packs, npacks, d
  pkv=1356; values in a packet
  ;
    d=uintarr(pkv,npacks)
    ; nondivided bytes
    d(*,*)= packs.dat(8+(lindgen(pkv)+1)*3/2-1)
    ; adding 4bits from the right
    d(lindgen(pkv/2)*2,*)=d(lindgen(pkv/2)*2,*)*16+ $
                     packs.dat(8+ lindgen(pkv/2) *3+1 )/16
    ; adding 4bits from the left
    d(lindgen(pkv/2)*2+1,*)=d(lindgen(pkv/2)*2+1,*)+ $
              256* (packs.dat(8+ lindgen(pkv/2) *3+1 ) mod 16)
end

pro get_waveforms,packs, npacks, d
  pkv=1360; values in a packet
  ;
    d=intarr(pkv,npacks)
    ; nondivided bytes
    d(*,*)= packs.dat(8+(lindgen(pkv)+1)*3/2-1)
    ; adding 4bits from the right
    d(lindgen(pkv/2)*2,*)=d(lindgen(pkv/2)*2,*)*16+ $
                     packs.dat(8+ lindgen(pkv/2) *3+1 )/16
    ; adding 4bits from the left
    d(lindgen(pkv/2)*2+1,*)=d(lindgen(pkv/2)*2+1,*)+ $
              256* (packs.dat(8+ lindgen(pkv/2) *3+1 ) mod 16)
   ; adding sign
    ineg=where(d and 2048,nneg)
    if nneg gt 0 then d(ineg)=-(d(ineg) xor 4095)-1
end




pro reltime, XO, sec, c10MHz_32, rt,corr,bits31=bits31,a=a
 ;
 ; returns time in seconds rel to XP start, based on onboard XO corrected by 1sec pulses
 ;
 ; XO...input structure for XO acuracy corrections: n, sec, c10MHz, corr (obtained ideally from SR or SA)
 ; sec... input scalar or array of seconds
 ; c10MHz_32... input scalar or array of 32 ls bits of 10 MHz counter
 ;    if bits31 keyword is set then 31 ls bits are assumed
 ; rt... returns rel time in seconds, a double array of the same size
 ; corr ...returns correction factor for sampling frequency (fs=8e7*corr), a double array of the same size
 ;
    rt=-333d0 & corr=1d0
    if (XO.n gt 0)and(n_elements(sec) gt 0)and(n_elements(c10MHz_32) eq n_elements(sec)) then begin
     rt=dblarr(n_elements(sec)) & corr=dblarr(n_elements(sec))+1d0
     for i=0L,n_elements(sec)-1L do begin
      is=(sort(abs(long(XO.sec) - sec[i])))[0]
      if i eq 0 then is0=is
      corr[i]=1d0/XO.corr[is]
      if keyword_set(bits31) then begin
        infs=' 31bit '
        c32sec=long64(XO.c10MHz[is] mod 2ll^31)
        ; if c32 already overflown, but c32 sec not yet
        if (c10MHz_32[i] lt c32sec)and(c10MHz_32[i] lt 2ll^27)and(c32sec gt 2ll^31-2ll^27)  $
         then c32=long64(c10MHz_32[i])+2ll^31 else c32=long64(c10MHz_32[i])
      end else begin
        infs=' 32bit '
        c32sec=long64(XO.c10MHz[is] mod 2ll^32)
        ; if c32 already overflown, but c32 sec not yet
        if (c10MHz_32[i] lt c32sec)and(c10MHz_32[i] lt 2ll^27)and(c32sec gt 2ll^32-2ll^27)  $
         then c32=long64(c10MHz_32[i])+2ll^32 else c32=long64(c10MHz_32[i])
      end
      rt[i]=double(XO.sec[is]) + double(c32-c32sec)/ (corr[i]*1d7)
     end
      if keyword_set(a) then a=[a,'   --- Rel. time correction for '+strtrim(string(n_elements(sec)),2)+$
       ' points based on '+strtrim(string(XO.n),2)+' XO calibration points at 1Hz counter from '+$
       strtrim(string(min(XO.sec)),2)+ ' to '+strtrim(string(max(XO.sec)),2),$
       '   --- First point: 1Hz counter = '+strtrim(string(sec[0]),2)+','+infs+'10MHz counter = '+$
       strtrim(string(c10MHz_32[0]/1d7,form='(f30.9)'),2)+' corrected by 40bit 10MHz counter of '+$
       strtrim(string(XO.c10MHz[is0]/1d7,form='(f30.9)'),2)+' obtained at 1Hz counter of '+$
       strtrim(string(XO.sec[is0]),2)+'. Result: '+strtrim(string(rt[0],form='(f30.9)'),2)]
    end else if keyword_set(a) then a=[a,'   --- Rel. time correction failed for '+strtrim(string(n_elements(sec)),2)+$
       ' points and '+strtrim(string(XO.n),2)+' XO calibration points.']
end




pro check_headers,file,npacks,packs, apid, tns, t100nsc, a, err, XO=XO
    ;
    ; if set: XO...output structure for XO acuracy corrections: n, sec, c10MHz, corr
    ; a ... string array
    ; err...error flag
    ;t100nsc... 10 MHz counter at the last 1 Hz pulse
    ;

  a=['-----------------> BLESK ver.'+ver()+' run on '+systime(/utc)+' UTC','-- File '+file+', '+strtrim(string(npacks),2)+$
    '  packets found.']
 if npacks gt 0 then begin
  a=[a, 'APID:']
  id=lindgen(npacks)
  histo,[packs.id], id, a, vald,fred, text=[replicate('not used',1224),$
    'INIT','SA','SR','SL','EB','EW','HFPE','RW','FIN', replicate('not used',64303)]
  if ((n_elements(vald) ne 1) or (fred[0] ne npacks) or (vald[0] ne apid)) then begin
   err=1 & a=[a,'!!! Inconsistent APID numbers: should be = '+strtrim(string(apid),2)+' in all '+strtrim(string(npacks),2)+' packets.']
  end
  a=[a,'Mode:']
  histo,[packs.iinf mod 16], id, a, vald,fred, text=['Standby','EBEW-SRSL','EBEW-SASL','EBEW-SASR',$
         'RW-EW-SRSL','RW-EW-SASL','RW-EW-SASR',$
         'HFPE-EW-SRSL','HFPE-EW-SASL','HFPE-EW-SASR',$
         'TestEBEW-SRSL','TestEBEW-SASL','TestRW-EW-SASR','TestHFPE-EW-SRSL','!wrong-mode','!wrong-mode']
  if ((n_elements(vald) ne 1) or (fred[0] ne npacks) or (vald[0] ge 14) ) then begin
   err=1 & a=[a,'!!!ERROR: Must be constant and <14 in all '+strtrim(string(npacks),2)+' packets.']
  end
  ;
  a=[a,'Total number of packets: '] & histo,packs.totno, id,  a, vald,fred
    if ((n_elements(vald) ne 1) or (fred[0] ne npacks)) then begin
     err=1 & a=[a,'!!!ERROR: Must be constant in all '+strtrim(string(npacks),2)+' packets.']
    end
  ;
  a=[a,'Packet No from '+strtrim(string(min(packs.no)),2)+' to '+strtrim(string(max(packs.no)),2)]
  if (min(packs.no) ne 1 ) then begin
     err=1 & a=[a,'!!!ERROR: Minimum packet number must be = 1']
    end
  if (max(packs.no) ne packs[0].totno) then begin
     err=1 & a=[a,'!!!ERROR: Maximum packet number must be = '+strtrim(string(packs[0].totno),2)]
    end
  if npacks gt 1 then begin
    dno= 0+packs[1:npacks-1].no - packs[0:npacks-2].no
    if ((apid eq 1229)or(apid eq 1228)) then begin
      ibadstep=where((dno ne 1),nbadstep)
    end else begin  
      ibadstep=where((dno ne 1)and(dno ne -packs[0:npacks-2].totno+1),nbadstep)
    endelse 
    if nbadstep gt 0 then begin
      a=[a,'   !!!!WARNING: Unexpected sequence of packet numbers detected.',$
           '   -> TRYING TO REORDER THE SEQUENCE ACCORDING TO THE PACKET NUMBER']
      ireo=sort(packs.no)
      packs=packs[ireo] & tns=tns[ireo]
      dno= 0+packs[1:npacks-1].no - packs[0:npacks-2].no
      ibadstep=where((dno ne 1),nbadstep)
      if nbadstep gt 0 then begin
        err=1 & a=[a,'!!!ERROR: Reordering unsuccesful. Number of non-unitary steps = '+strtrim(string(nbadstep),2)]
        a=[a,'   Differences of packet numbers in neigboring packets: ']
        histo,dno, id, a, vald,fred
      end else a=[a,'   -> REORDERING SUCCESSFUL: A SEQUENCE WITH INCREMENT 1 OBTAINED']
    end else a=[a,'   Sequence of packet numbers -- OK -- ']
  end
  if (packs[0].no ne 1 ) then begin
     err=1 & a=[a,'!!!ERROR: First packet number must be = 1']
   end
  if (packs[npacks-1].no ne packs[npacks-1].totno) then begin
     err=1 & a=[a,'!!!ERROR: Last packet number must be = '+strtrim(string(packs[npacks-1].totno),2)]
   end
    ;
  mhzc=256ull*packs.mhzcnt256+packs.mhzcnt
  tmsc=1000000ull*packs.hzcnt+mhzc
  a=[a,'Hz and MHz counters from '+ strtrim(string(1d-6*min(tmsc), form='(f30.6)'),2) +$
     ' to '+strtrim(string(1d-6*max(tmsc), form='(f30.6)'),2)]
  a=[a,   '-each separately from '+$
     strtrim(string(min(packs.hzcnt)),2)+'s to '+strtrim(string(max(packs.hzcnt)),2)+$
     's and from '+strtrim(string(min(mhzc)),2)+'us to '+strtrim(string(max(mhzc)),2)+'us.']
  if npacks gt 1 then begin
   dtms= long64(tmsc[1:npacks-1] - tmsc[0:npacks-2])
   a=[a,'   Differences in neigboring packets: ']
   histo,dtms, id, a, vald,fred,point=6
  end
   ;
  tns_str,min(tns),mintns
  tns_str,max(tns),maxtns
  a=[a,'MIU packet time from '+mintns+' to '+maxtns]
  if npacks gt 1 then begin
   dtms= long64(tns[1:npacks-1] - tns[0:npacks-2])/1000ll
   a=[a,'   Differences in neigboring packets: ']
   histo,dtms, id,  a, vald,fred,point=6
  end
   ;
  a=[a,'MIU orbit number: ']
  histo,[packs.orb], id,  a, vald,fred
    ;
  a=[a,'MIU orbit information / directon: ']
  histo,[packs.orbinf/128b], id,  a, vald,fred,text=['DOWN','UP']
  a=[a,'MIU orbit information / number of events: ']
  histo,[packs.orbinf mod 64], id, a, vald,fred
  ;
  a=[a,'Waveform overflows in 306us: ']
  histo,[packs.iinf/16b], id, a, vald,fred, text=['none','<16','<64','<128','<256','<512','<1024','>2048']
  ;
  a=[a,'Instrument anomalies: ']
  anom=['PE or RW memory full','No "data ready" signal from ADC','No "second" and/or "microsecond" clock signal from FPGA XP',$
        'Inconsistent event memory pages','Latch-up in antenna 3','Latch-up in antenna 4','Survey memory overflow','Not used']
  ans=strarr(256) & ans[0]='OK' & for k=0b,255 do for l=0b,7 do if (k and 2b^l) ne 0 then ans[k]=ans[k]+'--'+anom[l]
  histo,[packs.ianom], id, a, vald,fred,text=ans
   if (n_elements(vald) gt 1) or (vald[0] ne 0) then begin
     err=1
     a=[a,'      -> ERROR']
   end
  ;
  a=[a,'FPGA configuration version: ']
  histo,[packs.ver], id,  a, vald,fred
  if ( (n_elements(vald) ne 1) or (fred[0] ne npacks)) then begin
   err=1 & a=[a,'!!!ERROR: Must be constant in all '+strtrim(string(npacks),2)+' packets.']
  end
  ;
  t100nsc=256ull*packs.tenMHzcnt256+packs.tenmhzcnt
  a=[a,'10 MHz counter at the last 1Hz pulse from '+ strtrim(string(1d-7*min(t100nsc), form='(f30.7)'),2) +$
     ' to '+strtrim(string(1d-7*max(t100nsc), form='(f30.7)'),2)]
;---
;  histo,[t100nsc], id,  a, vald,fred,sfunc='hex_sfunc'
;---
 if npacks gt 1 then begin
  a=[a,'   Differences in neigboring packets: ']
  dt100ns= long64(t100nsc[1:npacks-1] - t100nsc[0:npacks-2])
  dts= packs[1:npacks-1].hzcnt - packs[0:npacks-2].hzcnt
  histo,[dt100ns], id,  a, vald,fred,point=7
  inz=where(dt100ns gt 0,nnz)
  if keyword_set(XO) then begin
      XO={n:0} ; default, no data
      if nnz gt 1 then begin
         XO={n:nnz, sec:packs[inz].hzcnt, c10MHz:t100nsc[inz], corr:double(dts[inz])*1d7/double(dt100ns[inz])}
         ppmz=dt100ns[inz]/dts[inz]-10000000ll
         a=[a,'   XO PPM: ']
         histo,[long(1d7/XO.corr-1d7)], id[inz],  a, vald,fred,point=1
      end else  a=[a, '-ppm sequence is not obtained']
  end
 end else if keyword_set(XO) then begin
                   XO={n:0} ; default, no data
                   a=[a, '-ppm sequence is not obtained']
          end
;
  a=[a,'Temperature FPGA: ']
  histo,[packs.tf], id,  a, vald,fred,sfunc='tfpga_sfunc'
   a=[a,'Temperature XO: ']
  histo,[packs.tb], id, a, vald,fred,sfunc='txo_sfunc'
  a=[a,'Temperature ant 3: ']
  histo,[packs.t3], id, a, vald,fred,sfunc='tsens_sfunc'
  a=[a,'Temperature ant 4: ']
  histo,[packs.t4], id,  a, vald,fred,sfunc='tsens_sfunc'
  ;
  a=[a,'Number of external events: ']
  histo,[packs.miu_ev], id,a, vald,fred
  a=[a,'Number of generated alerts: ']
  histo,[packs.hfa_al], id,  a, vald,fred
  a=[a,'Survey memory error: ']
  histo,[(packs.rw_gain / 128b) mod 2], id,  a, vald,fred
   if (n_elements(vald) gt 1) or (vald[0] ne 0) then begin
     err=1
     a=[a,'      -> ERROR']
   end
  a=[a,'Event memory error: ']
  histo,[(packs.rw_gain / 64b) mod 2], id, a, vald,fred
   if (n_elements(vald) gt 1) or (vald[0] ne 0) then begin
     err=1
     a=[a,'      -> ERROR']
   end
  a=[a,'Configuration error: ']
  histo,[(packs.rw_gain / 32b) mod 2], id, a, vald,fred
   if (n_elements(vald) gt 1) or (vald[0] ne 0) then begin
     err=1
     a=[a,'      -> ERROR']
   end
  a=[a,'SDRAM CRC error: ']
  histo,[(packs.rw_gain / 16b) mod 2], id, a, vald,fred
   if (n_elements(vald) gt 1) or (vald[0] ne 0) then begin
     err=1
     a=[a,'      -> ERROR: SDRAM CRC']
   end
  a=[a,'Radio gain: ']
  histo,[packs.rw_gain mod 16b], id,  a, vald,fred
  a=[a,'Radio bandwidth: ']
  histo,[packs.gains_bw mod 4], id, a, vald,fred
  a=[a,'Digital gain: ']
  histo,[(packs.gains_bw mod 16b)/4], id, a, vald,fred
  a=[a,'Analog gain: ']
  histo,[packs.gains_bw/16], id, a, vald,fred
 end
end

pro write_rep,file,a,err,global_err
  openw,1,strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.txt'
    errstr=['-- OK --','-- ERROR --']
    if err eq 1 then global_err=1
    printf,1,errstr[err]
    for i=0,n_elements(a)-1 do printf,1,a[i]
  close,1
end

function tfpga_sfunc,b
  if b le 127 then tdegC = ((-b * 1.3) + 55) else tdegC = (((256-b) * 1.3) + 55)
  return, string(b,form='(Z2.2)')+' -> '+strtrim(string(tdegC,form='(f30.1)'),2)+' degC'
end

function txo_sfunc,b
  tdegC=b* 1.608 - 273.15
  return, string(b,form='(Z2.2)')+' -> '+strtrim(string(tdegC,form='(f30.1)'),2)+' degC'
end

function tsens_sfunc,b
  x=float(b)/32.
  tdegC= +115.858 -262.94*x +338.667*x^2 -239.3074*x^3 $
          +92.6344*x^4 -19.7001*x^5 +2.1562*x^6 -0.09486*x^7
  return, string(b,form='(Z2.2)')+' -> '+strtrim(string(tdegC,form='(f30.1)'),2)+' degC'
end

function hex_sfunc,b
  return, strtrim(string(b,form='(Z30)'),2)
end


pro histo,d, id,a, vald,fred, text=text, point=point,pack=pack,sfunct=sfunct
  ;
  ; calculates a histogram of values with backward indices
  ; d.. data, already selected from an original array at indices id
  ; a..output string array
  ; vald.. output array of values
  ; fred ..output array of corresponding frequencies
  ; text .. array of explanation strings for every value, to be placed between ()
  ; point.. show as floating point
  ; sfunct.. function to transform values into an alternative string
  ;
    rest=n_elements(d)
    ind=lindgen(rest)
    if not keyword_set(pack) then pack='packets'
    while rest gt 0 do begin
      j=where(d[ind] eq d[ind[0]],nj)
      if rest eq n_elements(d) then begin ; first pass
        vald=d[0] & fred=nj
      end else begin
        vald=[ vald,d[ind[0]] ] & fred=[ fred,nj ]
      end
      if nj eq n_elements(d) then  nos=' (all '+pack+')' else begin
        nos=': '
        for l=0,17< (nj-1) do $
         nos=nos+strtrim(string(id[ind[ j[l] ]] + 1),2)+', '
        if 17 lt nj-1 then nos=nos+'...' else nos=strmid(nos,0,strlen(nos)-2)
      end
      if keyword_set(text) then explan=' ('+text[d[ind[0]]]+')' else explan=''
      if keyword_set(sfunct) then funct=call_function(sfunct,d[ind[0]])  $
                             else begin
                                if keyword_set(point) then begin
                                  form='(f30.'+strtrim(string(point),2)+')'
                                  data=d[ind[0]]/10d0^point
                                end else begin
                                  form='(I30)' & data=d[ind[0]]
                                end
                                funct=strtrim(string(data, form=form),2)
                             end
      a=[a, '      '+funct+explan+' in '+strtrim(string(nj),2)+' '+pack+nos]
      in=where( d[ind] ne d[ind[0]], rest)
      if rest ne 0 then ind=ind[in]
    end
end



pro debug, file, event=event
  read_bin,file,npacks,packs
    C10Mhz=packs.dat(0)*16777216ul +packs.dat(1)*65536ul+packs.dat(2)*256ul+packs.dat(3)
    pagno=packs.dat(4)*1024u+packs.dat(5)*4u+packs.dat(6)/64u
    evno=fix((packs.dat(6) mod 64b)/2b)
    evpagno=fix((packs.dat(6) mod 2b)*128b+packs.dat(7)/2b)
    miuev=fix(packs.dat(7) mod 2b)
   nfile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'_debug.txt'
   openw,1,nfile
   if keyword_set(event) then $
     for i=0,n_elements(packs)-1 do $
       printf,1,i+1,packs[i].id, packs[i].no, packs[i].totno, packs[i].hzcnt,$
                packs[i].mhzcnt, packs[i].dat[0:3],pagno[i],evno[i],evpagno[i],miuev[i],$
              form='(i6,"|",i4,4("|",i6),"||",4Z3.2,"|",I5,I3,I4,I2,"|")' $
   else $
     for i=0,n_elements(packs)-1 do $
       printf,1,i,packs[i].id, packs[i].no, packs[i].totno, packs[i].dat[0:7], $
              form='(i6,"|",i4,"|",i6,"|",i6,"||",4Z3.2,"|",4Z3.2,"|")'
   close,1
  print,'-> File '+nfile+' generated.'
end


;---------------------------------------------------------------
;------------     BLESKA log file, time headers    -------------
;---------------------------------------------------------------



pro read_log,logfile, t0,infos,err
  ; read bleska log file
  ; to get start time string t0
  ;
  s=file_info(logfile)
  if s.read then begin
   openr,1,logfile
   a='' & t0='' & lat=-3333. & lon=-3333.
   while not eof(1) do begin
    readf,1,a
    if strmid(a,0,24) eq '>XPStart: Started OK at ' then t0=strmid(a,24,23)
    if strmid(a,0,5) eq ' POS:' then reads,strmid(a,5,strlen(a)-8),lat,lon,alt
   end
   close,1
   if lat ge 0 then lats='N' else lats='S'
   if lon ge 0 then lons='E' else lons='W'
   infos='Start: '+t0+'UT, Pos: '+strtrim(string(abs(lat),form='(f6.2)'),2)+lats+$
        ' '+strtrim(string(abs(lon),form='(f7.2)'),2)+lons; +'  Alt: '+strtrim(string(alt,form='(I60)'),2)+'m'
   print
   if t0 ne '' then  print,'Start time UT: ',t0,' found in ',logfile $
               else  print,'No start time found in ',logfile
  end else print,'Unable to read'+logfile
end

pro fill_time,file,t0
; fills time into records instead of mexic, saves as .dat
;
  read_bin,file,npacks,packs
  if (npacks gt 0) and (t0 ne '') then begin
    str_tns, t0+'.000.000',tns0
    Hz_cnt=packs.hzcnt
    MHz_cnt=256ull*packs.mhzcnt256+packs.mhzcnt
    tns=tns0(0)+1000000000ull*Hz_cnt+1000ull*MHz_cnt
    from_tns, tns, nanosecond,second,minute,hour,day,month,year
    ;
    packs.usec16=nanosecond / 16000ul
    packs.usec_ind= ((nanosecond/1000ul) mod 16ul)*16ul
    packs.Year=year & packs.Month=month & packs.Day=day
    packs.Hour=hour & packs.Minute=minute & packs.Sec= second
    ;
    packs.Orb=60000 & packs.Orbinf=0
    ;
    nfile=strmid(file,0,strpos(file,'.',/REVERSE_SEARCH))+'.dat'
    openw,1,nfile,/SWAP_IF_LITTLE_ENDIAN
    writeu,1,packs
    close,1
    print,'-- Time filled in '+strtrim(string(npacks),2)+' packets -> '+nfile
      tns_str,min(tns),mintns
      print,'    Min time:'+string(min(Hz_cnt),form='(i6)')+'s'+$
              string(min(MHz_cnt),form='(i9)')+'us -> '+mintns
      tns_str,max(tns),maxtns
      print,'    Max time:'+string(max(Hz_cnt),form='(i6)')+'s'+$
               string(max(MHz_cnt),form='(i9)')+'us -> '+maxtns
end else print,'-- File',file,' doesnt exist or t0 is undefined.'
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
 ;
 ; open lanscape ps plot
 ; close by device,/close
 ;
 common plot_ini,yfraction
 set_plot,'ps'
 if not keyword_set(yfrac) then yfrac=1.
 yfraction=yfrac
 device,file=file,xsize=25.,ysize=18.*yfrac,/landscape,xoffset=1.,yoffset=28.,$
             bits_per_pixel=8,/color, /helvetica
 !p.font=0
 !p.charsize=1.
 !p.multi=[0,3,3]
 !p.thick=1
 !x.thick=1
 !y.thick=1
 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'
end


pro plot_wf,wf,time,i1,i2, der=der, int=int, ytit=ytit, tit=tit, pos=pos,$
            yran=yran, cent=cent, tunits=tunits,noisescale=noisescale,fs=fs,fhp=fhp,$
            out_t=out_t,out_wf=out_wf, out_start=out_start,out_yran=out_yran,$
            notimeaxis=notimeaxis,noyaxis=noyaxis,noerase=noerase, overplotcolor=overplotcolor,overplotsym=overplotsym
;
; plots a  waveform wf
;  selects data from index i1 to index i2 (if absent plots all data)
; keywords: int - time integral  ---- 1st order high- pass filtered at fhp
;           der - time derivative
;           ytit - y title
;           tit - top title
;           yran - y axis range, autoscaled if absent
;           cent - center x axis in microsec at 0 between i1 and i2
;                  if not set:  x axis in ms, annotated
;                  by 't (milliseconds) from '+time
;           tunits - in seconds, changes default units on the time axis, (1e-6 ..us)
;           noisescale - if set, this minimum scale is used when autoscale
;           fs - sampling frequency in Hz
;           fhp - corner frequency of the 1st order high-pass filter applied after integration, in Hz
;           out_t,out_wf,out_yran,out_start  - plotted waveform
;                          (after integration or derivative but before compression, if applicable)
;			notimeaxis,noyaxis,noerase - switches to combine plots
;           overplotcolor - switch to overplot an existing b&w plot with a given color
;           overplotsym - to overplot with a symbol
;
  if not keyword_set(fs) then fs=8e7 ; sampling frequency
  if keyword_set(fhp) then nsm=round(fs/fhp) else nsm=0
       ;  width of the inverse boxcar high-pass filter applied after integration
;
  n=n_elements(wf)
 if n gt 0 then begin
  y=wf
  ;
  ; plot all data if not selected
  if not keyword_set(i2) then i2=n-1
  if not keyword_set(i1) then i1=0
  ;

  ; 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 begin
       y=[reform(y),replicate(!values.f_nan,i2-n+1)]
       n=i2+1
      end
  ;
  ; axis and plot annotations
  if not keyword_set(tunits) then $
     if keyword_set(cent) then tunits=1e-6 else tunits=1e-3
  if tunits eq 1e-9 then tunitss='ns' else $
   if tunits eq 1e-6 then tunitss='!mm!3s' else $
    if tunits eq 1e-3 then tunitss='ms' else $
      if tunits eq 1 then tunitss='s' else $
         if tunits eq 60 then tunitss='minutes' else $
            tunitss=strtrim(string(tunits),2)+' seconds'
  xtit='Time ('+tunitss+') from '+time
  if keyword_set(cent) then begin
    x=(dindgen(i2-i1+1)-(i2-i1+1)/2.)/fs/tunits
  end else begin
      x=(dindgen(i2-i1+1)+i1)/fs/tunits
  end
  if not keyword_set(ytit) then ytit='y (a.u.)'
  if not keyword_set(tit) then tit=''
  ;
  common plot_ini,yfraction
  if keyword_set(pos) then begin
    !y.ticklen=0.005/(pos(2)-pos(0))
    !x.ticklen=0.005/(pos(3)-pos(1))/yfraction
  end else begin
    !y.ticklen=0.01
    !x.ticklen=0.01/yfraction
  end
;
  ; integration
  if keyword_set(int) then begin
;     y=Ir_filter([1.,-1.],[whp*(1./whp+1.),1.], $
;      total(y(i1:i2)-total(y(i1:i2),/nan)/float(i2-i1+1),/CUMULATIVE,/nan)/fs ,/doub)$
;     ; cumulative sum of all data up to the current index, high pass filtered at fhp: whp=fhp*2.*!pi/fs

     ;_ cumulative sum of all data up to the current index, added nsp/2 on both sides, if possible
     if i1 gt nsm/2 then begin
       ii1=i1-nsm/2 & iii1=nsm/2
     end else begin
       ii1=0 & iii1=i1
     end
     if i2 lt n-1-nsm/2 then ii2=i2+nsm/2 else ii2=n-1
     ;_ integration, offset substracted
     y=total(y[ii1:ii2]-total(double(y[ii1:ii2]),/nan)/double(ii2-ii1+1),/CUMULATIVE,/nan)/fs
     ;_ high pass filtered at fhp
     if (nsm gt 1)and(nsm lt n_elements(y)) then y=y-smooth(y,nsm,/nan,/EDGE_TRUNCATE)
     y=y[iii1:iii1+i2-i1]
  end else if keyword_set(der) then begin
     savgolFilter = SAVGOL(4, 4, 1, 2)*fs ; use 16,16,1,4 for more smoothing
     y=CONVOL(y[i1:i2], savgolFilter, /EDGE_TRUNCATE,/nan)
     ;y=deriv(y(i1:i2))*fs ; for raw derivatives, little smoothing
  end else y= y[i1:i2]  ;  copy the input data
  ;
  ;wf output
    out_t=x & out_wf=y & out_start=xtit
  ;
  ; 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)) then begin
   totmaxy=max(y, min=totminy,/nan)
   if keyword_set(int) or keyword_set(der)  $
      or (not keyword_set(noisescale)) then noisesc=0 else noisesc=noisescale
   yran=[-noisesc<totminy*1.05,noisesc>totmaxy*1.05 ]
  end
  out_yran=yran
  if keyword_set(notimeaxis) then begin
    xtickname=replicate(' ',30)
    xti=''
  end else begin
    xtickname='' & xti=xtit
  end
  if keyword_set(noyaxis) then begin
    ytickname=replicate(' ',30)
    yti=''
  end else begin
    ytickname='' & yti=ytit
  end
  ;
  if keyword_set(overplotcolor) or  keyword_set(overplotsym) then begin
    if not keyword_set(overplotcolor) then overplotcolor=255
    if keyword_set(overplotsym) then begin
      psym=8
      usersym,overplotsym*[-1,0,1,0],overplotsym*[0,-1,0,1],color=overplotcolor,/fill
    end else psym=0
    oplot,x,y, color=overplotcolor,psym=psym
  end else plot,x,y, yrange=yran,ytitle=yti+compr,xtitle=xti,title=tit,xstyle=1,color=255,$
       pos=pos,xtickname=xtickname,ytickname=ytickname,noerase=noerase
 endif ; waveform not empty
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,i1,i2, int=int,der=der,ytit=ytit,tit=tit,pos=pos,$
   yran=yran,xran=xran, fs=fs,logf=logf,col=col,overplot=overplot, $
   dataonly=dataonly,psd=psd, freq=freq, spec=spec
;
; plots a power spectrum of a waveform wf
;  selects data from index i1 to index i2 (if absent plots all data)
; keywords: int - time integral
;           der - time derivative
;           ytit - y title
;           tit - top title
;           yran - y axis range, autoscaled if absent
  ;
  ; 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
  if not keyword_set(ytit) then psdtit='Power (dB)' else psdtit=ytit
  if not keyword_set(tit) then tit=''
  if not keyword_set(col) then col=255
  if not keyword_set(fs) then fs=8e7 ; sampling frequency
    xtit='Frequency (Hz)'
    fscale=1
  if not keyword_set(logf) then $
   if fs gt 1e6 then begin
    xtit='Frequency (MHz)'
    fscale=1e-6
   end else if fs gt 1e3 then begin
    xtit='Frequency (kHz)'
    fscale=1e-3
   end

  common plot_ini,yfraction
  if keyword_set(pos) then begin
    !y.ticklen=0.005/(pos(2)-pos(0))
    !x.ticklen=0.005/(pos(3)-pos(1))/yfraction
  end else begin
    !y.ticklen=0.01
    !x.ticklen=0.01/yfraction
  end


  ; substraction of the DC offset
  y= wf(i1:i2)-total(double(wf(i1:i2)))/double(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/4    ; final number of points in the averaged spectrum
   nn=7 ; final no of averaged ffts
  end else begin
   nsmooth=1; moving average
   npoi=1024 ; final frequency average of the spectrum to max npoi points, + 5 point moving average
   nfft=32 ; no of averaged ffts will be from nfft to 2*nfft,
           ; nfft must be power of 2  (Hanning window, approx 50% overlap)
   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[long(i*(n2-ovl)):long(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.
 ;
 ; integration and derivative done in the frequency domain

   if keyword_set(int) then s=s/(2*!pi*x/fscale)^2
   if keyword_set(der) then s=s*(2*!pi*x/fscale)^2
  ;
  ; plotting

   if keyword_set(logf) then begin
     nl=1024 < nspec/2
     logav,x,s,nspec-1,xl,sl,nl
     x=xl &s=sl
     xform='logticks'
   end else xform=''
  ; plotting
  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.)
       yform=''
     end
   if not keyword_set(dataonly) then $
    if keyword_set(overplot) then oplot,freq,spec,color=col $
    else plot,freq,spec,yrange=yran,xtickform=xform,$
     ylog=keyword_set(psd),ytickform=yform,$
     ytitle=psdtit,xtitle=xtit,xtickn=xtickn, $
     pos=pos,title=tit,xstyle=1,xlog=logf,color=col
  end
end



pro plot_spectrogram,sp,x,y,xtit=xtit,ytit=ytit,ztit=ztit,tit=tit,min_fr_Hz=min_fr_Hz,max_fr_Hz=max_fr_Hz,$
   zran=zran,logy=logy,logz=logz,nox=nox,pos=pos,avex=avex,markx=markx
;
; plots a power spectrogram of a matrix sp
;           ytit - y title
;           tit - top title
;           yran - y axis range
;           xave - averaging interval on x
;           markx - mark x axis by triangles, boolean vector, same size as x
;           min_fr_Hz,max_fr_Hz  - if set then force the min, max of the frequency axis 
  ; 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)'
  ;
  if not keyword_set(xtit) then xtit='Time (h)'
  ;
 bottom=64
 if not keyword_set(pos) then pos=[0.07,0.08,0.91,0.95]
 posi=[pos[0],pos[1],pos[2],pos[3]]
 poss=[pos[2]+0.070,pos[1],pos[2]+0.077,pos[3]]
 ;
 ; -------------- color scale
 if logz then begin
   scalez=alog10(zran)
   zform='logticks'
 end  else begin
   scalez=zran
   zform=''
 end
;
 !x.ticklen=-0.0001 & !y.ticklen=-0.5
 plot,[0,1],zran, /noerase,/nodata,ylog=logz,$
      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
 ;
 ; ------------------ time axis corrections
 if keyword_set(nox) then begin
   xtits=''
   xtickn=replicate(' ',30)
 end else begin
   xtits=xtit
   xtickn=''
 end
 ;
 if not keyword_set(markx) then markx=bytarr(n_elements(x))
 ;
 siz=size(sp) & if siz[0] ne 2 then navey=1 else navey=siz[2]
 minx=min(x, max=maxx,/nan)
 if keyword_set(avex) then begin
   navex=long((maxx-minx)/avex)+1 & avex_=float(maxx-minx)/float(navex)
   if navex gt 3000l then navex=3000l
   sp_=replicate(!values.f_nan,navex,navey)
   markx_=bytarr(navex)
   x_=minx+(findgen(navex)+0.5)*avex_
   for iave=0,navex-1 do begin
    indave=where( (x ge minx+iave*avex_)and(x le minx+(iave+1)*avex_), nindave)
    if (nindave gt 0) then begin
     if total(/nan,markx[indave]) gt 0 then markx_[iave]=1
     ifini=where(finite(sp[indave,0]),nfini)
     if nfini gt 0 then sp_[iave,*]= reform(total(/nan,sp[indave,*],1))/float(nfini)
    end
   end
 end else begin
   sp_=sp
   markx_=markx
   x_=x
 end
 ;
 ; --------------------frequency axis corrections
 miny=min(y, max=maxy,/nan)
 if n_elements(min_fr_Hz) eq 1 then miny=min_fr_Hz
 if n_elements(max_fr_Hz) eq 1 then maxy=max_fr_Hz
 siz=size(sp_) & nx=siz[1] & if siz[0] ne 2 then ny=1 else ny=siz[2]
 if logy then  begin
    dy=min( alog10(y[1:ny-1]/y[0:ny-2])) & nny=ceil(alog10(maxy/miny)/dy)
    v=alog10(miny)+alog10(maxy/miny)/(nny-1)*findgen(nny)
    sp__=replicate(!values.f_nan,nx,nny)
    for ix=0,nx-1 do sp__[ix,*]=interpol(sp_[ix,*],alog10(y),v)
  endif else begin
    dy=min(y[1:ny-1]-y[0:ny-2]) & nny=ceil((maxy-miny)/dy)
    v=miny+(maxy-miny)/(nny-1)*findgen(nny)
    sp__=replicate(!values.f_nan,nx,nny)
    for ix=0,nx-1 do sp__[ix,*]=interpol(sp_[ix,*],y,v)
  endelse
 ;
 ; --------------------- plotting
 plot,/noerase,[minx,maxx],[miny,maxy],/nodata,$
      ylog=logy,xst=4,yst=4,pos=posi
 TVImage,sp__,posi, zran(0), zran(1),$
         top=253,bottom=bottom,invalid=0
 common plot_ini,yfraction
    !y.ticklen=-0.006/(pos(2)-pos(0))
    !x.ticklen=-0.006/(pos(3)-pos(1))/yfraction
 plot,/noerase,[minx,maxx],[miny,maxy],tit=tit,/nodata,$
      xtit=xtits,xtickn=xtickn,ytit=ytit,ystyle=1,ylog=logy,$
      pos=posi,ythick=3, xthick=3,xstyle=1,color=255
 for i= 0, n_elements(markx_)-1 do if markx_[i] then plots,x_[i],miny,psym=5,color=255,thick=3
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





;---------------------------------
;processing routines.
;----------------------------------



pro ascii_output,ascfile,datf,id,block_ev,nd_nev,dtimes,dgains,t0,t1,tunits,fs,fhp,d,$
                 wf=wf, noint=noint
 case  id of
   '1227':dtext='SL block '+strtrim(string(block_ev),2)+'/'+strtrim(string(nd_nev),2)
   '1229':dtext='EW event '+strtrim(string(block_ev),2)+'/'+strtrim(string(nd_nev),2)
   else:dtext='UNKNOWN'
 end
   tit='TM File '+datf+id+'.dat  Gain='+dgains[block_ev-1]+' BLESK-'+ver()
   ytit='Integrated waveform (nT)'
   if keyword_set(noint) then int=0 else int=1
   ;
   plot_init,ascfile+'.ps'
   plot_wf,d,dtimes(block_ev-1),round(8e7*tunits*t0),round(8e7*tunits*t1),$
               tunits=tunits,int=int,der=0,fs=fs,fhp=fhp, ytit=dtext+$
               '!c !c '+ytit,noisescale=20,$
               tit=' !d'+tit+'!n', out_t=t, out_wf=wf
  plot_close
  nwf=n_elements(wf)
  openw,1,ascfile+'.txt'
  printf,1,'ASCII data file extracted from '+tit
  printf,1,dtext+' :   '+strtrim(string(nwf),2)+'  data points '
  str_tns, dtimes(block_ev-1),tns0
  tns1=tns0+ulong((1d9*tunits)*t[0]+5d-1) &  tns_str,  tns1,tstart
  printf,1,'Start time : '+tstart
  tns2=tns0+ulong((1d9*tunits)*t[nwf-1]+5d-1) &  tns_str,  tns2,tstop
  printf,1,'Stop  time : '+tstop
  printf,1,ytit+', high-pass inv. boxcar filter at '+strtrim(string(fhp),2)+$
             ' Hz, sampled at '+strtrim(string(fs),2)+' Hz'
  printf,1,wf
  close,1
end



pro martin20120830
 off=4.9 ; correction of digitization offset in A/D units
 cal=1.28e6   ; offset corrected A/D units -> nT/s;   nT/s= cal*(A/D +offs)
 dir='C:\sant\taranis\bleska\data\20120828_rustrel\martin'
 datf='20120830_033104_'
   fs=8e7 ; sampling
   fhp=2e3 ; high pass
  ;
  cd,dir
  read_ini,datf+'1224.dat',  conf,rwconf,nrwconf
    ;
  ; ____SL
  read_Survey,datf+'1227.dat', ppm, 1227, d,nd_nev, dind,dtimes,dgains, a,err, conf
;
  block_ev=11
  t0=49 & t1=119 & tunits=1e-6; microsec
  dd=calib_waveform(d[dind[block_ev-1]:dind[block_ev]-1]+off,cal=cal)
  ascii_output,'martin20120830_sl',datf,'1227',block_ev,nd_nev,dtimes,dgains,t0,t1,tunits,fs,fhp,dd,$
      wf=wf1
   ;
   ; ____EW
  read_Event,datf+'1229.dat', ppm, 1229, d,nd_nev, dind,dtimes,dtypes,dgains, a, err, conf
  ; one page correction -BAD COMMANDING
  str_tns, dtimes,tns & tns=tns+3468000ul &  tns_str,  tns,dtimes
  ;
   block_ev=1
  dd=calib_waveform(d[dind[block_ev-1]:dind[block_ev]-1]+off,cal=cal)
    t0=10.026d0 & t1=11.026d0 & tunits=1e-3; msec
   ascii_output,'martin20120830_ew',datf,'1229',block_ev,nd_nev,dtimes,dgains,t0,t1,tunits,fs,fhp,dd,$
     wf=wf2
   ;
   str_tns, dtimes,tns & tns=tns+ulong(t0*tunits*1d9+5d-1) &  tns_str,  tns,dtimes
   t0=0.270d0 & t1=0.340d0 & tunits=1e-3; msec
    ascii_output,'martin20120830_ew_det',datf,'1229',block_ev,nd_nev,dtimes,dgains,t0,t1,tunits,fs,fhp,wf2,$
     noint=1
 end

