program getarrtime

c computes pulse arrival time for Crab, given published ephem data
c derived from "werner.f" by MK & CAJ  June 2004
c
      implicit none

      real*8 epoch, secs, v, vdot
      integer i/0/
      real*8 treq, ti, p, pdot, vddot, freq
      real*8 turns, pmean, tarr
      character*20 eph

      do i=1,2
         if (i .eq. 1) then
            eph ='Aug 03'
            epoch=52866d0+0.018825/86400d0
            v = 29.8034282349d0
            vdot = -373501.93d-15            
            
         elseif (i .eq. 2) then
            eph ='Sep 03'
            epoch=52897d0+.002531/86400d0
            v = 29.8024278709d0
            vdot = -373484.63d-15
         endif
      
         treq = 52897.d0
                  
         ti = (treq-epoch)*86400.d0
         p = 1.0d0/v
         pdot = -vdot/(v*v)
         vddot = 2.0d0*pdot*pdot/(P*P*P)
         
         freq = v + ti*vdot + ti*ti/2.0d0*vddot
         turns = v*ti + vdot*ti*ti/2.0d0 + vddot*ti*ti*ti/6.0d0
c     write(*,*) turns
         turns = dmod(turns,1.0d0)
c     write(*,*) turns
         if (treq-epoch.lt.0.0d0) turns=turns+1.0
c     write(*,*) turns
         
         pmean = 1.0d0/freq
         
         tarr = treq - dmod(turns,1.0d0)*pmean/86400.0d0
         
         write(*,*) ' Eph  : ',eph
         write(*,*) ' P    : ',p
         write(*,*) ' Pmean: ',pmean
         write(*,*) ' nu   : ',freq
         write(*,*) ' turns: ',turns
         write(*,*) ' t_arr: ',int(tarr), (tarr-int(tarr))*86400d0
      enddo
      end