;--------------
; processing of handyscope files
;   compile: blesk.pro, rdh.pro
;   run: process, '*.bin' or process,'filename.bin'
;-------------

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



pro rdh,file
  ;
  ; reads one handyscope  file and an associated auxilliary txt file
  ; outputs *_sp#.ps and *_overview#.ps
  ; copies them to sp#.ps and overview#.ps
  ; adds text information to overview.txt
  ; adds aux information necessary for generation of overview#.ps to overview#.idl
  ;
  ;
      durus=52 ; duration of trigger detail in us
      fhp=2.99e3 ; frequency of the inverse boxcar high-pass filter applied after EW integration, in Hz
           ; lower cutoff frequency for spectra
      ifus=2. ; duration of interferometry trigger detail in us, additional inverse boxcar high-pass filter
             ; at 2/ifus is used for intergrated interferometry signal
      ifint=0 ; =1 if integrated waveform is used for interferometry
   ;
   ; Lomnicky stit
   ; calibration since 2016/08/06
      gaindB=[23.,29.,35.,0.,41.,0.,0.,0.,47.,30.]
      AntSurface = 0.228 ; = 0.435 * 0.525 m^2
      nant=4 ; number of antennas

      azim=[325.,55.] ; degrees, azimuth of the ant plane, 1,2, to the right when looking along dB/dt>0
      sep=89; interferometry, separation of antennas in metres
      sepazim=-17. ; degrees, interferometry separation azimuth 1->2
;
      coordinates='43.9410N,5.4836E'
  ;
  ; old settings
       ; old calibration till 2013/08/30, without antialiasing filters
       ; heliport position, close perp antennas
       ;
      azim_1308=[49.,-41] ; degrees, azimuth of the ant plane, 1,2
      sep_1308=1; interferometry, separation of antennas in metres
      sepazim_1308=94. ; degrees, interferometry separation azimuth 1->2
   ;           0=30dB,1=36dB,2=42dB,4=48dB,8=54dB,  9=37dB loGain;
      gaindB_1308=[30.,36.,42.,0.,48.,0.,0.,0.,54.,37.]
  ;
  ; old settings
   ; calibration till 2013/09/26
   ; wall 30m interferometry position
       ;
      azim_1309=[49.,-41] ; degrees, azimuth of the ant plane, 1,2
      sep_1309=1; interferometry, separation of antennas in metres
      sepazim_1309=94. ; degrees, interferometry separation azimuth 1->2
   ;           0=30dB,1=36dB,2=42dB,4=48dB,8=54dB,  9=37dB loGain;
      gaindB_1309=[30.,36.,42.,0.,48.,0.,0.,0.,54.,37.]
  ;
  ; _ auxilliary txt file
   rdth,file,trigtimes,waitt, $; trig time string, wait time in seconds
              nsam, $  ; number of samples in each channel
              postsam, $; number of post-trigger samples in each channel
              terr,$ ; time accuracy
              t1,t2,gain,poe,fs,$; sampling fr in Hz
              bits,$
              range,$ ;fltarr[2] oscilloscope range in each channel +-range V
              trigl0,trigl1,$ ;fltarr[2] trig level in each channel, fraction of -range..range
              timeout ;0 if triggered
   ;
   ;_trigtime -> first sample
  str_tns, trigtimes,trigtime
  tns_str, trigtime-double(nsam-postsam)/fs*1d9, times
   ;
   ;_trigtime -> hist settings
  times_1308='2013-08-30T20:00:00.000.000.000'
  times_1309='2013-09-26T18:00:00.000.000.000'
  str_tns,times_1308 ,time_1308
  str_tns,times_1309 ,time_1309
  case 1 of
   trigtime lt time_1308: begin
      azim=azim_1308
      sep=sep_1308
      sepazim=sepazim_1308
      gaindB=gaindB_1308
     print,'      using settings valid before '+strmid(times_1308,0,10)
    end
   (trigtime lt time_1309)and (trigtime ge time_1308): begin
      azim=azim_1309
      sep=sep_1309
      sepazim=sepazim_1309
      gaindB=gaindB_1309
     print,'      using settings valid before '+strmid(times_1309,0,10)
    end
   else: begin
     print,'      current settings used.'
   end
  endcase
  ;
  ;_data reading
  a=uintarr(nsam,nant)
  openr,1,file &  readu,1,a &  close,1
  wf=float(a)-2.^(bits-1) & a=0
  cal=fltarr(nant)
  ;
  ;_antenna cycle
  intwf=fltarr(nsam,nant)
  for i=0,1 do begin
   ;_calibration  A/D units -> nT/s;   nT/s= cal*(A/D +offs)
    cal[i]=range[i]/2d0^(bits-1)/AntSurface/1d1^(gaindB[gain]/20.)*1d9
   ; _plot each separately
    plot_one_antenna, reform(wf[*,i]),  cal[i], intwfi, $
    	i, azim[i], file, fs, fhp,times, durus,$ ; ant no, azimuth, file,...
       trigtimes,waitt,nsam,postsam,terr,$ ; configuration...
       t1,t2,gaindB[gain],poe,bits,range[i],trigl0[i],trigl1[i],timeout
    intwf[*,i]=intwfi
  end

  for i=2,3 do begin
    plot_one_sevan, reform(wf[*,i]),i,file,fs,fhp,times,nsam,postsam,bits,trigl0[i],trigl1[i],timeout
  end 
   ;
  ;_ interferometry of 2 perp antennas
  ;if timeout eq -1 then begin
   ;interferometry window position, at max
   ;val=0. & calwf=fltarr(nsam,2)
   ;for i=0,1 do begin
   ;  calwf[*,i]=wf[*,i]*cal[i]/1e6 ; calibrated to mT/s
   ;  vali=max(abs(calwf[*,i]),indi)
   ;  if vali gt val then ind=indi
   ;end
    ;ind=nsam-postsam ; at trigger

    ;Interferometry_2perp, calwf,intwf,fs, nsam, ind, ifus, ifint,$
    ;              sep,sepazim,azim,$
    ;             ind1,ind2,ifhp,ifintwf, iffs,ninter, finewf, fineintwf, fineif,$
    ;              scale2,lag2, maxcorr, sazim, sazim_,selev,ifver
    ;
    ;plot_Interferometry_2perp, file, cal,$
    ;              terr,t1,t2,poe,gaindB[gain],range,bits,trigl0,trigl1,timeout,$
    ;              calwf,intwf,fs, ind, ifus, ifint,$
    ;              trigtime, nsam,postsam,$
    ;              sep,sepazim,azim,$
    ;              ind1,ind2,ifhp,ifintwf, iffs,ninter, finewf, fineintwf, fineif,$
    ;              scale2,lag2, maxcorr, sazim, sazim_,selev,ifver
  ;end
 ;
  textoverview,file,wf,nant,intwf,cal,waitt,terr,timeout,scale2,lag2, maxcorr, sazim, sazim_,selev

