c-------------------------------------------------------------------- 
c         Find Omega_bar from tables given in Stat.dat
c-------------------------------------------------------------------- 
      PARAMETER (Nmax = 5000) ! length of table 
      PARAMETER (Nave = 5   ) ! number of snapshots to average periods 
      PARAMETER (aJump =-20.   ) ! jump in angle, which needs addition of 180 degrees
      COMMON/DATA/ Time(Nmax),Angle(Nmax),Period(Nmax),A2(Nmax),
     &       Aexpn(Nmax),Nstep(Nmax)   
      COMMON/WEI/Weights(0:100),iFlag
c      DATA Time,Angle,Period,A2/Nmax*0.,Nmax*0.,Nmax*0.,Nmax*0./

      Open(1,file='Stat2Clean.dat')
      Open(2,file='omegap2.dat')
      Nn       =1
      Nstep0   =0
      Nleft    =(Nave-1)/2
      Nright   = Nleft +Nave
      iFlag    = 0
 10   Read(1,*,err=20,end=20)Aexpn(Nn),Time(Nn),Nstep(Nn),
     &                      aL,Angle(Nn),Rad,A2(Nn)
c      write (*,*)  Aexpn(Nn),Time(Nn),Nstep(Nn)
           Nn = Nn +1
           if(Nstep(Nn).lt.Nstep0)Then
              write (*,*)  ' wrong line in Stat.dat:',Nn,Nstep(Nn) 
              Stop
           EndIf 
           Nstep0 = Nstep(Nn)
           goto 10
 20   Nn = Nn -1
      write (*,*)  ' Number of snapshots =',Nn
c      itest =70
c      do i=1,10
c      write (*,'(i4,6g11.4)')itest,Aexpn(itest),Time(itest),Angle(itest)
c      EndDo 

      Do i=1,Nn-2            ! find periods
         ang0 = Angle(i)        ! starting Angle
         per  = 0.              ! period
         ang1 = ang0            !  angle on previous snapshot
          da  = 0.              ! offset for angle
         If(i.eq.iTest)write (*,'(" i=",i5," Angle0=",f8.4)')i,ang0
         Do j=i+1,Nn
            ang = Angle(j)
            if(Angle(j)-Angle(j-1).lt.aJump)da =da +180.
             ang=ang+da
c       If(i.eq.iTest)write(*,'(5x,i4,5f9.4)')j,ang,da,Angle(j),ang1,ang0
            If(ang-ang0.gt.180.)Then   ! found period
               per = Time(j-1)-Time(i)+
     &               (ang0+180.-ang1)/(ang-ang1)*(Time(j)-Time(j-1))
               goto 30
            EndIf 
            ang1 = ang
         EndDo 
 30      Period(i) =per
c         If(i.eq.iTest)Then
c            write (*,*)  ' Period=',per
c            stop
c         EndIf 
      EndDo 
      iSmooth =(Nave+1)/2
      width = iSmooth/2.
      write (2,'(" Naverage=",i3," tolerance angle=",f8.2)')Nave,aJump 
      write (2,'("  Aexpn   T/Gyrs Step <Period>/Gyr Omega*Gyr",
     &           " Period/Gyr  a2      <a2>")')
         pr = 0.
      Do i=1+Nleft,Nn-2
         pp = 0.
         ip = 0
         do j=i-Nleft,min(i+Nright,Nn)
            If(Period(j).gt.0.)Then
               ip = ip+1
               pp = pp + Period(j)
            EndIf 
         EndDo
         pp = 2.*pp/max(ip,1)
         if(pp.gt.1.e-5)pr=pp
         write (*,*)  pp,pr
         Om = 2.*3.1415926/max(pp,pr,1.e-5)
         ss = Smooth(A2,Nn,i,iSmooth,width)
         write (2,50)Aexpn(i),Time(i),Nstep(i),pp,Om,Period(i),
     &                A2(i),ss
 50      format(f8.5,f8.4,i5,f9.4,5f10.4)
      EndDo 
      stop
      end
c-------------------------------------------------------------------- 
c               smooth A at point iPoint over iSmooth points
      FUNCTION Smooth(A,N,iPoint,iSmooth,width)
c-------------------------------------------------------------------- 
      DIMENSION A(N)
      COMMON/WEI/Weights(0:100),iFlag
      real*8    sum
      
c      write (*,*)  ' Flag=',iFlag,iSmooth,width
      If(iFlag.eq.0)Then
         iFlag =1
         sum = 0.
         do i=1,iSmooth
            sum = sum +exp(-0.5*(i/width)**2)
         EndDo 
         sum = 2.*sum +1. ! i=1 goes with weight =1
         do i=0,iSmooth
            Weights(i) = exp(-0.5*(i/width)**2)/sum
         EndDo 
      EndIf

      If(iPoint.le.iSmooth+1)Then  ! make a bit of smoothing on the edges
         Smooth =(A(Ipoint)+A(Ipoint+1))/2.
         Return
      EndIf       
      If(iPoint.ge.N-iSmooth)Then
         Smooth =(A(Ipoint)+A(Ipoint-1))/2.
         Return
      EndIf 

      sum = 0.
c      gg  = 0.
      Do i=-iSmooth,iSmooth
         sum =sum +A(i+iPoint)*Weights(iabs(i))
c         gg  =gg  +    Weights(iabs(i))
      EndDo 
      Smooth = sum
c      write(*,*) Smooth,A(iPoint),gg

      return
      end