program winselect
c       USE DFLIB
c  ____________________________________________________________________
c |    VERSION PYB - LISBON/GRENOBLE, 13/07/2004                                |
c | This subroutine selects, on 3D noise recordings,                   |
c |         the most stationary windows on the basis of an             |
c |        "antitrigger" algorithm (sta/lta ratio                      |
c |         should remain between given thresholds)                    | 
c |     with two additional, optional, criteria                         |
c |          - non-saturation (lsat)                                   |
c |          - relatively quiet windows (lnoisy)                       |
c |      moderated by a "tolerance" criterion: the non-fulfillment     |
c |      of all these criteria is "tolerated" as long as the total     |
c |      associated duration does not exceed a prescribed length rtol  |
c |____________________________________________________________________|
c
c
c   Input : Option (selection parameters), signal (3 components)
c   Output : number of windows nwin + begin and end times (lbeg, lend)
c
      implicit none
c      
	integer ntime, nwinmax
      parameter(ntime=600001)
      parameter(nwinmax=1000)
c
      real rup(ntime),rns(ntime),rew(ntime)
      real s1(ntime),s2(ntime),s3(ntime),s(ntime),ssat(ntime)
      real sta(ntime),lta(ntime),ratio(ntime)
	real wi(ntime)
	integer lsumcrit(ntime),nlongtol
      integer lbeg(nwinmax),lend(nwinmax)
      integer ndata, ntlong
      integer lnoisy,lsat,lvalinit,noverlap,nstep
      integer lout,lk,lnoise,nt,nsta,nlta,ncos,nzero,nexpl
      integer ncz,ndatacz
      integer ind,ind1,ind2,ind3,ind4,ind41,ind42
      integer lfen,nfen,i0,ideb,lk3
      integer i,j,it,i1,ista,ilta,iwin
      integer longcar
      real  deltat,tstep,tlong,xtlong
      real ltamin,ltamax,ltamoy,ltamean,roverlap,toverlap,rtol
	real rvalinit,tsta,tlta,seuilmin,seuilmax
	real smoy,smoy1,smoy2,smoy3,ss,ss1,ss2,ss3,smin,smax
	real coefew,coefns,ratmax,ratmin,tlongtol,tbeg,tend
	real alog10,sqrt
      character*200 noisefilename, noisefilenam2
      character kup*2,kns*2,kew*2,csuf*3
      character chz*3,chn*3,che*3
      character station*10
c
        character*200 fnlst,fnpar,fnout
        integer ier, iext
        real datath(ntime,3)
        equivalence(datath(1,1),rup(1))
        equivalence(datath(1,2),rns(1))
        equivalence(datath(1,3),rew(1))
ccccccc
c
c     standard output
      lout = 6
c_______________________________________________________________________
c    Initialisation: files to be read or written
        call  gcmlnw(fnlst,fnpar,fnout,ier)
        if ( ier  .eq.  1 ) stop
        open(10, file=fnlst,status='old')
        open(20, file=fnpar,status='old')
        open(30, file=fnout,status='unknown')
c_______________________________________________________________________
c    Initialisation: Read window selection parameters
c
c links with GUI to enter parameters
c    'tsta = length of the window on which the SHORT term
c            average is computed (standard value: 1.0 s)'
c    'tlta = length of the window on which the LONG term',
c             average is computed (standard value: 25 s)'
c
c  Thresholds
c      'sta/lta thresholds: minimum value = seuilmin ',
c      'sta/lta thresholds: maximum value = seuilmax ',
c       to be fulfilled on a minimum length tlong (at least 10 times the
c        longest period of interest, 20 times is better)
c
c  Other criteria
c      'No saturation condition : signal < 99.5% (Max-signal)',
c      'Quiet window condition : lta < median lta value',
c
c  Tolerance
c      'Criteria may not be respected for a limited duration rtol',
c
c   Overlap percentage over two next windows: roverlap (from 0.0 to 0.5)
c
       rvalinit=1.005e+20
       tsta=rvalinit
       tlta=rvalinit
       seuilmin=rvalinit
       seuilmax=rvalinit
       tlong=rvalinit
       roverlap=rvalinit
 	 rtol=rvalinit
       lvalinit=1000
       lsat=lvalinit
       lnoisy=lvalinit
       read(20,*) tsta, tlta, seuilmin, seuilmax, tlong, roverlap, lsat,
     &            lnoisy,rtol
       
       write(6,*)'I have read input_selectauto'
c
c      Default values
c
       if(tsta.eq.rvalinit) tsta=1.0
       if(tlta.eq.rvalinit) tlta=25.
       if(seuilmin.eq.rvalinit) seuilmin=0.5
       if(seuilmax.eq.rvalinit) seuilmax=2.0
       if(tlong.eq.rvalinit) tlong=20.
       if(roverlap.eq.rvalinit) roverlap=0.2
       if(lsat.eq.lvalinit) lsat=1
       if(lnoisy.eq.lvalinit) lnoisy=0
	 if(rtol.eq.rvalinit) rtol=0.
