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

      Open(1,file='Stat.dat')
      Open(2,file='omegap2.dat')
      Nn       =1
      Nstep0   =0
      Nleft    =(Nave-1)/2
      Nright   = Nleft +Nave
      iFlag    = 0
      iTest   = 79
 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)
c      Time(Nn) =Aexpn(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
      Do i=2,Nn
         If(Time(i).lt.Time(i-1))Then
            write(*,'(a,2i5,3f9.3)')' ERROR: wrong order of snapshots:',
     &                  i,Nstep(i),Time(i)
            do j=0,5
               if(i-j.gt.1)write (*,'(a,i5,a,g12.4)')
     &               ' iStep=',Nstep(i-j),' Time=',Time(i-j)
            EndDo 
            Stop
         EndIf 
      EndDo 


      Phi(1) = Angle(1)
      phase  = 0.
      Do i=2,Nn            ! get phi(t)
         if(Angle(i)-Angle(i-1).lt.aJump) phase= phase+1.
      EndDo 
      Om_aver = phase*180./Time(Nn)
      phase =0
      Do i=2,Nn            ! get phi(t)
         if(Angle(i)-Angle(i-1).lt.aJump) phase= phase+1.
         Phi(i) =Angle(i) + phase*180.-Om_aver*Time(i)
      EndDo 

         
      Do i=2,Nn-2            ! find periods
         jstart = max(i-5,1)  ! try to center time around the moment
         ang0 = Angle(jstart)        ! starting Angle
         if(Angle(jstart)-Angle(jstart-1).lt.aJump) ang0= ang0+180.
         per  = 0.              ! period
         ang1 = ang0            !  angle on previous snapshot
          da  = 0.              ! offset for angle
         !If(i.eq.iTest)
c         write (*,'(" i=",2i5," Angle0=",f8.4)')i,jstart,ang0
         Do j=jstart,Nn
            ang = Angle(j)
            if(Angle(j)-Angle(j-1).lt.aJump)da =da +180.
             ang=ang+da
c             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(jstart)+
     &               (ang0+180.-ang1)/(ang-ang1)*(Time(j)-Time(j-1))
c      If(i.eq.iTest)write(*,'(15x,5f9.4)')per,Time(j-1),Time(jstart)
               goto 30
            EndIf 
            ang1 = ang
         EndDo 
         
         per = (Time(j-1)-Time(jstart))*(180./(ang-ang0))
 30      Period(i) =per
         If(per.le.0.)write (*,*)  ' err:',i,ang0,per
      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
c         write (*,*)  pp,pr
         Om = 2.*3.1415926/max(pp,pr,1.e-5)
         ss = Smooth(A2,Nn,i,iSmooth,width)
         Om2 =OmegaP(Phi,Time,Nn,i,150)+Om_aver*pi180
         write (2,50)Aexpn(i),Time(i),Nstep(i),pp,Om,2.*Period(i),
     &                A2(i),ss,Om2,Phi(i),Angle(i)
 50      format(f8.5,f8.4,i5,f9.4,5f10.4,2g16.6)
      EndDo 
      stop
      end
c-------------------------------------------------------------------- 
c               least square fit for Omega_p
      FUNCTION OmegaP(Phi,Time,N,iPoint,iSmooth)
c-------------------------------------------------------------------- 
      PARAMETER(pi180 =3.141592653/180.)
      Real*8 Phi(N),st,st2,ps,stp
      Real   Time(N) 

      st  =0.
      st2 =0.
      ps  =0.
      stp =0.
      iStart  = max(iPoint-iSmooth,1)
      iFinish = min(iPoint+iSmooth,N)
      Npoints = iFinish-iStart+1

      Do i=iStart,iFinish
         st  = st +Time(i)
         st2 = st2+Time(i)**2
         ps  = ps +Phi(i)
         stp = stp+Time(i)*Phi(i)
      EndDo 
      
      Omegap = (stp/Npoints - ps*st/Npoints**2)
     &         / (st2/Npoints -(st/Npoints)**2)*pi180

c      write(*,'(5x,3i5,9g12.4)')iPoint, iStart,iFinish,Time(iPoint),
c     &  Phi(iPoint),st/Npoints,ps/Npoints,Omegap,
c     &  (Phi(iFinish)-Phi(iStart))/(Time(iFinish)-Time(iStart))*pi180
      Return
      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