end




pro rdth,file,trigtimes,waitt,nsam,postsam,terr,$
              t1,t2,gain,poe,fs,bits,range,trigl0,trigl1,timeout
  ;
  ; reads auxilliary txt file
  ;
  li=''
  times='2013-08-09T00:00:00.000.000.000'
  waitt='0ns' & nsam=33553920l & postsam=nsam
  terr='0ns' &  t1='0 C'&  t2='0 C' & gain=4b & poe='0V'
  fs=2e8 & bits=12 & range=[4.,4.,4.,4.]& trigl0=[0.,0.,0.,0.] & trigl1=[0.,0.,0.,0.]
  timeout=0
  range0=1 & trig00=1 & trig10=1 ; switching between channels
  ;
 txtfile=strmid(file,0,strlen(file)-4)+'.txt'
 if file_test(txtfile) then begin
  openr,1,txtfile
  while (not eof(1)) do begin
    readf,1,li
    if strmid(li,0,15) eq 'Est time error ' then $
        terr=strmid(li,15,strlen(li)-15)
    if strmid(li,0,7) eq 'Core 0:' then $
        t2=strtrim(strmid(li,7,strlen(li)-7),2)
    if strmid(li,0,7) eq 'temp2: ' then $
        t2=strtrim(strmid(li,7,strlen(li)-7),2)
    if strmid(li,0,7) eq 'temp1: ' then $
        t1=strtrim(strmid(li,7,strlen(li)-7),2)
    if strmid(li,0,12) eq 'gain set to '  then $
        reads,strmid(li,12,strlen(li)-12),gain
    if strmid(li,0,15) eq 'Preamp voltage:'  then $
        poe=strmid(li,15,strlen(li)-15)
    if strmid(li,0,12) eq 'PoE voltage:'  then $
        poe=strmid(li,12,strlen(li)-12)
    ;
    if strmid(li,0,11) eq 'Resolution:' then $
        reads,strmid(li,11,strlen(li)-11),bits
    if strmid(li,0,14) eq 'Sampling freq:'   then $
        reads,strmid(li,14,strlen(li)-14),fs
    aux=0.
    if strmid(li,0,7) eq ' Range:' then begin
        reads,strmid(li, 7,strlen(li)-7),aux
        if range0 then begin
          range0=0 & range[0]=aux
        end else range[1]=aux & range[2]=aux & range[3]=aux
    end
    if strmid(li,0,15) eq '  Trig 0 level:' then begin
        reads,strmid(li,15,strlen(li)-15),aux
        if trig00 then begin
          trig00=0 & trigl0[0]=aux
        end  else trigl0[1]=aux
    end
    if strmid(li,0,15) eq '  Trig 1 level:' then begin
        reads,strmid(li,15,strlen(li)-15),aux
        if trig10 then begin
          trig10=0 & trigl1[0]=aux
        end  else trigl1[1]=aux
    end
    ;
    if strmid(li,0,12) eq 'Trigger at: ' then $
        trigtimes=strmid(li,12,23)+'.'+strmid(li,35,3)+'.000'
    if strmid(li,0,11) eq 'Waited for ' then $
        waitt=strmid(li,11,strlen(li)-11)
    if strmid(li,0,13) eq 'Postsamples: '  then $
        reads,strmid(li,13,strlen(li)-13),postsam
    if strmid(li,0,15) eq 'Valid samples: '  then $
        reads,strmid(li,15,strlen(li)-15),nsam
    if strmid(li,0,18) eq 'SW trigger timeout' then timeout=1
  end
  close,1
 end else print,'  ***'+txtfile+' NOT FOUND ***'