ccc    
ccc      tsta=1.0
ccc      tlta=25.0
ccc      seuilmin=0.5
ccc      seuilmax=2.0
ccc      tlong=10.0
ccc      roverlap=0.0
ccc      lsat=1
ccc      lnoisy=0
ccc      rtol=0.
c
c ______________________________________________________________________
c    Initialisation: Read noise files 
c      Authorized formats: gse and saf
c
c     Needed : rup, rns, rew, ndata, deltat
c     write(6,*) 'Name of the noise file (format a40)'
100   continue
      read(10,1000,end=80) noisefilename
      write(6,*) 'File to read: ',noisefilename
      lk = longcar(noisefilename)
      lk3=lk-2
      csuf=noisefilename(lk3:lk)
      coefew=1.
      coefns=1.
c     
      lnoise=0
        iext = 0
      if(csuf.eq.'gse' .or. csuf.eq.'GSE') iext=1
      if(csuf.eq.'saf' .or. csuf.eq.'SAF' .or. csuf.eq. '_fi') iext=2
        call  readerthw(iext,noisefilename,datath,ntime,
     &  ndata,deltat,chz,chn,che,station,ier)
c                        write(*,*)' ndata ier  ', ndata, ier
        kup = 'ok'
        kns = 'ok'
        kew = 'ok'
       lnoise=1
      if(ier.eq.1) then
        write(6,*)' Following file can not be read:'
        write(6,'(1x,a200)')noisefilename
        close(10)
        close(20)
        close(30)
        stop 200
       endif
c
c  Still to be done/ decided:
c   With or without instrument correction: here no instrument correction
c   are applied
c
c ______________________________________________________________________
      xtlong=tlong/deltat
      nexpl=alog10(xtlong)/alog10(2.)
      nt=2.**(nexpl+1)
      if(nt.gt.ntime) then
        write(lout,*) 'Problem : tlong too large : nt larger than ntime'
        write(lout,*) 'tlong =', tlong, ' nt= ', nt,'  ntime =',ntime
        stop 210
                      endif
c
      nlta=tlta/deltat
      nsta=tsta/deltat
c
c___________________________________________________________________
c
c   Computing the sta and lta values all along the raw signal
c
c  a) Tapering with cosine window, 0.5 s long on each side
c
      ncos=0.5/deltat
      nzero=2
c
      call iniwin(wi,ndata,nzero,ncos)
      do it=1,ndata
          s2(it)=rns(it)
          s3(it)=rew(it)
          s1(it)=rup(it)
          enddo
c
      call window(s1,s,1,ndata,wi)
      smoy=0.
      do i=1,ndata
        smoy=smoy+s(i)
        enddo
      smoy1=smoy/float(ndata)
      do i=1,ndata
        s1(i)=s(i)-smoy1
        enddo
c
      if(kns.eq.'ok') then
      call window(s2,s,1,ndata,wi)
      smoy=0.
      do i=1,ndata
        smoy=smoy+s(i)
        enddo
      smoy2=smoy/float(ndata)
      do i=1,ndata
        s2(i)=s(i)-smoy2
        enddo
                      else
      do i=1,ndata
        s2(i)=0.
        enddo
                      endif
c
      if(kew.eq.'ok') then
      call window(s3,s,1,ndata,wi)
      smoy=0.
      do i=1,ndata
        smoy=smoy+s(i)
        enddo
      smoy3=smoy/float(ndata)
      do i=1,ndata
        s3(i)=s(i)-smoy3
        enddo
                      else
      do i=1,ndata
        s3(i)=0.
        enddo
                      endif

c
c  b) computing the root mean square ground motion
c

      smoy=0.
      smax=0.
      smin=1.e+20
      ncz=ncos+nzero
      ndatacz=ndata + 1 - ncz
      write(6,*)'NCZ: ',ncz,' NDATACZ: ',ndatacz
      do i=1,ndata
        ss1=abs(s1(i))
        ss2=coefns*abs(s2(i))
        ss3=coefew*abs(s3(i))
        ss=ss1**2 + ss2**2 + ss3**2
        s(i)=sqrt(ss)
        smoy=smoy+s(i)
        ssat(i)=ss1
        if(ss2.gt.ssat(i)) ssat(i)=ss2
        if(ss3.gt.ssat(i)) ssat(i)=ss3
        if(ssat(i).ge.smax) smax=ssat(i)
        if(ssat(i).eq.0.) ssat(i)=1.e-20
        enddo
      smoy=smoy/float(ndata)
      smax=0.995*smax
