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 /Gyr Omega*Gyr", & " Period/Gyr 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