end


pro plot_one_sevan,wfi, i, file, fs, fhp, times, nsam, postsam, bits,  trigl0i, trigl1i, timeout
    psfile=strmid(file,0,strlen(file)-4)+'_sp'+strtrim(string(i),2)+'.ps'
    plot_init,psfile,yfrac=1
    !p.charsize=1. & !p.multi=0 &  !p.thick=2 & !x.thick=3 & !y.thick=3


     trigs=['Trigger','Timeout']
     miny=min(wfi,max=maxy) & yr=[miny,maxy]+[-1.,1.]*(maxy-miny)*0.1
     plot_wf,wfi,times,fs=fs,fhp=fhp,pos=[0.09,0.06,0.995,0.93],/noerase,$
                yran=yr,ytit=' Raw data [TM]',tit='SEVAN CH'+strtrim(string(i-1),2)
               ;'!d '+trigs[timeout]+' at '+$
               ;strtrim(string(round((trigl0i*2d0-1d0)*2d0^(bits-1))),2)+' and '+$
               ;strtrim(string(round((trigl1i*2d0-1d0)*2d0^(bits-1))),2)+' units, '+$
               ;strtrim(string(round(double(nsam-postsam)/fs*1d3)),2)+' ms!n'
        xyouts,/norm,al=1.,0.606,0.086,col=255,$
        '!c!mm!3 = '+strtrim(string(total(double(wfi))/double(nsam),form='(f30.2)'),2)+$
        ', !ms!3 = '+strtrim(string(stddev(wfi,/double),form='(f30.2)'),2)+' units  '
        xyouts,/norm,al=1.,0.836,0.064,col=255,$
         strtrim(string(min(wfi),form='(i30)'),2)+' to '+$
         strtrim(string(max(wfi),form='(i30)'),2)+' units  '

     plot_close
end