c      write(lout,1002) smoy,smax
c
c  c)  Computing the short term average STA :
c       as long as i < nsta + ncos + nzero, s(i-nsta) is replaced by  smoy
c
      do i=1,ncz
        sta(i)= smoy*float(nsta)
        enddo
      do i=ncz+1, ncz + nsta - 1
        i1=i-1
        sta(i)=sta(i1) + s(i) - smoy
        enddo
      do i=ncz+nsta,ndata-ncz
        i1=i-1
        ista=i-nsta
        sta(i)=sta(i1) + s(i) - s(ista)
        enddo
      do i=ndatacz,ndata
        i1=i-1
        ista=i-nsta
        sta(i)=sta(i1) + smoy - s(ista)
        enddo
c
c  d)  Computing the long term average LTA :
c       as long as i < nlta + ncos + nzero, s(i-nlta) is replaced by  smoy
c
c
      do i=1,ncz
        lta(i)= smoy*float(nlta)
        enddo
      do i=ncz+1,ncz + nlta - 1
        i1=i-1
        lta(i)=lta(i1) + s(i) - smoy
        enddo
      do i=ncz+nlta,ndata-ncz
        i1=i-1
        ilta=i-nlta
        lta(i)=lta(i1) + s(i) - s(ilta)
        enddo
      do i=ndatacz,ndata
        i1=i-1
        ilta=i-nlta
        lta(i)=lta(i1) + smoy - s(ilta)
        enddo
c
c   e) Computing the STA/LTA ratio
c
      ratmin=100.
      ratmax=0.
      ltamin=1.e+20
      ltamax=0.
	ltamean=0.
      do i=1,ndata
        sta(i)=sta(i)/float(nsta)
        lta(i)=lta(i)/float(nlta)
 	  ltamean=ltamean+lta(i)
        ratio(i)=sta(i)/lta(i)
        if(ratmin.gt.ratio(i)) ratmin=ratio(i)
        if(ratmax.lt.ratio(i)) ratmax=ratio(i)
        if(ltamin.gt.lta(i)) ltamin=lta(i)
        if(ltamax.lt.lta(i)) ltamax=lta(i)
        enddo
      ltamean = ltamean/float(ndata)
      ltamoy=0.8*ltamax + 0.2*ltamin

c
c   f) Identification of working windows
c
c New version
c  Computes at each time sample i the value of sumcrit 
c           over the previous "tlong" (seconds) (ntlong samples)
c
      lsumcrit(1)=0
      lfen=0
      ntlong=nint(tlong/deltat)
      do i=2,ndata
        ind1=ratio(i)/seuilmin
        ind2=seuilmax/ratio(i)
        ind3=smax/ssat(i)
        ind41=ltamoy/lta(i)
	  ind42=ltamean/lta(i)
	  ind4=ind41+ind42
        ind=ind1*ind2
        if(lsat.eq.1) ind=ind*ind3
        if(lnoisy.eq.1) ind=ind*ind4
        if(ind.ge.1) lsumcrit(i)= lsumcrit(i-1) + 1
        if(ind.eq.0) lsumcrit(i)= lsumcrit(i-1)
        enddo
c
c
      do j=ntlong+1,ndata	  
	   i = ndata+ntlong+1-j
	   i1=i-ntlong
         lsumcrit(i)=lsumcrit(i) - lsumcrit(i1)
        enddo
c
	tlongtol=0.999*(tlong-rtol)
	toverlap=roverlap*tlong
	tstep=tlong-toverlap
	nstep=tstep/deltat - 1
	ideb=ntlong + 1
	nlongtol = ntlong - nint(rtol/deltat)
c
	lfen=0
      write(6,*) ideb,tlong,rtol,tstep,nstep
	do 200 iwin=1,nwinmax
C        write(6,*) iwin,ideb,lfen
        i0=ideb
	  if(i0.gt.ndata) go to 210
	  do 201 i =i0, ndata
	    if(lsumcrit(i).ge.nlongtol) then
            lbeg(iwin)=i-ntlong
            lend(iwin)=i
	      ideb=i+nstep
	      lfen=iwin
	      go to 202
	                               else
            ideb=i+1
                                     endif
201       continue
c
202     continue
c        tbeg=float(lbeg(iwin)-1)*deltat
c        tend=float(lend(iwin)-1)*deltat
c       write(6,*) iwin,tbeg,tend
        if(ideb.gt.ndata) go to 210
c
200     continue
c
210   write(40,*) lfen
c
c
c  g) Saving the results
c
      write(6,*) lfen
c        write(30,*) lfen
      lk = longcar(noisefilename)
        noisefilenam2=' '
        noisefilenam2(1:1)='"'
        noisefilenam2(2:lk+1)=noisefilename
        noisefilenam2(lk+2:lk+2)='"'
      do i=1,lfen
        tbeg=float(lbeg(i)-1)*deltat
        tend=float(lend(i)-1)*deltat
