c
c  program orbrep - finds average orbit repeat time for each GPS satellite 
c
c    reads GPS navigation message in RINEX format
c    outputs PRN number and repeat time (relative to 86400) in seconds
c
c     Reference
c
c     Agnew, D. C. and K. M. Larson 
c     Finding the Repeat Times of the GPS Constellation 
c     GPS Solutions,  Vol 10 (4), DOI 10.1007/s10291-006-0038-4, 2006 .
c
c
c     Kristine M. Larson
c     University of Colorado
c
c     implicit none
c     implicit undefined (a-z) 
c
      integer maxsat, maxeph, i, neph,k
      parameter (maxsat = 50)
      parameter (maxeph = 50)
      integer clk(maxeph,maxsat,5), it(maxeph,maxsat,5)
      real*8  bele(maxeph, maxsat, 28), 
     +    bclk(maxeph, maxsat,3), mu, pi,this_shift,
     +    clksec(maxeph, maxsat), a, n, shift

      pi= 4.d0*datan(1.d0)
      mu = 398600.5d0
c     initialize them so that first values goes in first save
      call read_broadcast3(bele, clk, bclk, clksec,it)
      do i =1, maxsat
c         if no messages at all
          if (bele(1,i,8).eq.0) goto 15
c         how many messages are there for this satellite?
          do k=1,50
            if (bele(k,i,8).eq.0) then
              neph = k-1
              goto 12
            endif 
          enddo 
12        continue
          a = 0
          n = 0
          this_shift = 0.0d0
          do k = 1, neph
            a = bele(k,i,8)*bele(k,i,8)/1000.d0
            n = dsqrt(mu/a/a/a) + bele(k,i,3)
c           print*, i, 86400.d0 -2*2*pi/n
            this_shift = this_shift + 86400.d0 -2*2*pi/n
          enddo
          shift = this_shift/neph
          write(6,'(i5, f10.2)') i, shift
15        continue
      enddo

      end
      subroutine read_broadcast3(bele, clk, bclk, clksec,it )
c     implicit none
c     implicit undefined (a-z) 
      integer maxsat, maxephem
      parameter (maxsat = 50)
      parameter (maxephem = 50)
      character*60 temp
      real*8 bele (maxephem,maxsat,28), bclk(maxephem,maxsat,3),
     .  clksec(maxephem, maxsat), rt1, rt2, rt3, rt4
      integer i, j, k, k1, clk(maxephem, maxsat,5),
     .  iprn, file(maxsat), iversion, ios,
     .  it(maxephem,maxsat,5), it1, it2, it3, it4, it5
      do i = 1, maxsat
        file(i) = 0
        do j = 1, maxephem
          do k=1,28
             bele(j,i,1) = 0
          enddo
        enddo
      enddo
      read(5, '(i6)') iversion 
c     print*, 'VERSION OF NAVIGATION FILE', iversion
      if (iversion.eq.2.or.iversion.eq.1) then
102     read(5, '(60x, a20)') temp
        if (temp(1:13).eq.'END OF HEADER'.or.
     .   temp(1:8).eq.'        ') then
c         print*, 'found end of navigation file header'
        else
          goto 102
        endif
      else
        print*, 'can only read version 1 or 2'
        call exit(1)
      endif
12    read( 5, '(i2, 5i3, f5.1, 3d19.12)',err=14,iostat=ios) iprn, it1,
     .  it2, it3, it4, it5, rt1, rt2, rt3, rt4
      if (ios.ne.0) goto 14
c     write(6, *) iprn, it1, it2 
      file(iprn) = file(iprn) + 1
      it(file(iprn), iprn, 1) = it1
      it(file(iprn), iprn, 2) = it2
      it(file(iprn), iprn, 3) = it3
      it(file(iprn), iprn, 4) = it4
      it(file(iprn), iprn, 5) = it5
      if (file(iprn).gt.50) then 
        print*, 'MORE THAN 50 ephemris for PRN', iprn
        goto 14
      endif
      clk(file(iprn), iprn, 1) = it1
      clk(file(iprn), iprn, 2) = it2
      clk(file(iprn), iprn, 3) = it3
      clk(file(iprn), iprn, 4) = it4
      clk(file(iprn), iprn, 5) = it5
      clksec(file(iprn), iprn) = rt1
      bclk( file(iprn), iprn, 1) = rt2
      bclk( file(iprn), iprn, 2) = rt3
      bclk( file(iprn), iprn, 3) = rt4
      do j = 1, 6
        k1 = 4*(j-1) + 1
        read( 5, '(3x, 4d19.12)', err=14) 
     .       ( bele(file(iprn), iprn , k), k = k1, k1+3)
      enddo
      if (iversion.eq.2) read(5,'(a60)') temp
      goto 12
14    continue
      return
      end