pro plot_one_antenna, wfi, cali, intwfi, $  ; data, calibration in, integrated out
          i, azimi, file, fs, fhp,times,durus,$ ; ant no, azimuth, file,...
              trigtimes,waitt,nsam,postsam,terr,$ ; configuration to report...
              t1,t2,gdB,poe,bits,rangei,trigl0i,trigl1i,timeout
    ;
    ; ps plot and overview spectrogram for one antenna
    ; outputs integrated waveform
    ;
    psfile=strmid(file,0,strlen(file)-4)+'_sp'+strtrim(string(i),2)+'.ps'
    plot_init,psfile,yfrac=1
    !p.charsize=1. & !p.multi=0 &  !p.thick=2 & !x.thick=3 & !y.thick=3

    ;
    ;_title and global log spectrum
    trigs=['TRIGGER','SURVEY']
    plot_spect,/logf,fs=fs,wfi,pos=[0.09,0.6,0.995,0.93],$
               0,nsam-1, ytit='Power spectrum [dB from '+strtrim(string(round(cali/1d3)),2)+' !mm!3T/s]!c'+$
              ' !d 0dB...Noise with sigma at 1 dig. unit, Bandwidth = '+strtrim(string(round(fs/2d6)),2)+' MHz!n',$
               tit='SLAVIA '+ strtrim(string(i+1),2)+' (loop azimuth '+strtrim(string(round(azimi)),2)+'!uo!n)  '+$
                    strtrim(string(long(nsam/fs*1e3)),2)+' ms from '+times+' after '+waitt+' '+trigs[timeout]+'            !c!d'+$
               ' !d'+file+', dt='+terr+', T1='+t1+', T2='+t2+', U='+poe+$
               ', G='+strtrim(string(round(gdB)),2)+'dB, R='+$
               strtrim(string(round(rangei*1e3)),2)+'mV, N='+strtrim(string(nsam),2)+$
               ' at '+strtrim(string(fs/1e6,form='(f30.2)'),2)+'MHz, bits='+strtrim(string(bits),2)+$
               ', BLESK'+ver()+'!n',col=255,freq=freq,spec=spec , xran=[fhp,fs/2.]
    ;
    ;_integrated global waveform
    plot_wf,wfi*cali,times,fs=fs,fhp=fhp,out_wf=intwfi,$
               /int,out_yran=yri,ytit='Integrated data [nT]!c!d Highpass filtered at '+$
               strtrim(string(round(fhp/1d3)),2)+' kHz',$
               pos=[0.09,0.293,0.536,0.52],/noerase,/notimeaxis
       xyouts,/norm,al=1.,0.536,0.519,col=255,$
        '!c!ms!3 = '+strtrim(string(stddev(intwfi,/double),form='(f30.2)'),2)+' nT  '
       xyouts,/norm,al=1.,0.536,0.297,col=255,$
        strtrim(string(min(intwfi),form='(f30.2)'),2)+' to '+$
        strtrim(string(max(intwfi),form='(f30.2)'),2)+' nT  '
    ;
    ;_raw global waveform
     trigs=['Trigger','Timeout']
     miny=min(wfi,max=maxy) & yr=[miny,maxy]+[-1.,1.]*(maxy-miny)*0.1
     plot_wf,wfi,times,fs=fs,fhp=fhp,pos=[0.09,0.06,0.536,0.287],/noerase,$
                yran=yr,ytit=' Raw data ['+strtrim(string(round(cali/1d3)),2)+' !mm!3T/s]!c'+$
               '!d '+trigs[timeout]+' at '+$
               strtrim(string(round((trigl0i*2d0-1d0)*2d0^(bits-1))),2)+' and '+$
               strtrim(string(round((trigl1i*2d0-1d0)*2d0^(bits-1))),2)+' units, '+$
               strtrim(string(round(double(nsam-postsam)/fs*1d3)),2)+' ms!n'
        xyouts,/norm,al=1.,0.536,0.286,col=255,$
        '!c!mm!3 = '+strtrim(string(total(double(wfi))/double(nsam),form='(f30.2)'),2)+$
        ', !ms!3 = '+strtrim(string(stddev(wfi,/double),form='(f30.2)'),2)+' units  '
        xyouts,/norm,al=1.,0.536,0.064,col=255,$
         strtrim(string(min(wfi),form='(i30)'),2)+' to '+$
         strtrim(string(max(wfi),form='(i30)'),2)+' units  '
     ;
     ;_details around trigger
     plot_wf,reform(intwfi),trigtimes,/cent,$
                nsam-postsam-round(1d-6*durus*fs),nsam-postsam+round(1d-6*durus*fs), fs=fs,fhp=fhp,$
               yran=yri,pos=[0.54,0.293,0.995,0.52],/noerase,/notimeaxis,/noyaxis
     plot_wf,reform(wfi),trigtimes,/cent,yran=yr,$
                nsam-postsam-round(1d-6*durus*fs),nsam-postsam+round(1d-6*durus*fs),$
                fs=fs,fhp=fhp,pos=[0.54,0.06,0.995,0.287],/noerase,/noyaxis
     plot_close
     ;
     ; overview spectrogram (from blesk.pro)
     file_copy,psfile,'sp'+strtrim(string(i),2)+'.ps',/overwrite
     overviewspect,file, times,double(nsam)/fs/3600.,freq,spec, $
           savfile='overview'+strtrim(string(i),2)+'.idl', $
           psfile='overview'+strtrim(string(i),2)+'.ps',/logy,trig=(timeout eq 0)
     ;
     ;100 us details
      pal=long(fs*100e-6) ; parts of 100us
      thrcsm=7l ; threshold for max abs dev in units of std
      maxnsm=100; maximum number of smallest waveforms to be plotted
      psfile=strmid(file,0,strlen(file)-4)+'_sm'+strtrim(string(i),2)+'.ps'
      ;
      ;
      nsm=long(nsam/pal)
      dsm=reform(wfi[0: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)
      if nprtsm gt 0 then begin
       ;
       ;exclude frequencies from int
       ; RMC 216 kHz+-10 kHz
        fexca=216e3 & fexbw=20e3
        fexind=(fexca + [-fexbw/2., +fexbw/2.])*nsam/fs
        nsam2=2L^ceil(alog10(nsam)/alog10(2))
        addsam1=(nsam2-nsam)/2L & addsam2=nsam2-nsam-addsam1
        if addsam1 eq 0 then $
            if addsam2 eq 0 then ffti=fft(reform(intwfi),-1) $
                            else ffti=fft([reform(intwfi),fltarr(addsam2)],-1) $
                        else $
            if addsam2 eq 0 then ffti=fft([fltarr(addsam1),reform(intwfi)],-1) $
                            else ffti=fft([fltarr(addsam1),reform(intwfi),fltarr(addsam2)],-1)
        ffti[fexind[0]:fexind[1]]=complex(0.,0.)
        ffti[nsam-fexind[1]:nsam-fexind[0]]=complex(0.,0.)
        intwfi_=(float(fft(ffti,1)))[addsam1:nsam+addsam1-1]
      ;
        plot_init,psfile,yfrac=1
        !p.charsize=1. & !P.MULTI = [0, 2, 2, 0, 0] &  !p.thick=2 & !x.thick=3 & !y.thick=3
        for iprtsm =0, (maxnsm < (nprtsm-1) ) do begin
          plot_wf,reform(wfi),times,fs=fs,tunit=1e-3, partsm[iprtsm]*pal, (partsm[iprtsm]+1)*pal-1,$
               ytit='Part '+strtrim(string(iprtsm+1),2)+'/'+$
               strtrim(string(nprtsm),2)+' |wf-'+strtrim(string(round(offss)),2)+'|>'+$
               strtrim(string(round(thrsm)),2)+' units '+$
               '!c !c  Raw data ['+strtrim(string(round(cali/1d3)),2)+' !mm!3T/s]'
          plot_wf,reform(intwfi_),times,fs=fs,tunit=1e-3, partsm[iprtsm]*pal, (partsm[iprtsm]+1)*pal-1,$
               ytit=strtrim(string(round(fexbw/1e3)),2)+' kHz around '+strtrim(string(round(fexca/1e3)),2)+$
               ' kHz excluded!c!c Integrated data [nT]'
        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 ----------------------'
         ;
         plot_close
      end
end

pro textoverview,file,wf,nant,intwf,cal,waitt,terr,timeout,scale2,lag2, maxcorr, sazim, sazim_,selev
       ;
    ; text overview
    ;for i=0,nant-1 do begin
      minraw=[min(wf[*,0]),min(wf[*,1])]
      maxraw=[max(wf[*,0]),max(wf[*,1])]
      sigraw=[stddev(wf[*,0],/double),stddev(wf[*,1],/double)]
      minint=[min(intwf[*,0]),min(intwf[*,1])]
      maxint=[max(intwf[*,0]),max(intwf[*,1])]
      sigint=[stddev(intwf[*,0],/double),stddev(intwf[*,1],/double)]
    ;en
    txt='overview.txt'
    if file_test(txt) then openu,1,txt,/append else begin
      openw,1,txt
      printf,1,'File             Min1  Min2  Max1  Max2  Sigma1  Sigma2'+$
		' Mn1[nT] Mn2[nT]'+$
                ' Mx1[nT] Mx2[nT]'+$
		' Sg1[nT] Sg2[nT] Cal1 Cal2'
               ;Wait  TimeErr'; TempBox TempCPU'+$
               ;' Mx1[nT] Mx2[nT] Sg1[nT] Sg2[nT] Cal1 Cal2 Trig1 Trig2 Wait'  TimeErr TempBox TempCPU'+$
               ;'   Voltage dB Nsamples MHz bit'
    end
    pos0=(strpos(file,'/',/reverse_s)>strpos(file,'\',/reverse_s))+1
    ;if timeout then interfs='' else interfs=$
    ;  string(scale2,lag2, maxcorr, sazim, sazim_,selev,form='(f8.2,f7.1,f5.2,3i5)')
    interfs=''
    printf,1,strmid(file,pos0,15)+$
      string(minraw,maxraw,sigraw,minint,maxint,sigint,form='(4i6,10f8.1)')+$
      string(round(cal/1d3),form='(2i5)')+interfs
    ;string(round(trigl0*2d0^(bits-1)),form='(2i6)')+$
    ;' '+waitt+' '+terr ;string(waitt,terr,t1,t2,poe,gaindB[gain],nsam,$
    ;       round(fs/1e6),bits ,form='(a5,a9,2a8,a10,i3,i9,i4,i4)' )
    close,1
end

pro Interferometry_2perp, wf,intwf,fs,  nsam,ind,ifus, ifint,$
                  sep,sepazim,azim,$
                  ind1,ind2,ifhp,ifintwf, iffs, ninter,finewf, fineintwf, fineif,$
                  scale2,lag2, maxcorr, sazim, sazim_,selev,ifver
   ;
   ; interferometry for 2 perp antennas
   ;
   ; wf,intwf,fs, ... input wf, integrated wf, sampling fr,
   ; nsam,ind, ifus,   ... input total number of samples, central index of the interferometry interval, its duration in us,
   ; ifint        ... switch between original or integrated wf to be used in iterferometry
   ; sep,sepazim,azim         ... input antenna separation [m], separation azimuth [deg], azim. of antenna planes arr(2) [deg]
   ; ind1,ind2,ifhp        ... output indices from-to, inverse boxcar highpass filter frequency
   ; ifintwf, iffs, ninter,finewf, fineintwf, fineif ... output filtered integ wf,
   ;                                               sampling fr and no of interpolated points in the if interv.
   ;                                               the if interv: wf, integ wf, selected one
   ; scale2,lag2, maxcorr               ... output scale of ant2, its lag, Pearson correlation coeff
   ; sazim, sazim_,selev                ... source azumith from scale, from lag, elevation (3333. fill value),
   ;
   ifver='IF20130903'
   ;
   ifhp=2e6/ifus ; two highpass inverse boxcar filter intervals per window
   ;               for the integrated signal
   iffs=10.*fs ; interferometry time resolution: 0.1  sampling interval
   ;
   ; no of orig samples in the interferometry interval, its boundaries
   nind=round(1d-6*ifus*fs)
   ind1=ind-nind/2 & ind2=ind+nind/2-1
   if ind1 lt 0 then begin
     ind1=0 & ind2= nind-1
   end
   if ind2 gt nsam-2 then begin
     ind2=nsam-2 & ind1= nsam-nind-1
   end
   ;
   ; highpass filter the integrated signal
   ; interpolate both to iffs sampling frequency
   ninter=round((nind-1)*iffs/fs)+1
   ifintwf=intwf & finewf=fltarr(ninter,2) & fineintwf=fltarr(ninter,2)
   for i=0,1 do begin
     ; inverse boxcar highpass filter
     ifintwf[*,i]=intwf[*,i]-smooth(intwf[*,i],round(fs/ifhp),/nan,/EDGE_TRUNCATE)
     ; interpolation to iffs
     finewf[*,i]=interpolate(reform(wf[ind1:ind2+1,i]),cubic=-0.5,0.5+findgen(ninter)/float(ninter-1)*float(nind-1))
     fineintwf[*,i]=interpolate(reform(ifintwf[ind1:ind2+1,i]),cubic=-0.5,0.5+findgen(ninter)/float(ninter-1)*float(nind-1))
   end;
   ;
   ; select integrated or derivative
   if ifint then fineif=fineintwf else fineif=finewf
   ;
   ; val to scale antennas, at the central window position
   vali=fltarr(2)
   for i=0,1 do begin
     vali[i]=max(abs(fineif[*,i]))
   end
   ;
   ; optimum lag and scale , in iffs samples
   maxlag=round(sep/3e8*iffs)  ;  in iffs samples,  sep[m] /c [m/s]
   maxcorr=0.
   for lag=-maxlag,maxlag do begin
     corr=(C_CORRELATE( fineif[*,0], vali[0]/vali[1]*fineif[*,1],lag,/double))[0]
     if abs(corr) gt abs(maxcorr) then begin
       maxcorr=corr
       lag2=lag
     end
    end
   scale2=vali[0]/vali[1]*maxcorr/abs(maxcorr)
   maxcorr=abs(maxcorr)
   ;
   ;_source azimuth from scale
   sazim=azim[1]+atan(scale2)/!dtor
   if cos((sepazim-sazim)*!dtor)*lag2 gt 0 then sazim=sazim+180.
   if sazim gt 180. then sazim=sazim-360.
   ;
   ;_source elevation and alternative source azimuth from lag
   if abs(lag2) gt abs(cos((sepazim-sazim)*!dtor)*sep/3e8*iffs) then begin
      selev=0.
      if abs(lag2) gt sep/3e8*iffs then  if lag2 lt 0 then sazim_b=0. $
                                                      else sazim_b=180. $
                                   else sazim_b=acos(-lag2*3e8/iffs/sep)/!dtor
      if tan((azim[0]-sazim_b-sepazim)*!dtor)*scale2 lt 0 then sazim_b=-sazim_b
      sazim_=sepazim+sazim_b
      if sazim_ gt 180. then sazim_=sazim_-360.
      if sazim_ lt -180. then sazim_=sazim_+360.
    end else begin
      if abs(cos((sepazim-sazim)*!dtor)) lt 0.1 then selev=3333. else $
       selev=acos( abs(lag2 / cos((sepazim-sazim)*!dtor)/sep*3e8/iffs) ) /!dtor
       sazim_=3333.
    end
end


pro plot_Interferometry_2perp, file, cal,$
                  terr,t1,t2,poe,gdB,range,bits,trigl0,trigl1,timeout,$
                  wf,intwf,fs, ind, ifus, ifint,$
                  trigtime, nsam,postsam,$
                  sep,sepazim,azim,$
                  ind1,ind2,ifhp,ifintwf, iffs, ninter, finewf, fineintwf, fineif,$
                  scale2,lag2, maxcorr, sazim, sazim_, selev,ifver
;
; antenna colors (black, red)
 ifcol=[0,250]
 tns_str, trigtime-double(nsam-postsam+ind)/fs*1d9, iftimes ;; central time of the if window
;
 psfile=strmid(file,0,strlen(file)-4)+'_if.ps'
    plot_init,psfile,yfrac=0.5
   !p.charsize=0.5 & !p.multi=0 &  !p.thick=2 & !x.thick=3 & !y.thick=3
   ;
   ; extremes in the windows
   minifwf=min(finewf,max=maxifwf)
   minifint=min(fineintwf,max=maxifint)
   fineif[*,1]=fineif[*,1]*scale2
   minif=min(fineif,max=maxif)
   ;
   if sazim_ eq 3333. then sazim_s='' else sazim_s=' or (from lag) '+strtrim(string(round(sazim_)),2)+'!uo!n'
   if selev eq 3333. then selevs='N/A' else selevs=strtrim(string(round(selev)),2)+'!uo!n'
   for i=0,1 do begin
     if i eq 0 then fitwf=fineif[*,0] $
               else if lag2 lt 0 then fitwf=[replicate(!values.f_nan,-lag2), fineif[0:ninter-1+lag2,1] ] $
               else if lag2 gt 0 then fitwf=[fineif[lag2:ninter-1,1],replicate(!values.f_nan,lag2) ] $
               else fitwf=fineif[*,1]
     tity=['mT/s','nT'] ; [  der or int]
     plot_wf,fitwf,iftimes,/cent,yran=[minif,maxif],$
                tunits=1e-9,/notimeaxis,$
                fs=iffs,pos=[0.05,0.526,0.995,0.93],overplotcolor=ifcol[i],$
                ytit='Interferometry data ['+tity[ifint]+']!c!dInterpolated to '+$
                        strtrim(string(1./iffs*1e9,form='(f30.1)'),2)+'-ns resolution!n',$
               tit='INTERFEROMETRY: SLAVIA 1 (black) - SLAVIA 2 (red, lagged and scaled), Baseline = '+$
                      strtrim(string(sep,form='(f30.1)'),2)+$
                     'm, Azim!d1-2!n = '+strtrim(string(round(sepazim)),2)+'!uo!n, Scale!d2!n = '+$
                      strtrim(string(scale2,form='(f30.2)'),2)+', Lag!d2!n = '+$
                      strtrim(string(lag2/iffs*1e9,form='(f30.1)'),2)+' ns, '+$
                      ' Corr = '+strtrim(string(maxcorr,form='(f30.2)'),2)+'  Source azim = '+$
                      strtrim(string(round(sazim)),2)+'!uo!n'+sazim_s+',  elev = '+selevs+'!c'+$
               ' !d SLAVIA 1 loop plane azimuth = '+strtrim(string(round(azim[0])),2)+' degrees,  '+$
               'SLAVIA 2 loop plane azimuth = '+strtrim(string(round(azim[1])),2)+' degrees,  File='+$
               file+', Time accuracy = '+terr+',  T1 = '+t1+',  T2 = '+t2+',  Voltage = '+poe+$
               ',  Gain='+strtrim(string(round(gdB)),2)+'dB,  Range1='+$
               strtrim(string(round(range[0]*1e3)),2)+'mV,  Range2='+$
               strtrim(string(round(range[1]*1e3)),2)+'mV,  Nsamp='+strtrim(string(nsam),2)+$
               ' at '+strtrim(string(fs/1e6,form='(f30.2)'),2)+'MHz,  Bits='+strtrim(string(bits),2)+$
				', BLESK'+ver()+ifver+'!n'
   end
   for i=0,1 do begin
     plot_wf,reform(fineintwf[*,i]),iftimes,$
               /cent,yran=[minifint,maxifint],/notimeaxis,tunits=1e-9,/noerase,$
               fs=iffs,$
               pos=[0.05,0.293,0.995,0.52],overplotcolor=ifcol[i],$
               ytit='Integrated data [nT]!c!d Highpass filtered at '+$
                     strtrim(string(round(ifhp/1d3)),2)+' kHz'
     plot_wf,reform(ifintwf[ind1:ind2,i]),iftimes,tunits=1e-9,$
               /cent, fs=fs,overplotcolor=ifcol[i],overplotsym=0.3
   end
   ;
   trigs=['Trigger','Timeout']
   for i=0,1 do begin
     plot_wf,reform(finewf[*,i]),iftimes+'!d '+trigs[timeout]+' at '+$
               strtrim(string(round(double(nsam-postsam-ind)/fs*1d3)),2)+' ms; '+$
               'SLAVIA1: outside '+strtrim(string( cal[0]/1e6*(trigl1[0]*2d0-1d0)*2d0^(bits-1),form='(f20.2)'),2)+$
                   ' and '+strtrim(string( cal[0]/1e6*(trigl0[0]*2d0-1d0)*2d0^(bits-1),form='(f20.2)'),2)+' mT/s;  '+$
               'SLAVIA2: outside '+strtrim(string( cal[1]/1e6*(trigl1[1]*2d0-1d0)*2d0^(bits-1),form='(f20.2)'),2)+$
                   ' and '+strtrim(string( cal[1]/1e6*(trigl0[1]*2d0-1d0)*2d0^(bits-1),form='(f20.2)'),2)+' mT/s',$
                /cent,yran=[minifwf,maxifwf], tunits=1e-9,$
                fs=iffs,pos=[0.05,0.06,0.995,0.287],overplotcolor=ifcol[i],/noerase,$
                ytit=' Calibrated data [mT/s]!c'+$
                '!d1:'+strtrim(string(round(cal[0]/1d3)),2)+'!mm!3T/s/unit; '+$
                  '2:'+strtrim(string(round(cal[1]/1d3)),2)+'!mm!3T/s/unit'
     plot_wf,reform(wf[ind1:ind2,i]),iftimes,tunits=1e-9,$
               /cent, fs=fs,overplotcolor=ifcol[i],overplotsym=0.3
   end
   plot_close
     file_copy,psfile,'if.ps',/overwrite
end