c        write(6,*) i,tbeg,tend
        write(30,3000) noisefilenam2,tbeg,tend,iext,chz,chn,che,station
        enddo
c
c________________________________________________________________
c
      go to 100
c
80    continue
      close(10)
      close(20)
      close(30)
c
3000  format(a200,3x,f8.3,3x,f8.3,1x,i1,3(3x,a3),3x,a10)
1000  format(a200)
c
c1050  format(i10,3x,f6.0,2f11.3)
c1001  format(i2,3x,3f10.2)
c1002  format('   Max:',f12.5,10x,'Moyenne:',f12.5)
c1003  format('   Moyenne des signaux:',f12.5)
c1005  format(a40)
c
      stop
      end
c *********************************************************
c
      function longcar(mot)
      character*(*) mot
      do 1 l=len(mot),1,-1
        if(mot(l:l).ne.' ') go to 2
1       continue
2     longcar=l
      return
      end

c *********************************************************
      subroutine liskonno(nx,y,dx,bexp)
c ** smoothing of a function y (equally-spaced, dx) with the "Konno-Ohmachi"
c ** function sin (alog10(f/fc)^exp) / alog10(f/fc)^exp) ^^4 
c **   where fc is the frequency around which the smoothing is performed
c ** exp determines the exponent 10^(1/exp) is the half-width of the peak
c 
c ** cf Konno & Ohmachi, 1998, BSSA 88-1, pp. 228-241

      dimension y(*)
      common/lis0/ylis(16385),freq(16385)
c *********************************************
      if(nx.gt.16385) then
         write(6,*) ' nx =',nx,' too large for lisprop (max = 16385)'
         stop
         end if
c *********************************************
      do ix=1,nx
         freq(ix)=float(ix-1)*dx
        enddo
        fratio=10.**(2.5/bexp)
      ylis(1)=y(1)
      do 1 ix=2,nx
         fc=freq(ix)
         fc1=fc/fratio
         fc2=fc*fratio
         ix1=fc1/dx
         ix2=fc2/dx +1
         if(ix1.le.1) ix1=2
         if(ix2.ge.nx) ix2=nx
         a1=0.
         a2=0.
         do 2 j=ix1,ix2
            if(j.ne.ix) then
              c1=bexp*alog10(freq(j)/fc)
              c1=(sin(c1)/c1)**4
              a2=a2+c1
              a1=a1+c1*y(j)
                        else
              a2=a2+1.
              a1=a1+y(ix)
                        endif
2          continue
         ylis(ix)=a1/a2
    1    continue
      do 10 ix=1,nx
10       y(ix)=ylis(ix)
      return
      end
c ********************************************************
      subroutine window(datin,datout,n1,nwin,w)
c ********************************************************
c     datin    = initial data
c     datout   = windowed data
c     n1       = index of the first point
c     nwin     = number of points of the window
c     w(>nwin) = apodisation function
c ********************************************************
      real datin(*),datout(*),w(*)
c
      xmean=0.
      j=n1
      do 10 i=1,nwin
         datout(i)=datin(j)
         xmean=xmean + datout(i)
10       j=j+1
      xmean=xmean/float(nwin)
      do 20 i=1,nwin
20       datout(i)=datout(i)-xmean
      xmean=0.
      do 30 i=1,nwin
         datout(i)=datout(i)*w(i)

30       xmean=xmean+datout(i)
      xmean=xmean/float(nwin)
      do 40 i=1,nwin
40       datout(i)=datout(i)-xmean
      return
      end
c ********************************************************
      subroutine iniwin(wi,nwin,nzero,ncos)
c ********************************************************
c     creates the window shape, and stores in wi(>nwin)
c ********************************************************
      real wi(nwin),pi
      pi=3.1415927
c ********************************************************
      if(nwin.eq.0) then
         write(6,*) ' nwin = 0 ; iniwin stop 1'
         stop
         end if
c ********************************************************
      if(2*nzero.ge.nwin) then
         nzero=(nwin-1)/2
         write(6,*) ' nzero too large in iniwin, reset to:',nzero
         end if
c ********************************************************
      if(2*(nzero+ncos).ge.nwin) then

         ncos=(nwin-2*nzero-1)/2
         write(6,*) ' ncos (iniwin) reset to:',ncos
         end if
c ********************************************************
      do 10 i=1,nwin/2
         if(i.le.nzero) then
            wi(i)=0.
         else if(i.le.(nzero+ncos)) then
               wi(i)=.5*(1. - cos(pi*(i-nzero)/(ncos+1)))
              else
               wi(i)=1.
              end if
         wi(nwin+1-i)=wi(i)
10       continue
      return
      end