program xbesselts
!
!     Exercise driver for subroutine besselts
!
!     B. Knapp, 2002-08-14
!
      implicit none
!
!     Variables
      integer*4 year, month, day, k, n, j, type, annular_total,
     &   hour, minute
      parameter (n=49)
      real*8 yd, yf, jdmax, gamma, u, magnitude, jd0, x(n), y(n),
     &   sind(n), cosd(n), mu(n), l1(n), l2(n), dd_dt, dmu_dt,
     &   tanf1, tanf2, jd, h
!
!     Externals
      real*8 ymd2yd
      external ymd2yd, besselts
!
      write(*,'(/a$)') ' Date of eclipse (y,m,d)? '
      read(*,*) year, month, day
      yd = ymd2yd(year,month,day)
      yf = dint(yd/1000.d0) + mod(yd,1000.d0)/365.25d0
      k = nint((yf-2000.d0)*12.3685d0)
      write(*,*) ' Computing eclipse for lunation # ',k
!
      call besselts(k, n, jdmax, gamma, u, type, annular_total,
     &   magnitude, jd0, x, y, sind, cosd, mu, dd_dt, dmu_dt, l1, l2,
     &   tanf1, tanf2)
!
      write(*, 5) '         Type:', type
      write(*, 5) 'Annular_Total:', annular_total
      write(*,10) '        JDmax:', jdmax
      write(*,10) '        Gamma:', gamma
      write(*,10) '            u:', u
      write(*,10) '    Magnitude:', magnitude
      write(*,10) '      JDEpoch:', jd0
      write(*,20) '        tanf1:', tanf1
      write(*,20) '        tanf2:', tanf2
      write(*,20) '        dD/dt:', dd_dt
      write(*,20) '       dMu/dt:', dmu_dt
      do j=1,n
         jd = jd0+(j-1)*(1.d0/144.d0)
         h = mod(jd+0.5d0, 1.d0)*24.d0
         hour = int(h)
         minute = nint(mod(h,1.d0)*60.d0)
         if (minute .eq. 60) then
            hour = hour+1
            minute = 0
         endif
         if (hour .eq. 24) then
            hour = 0
         endif
!        write(*,'(2f16.12)') h, mod(h,1.d0)*60.d0
         write(*,30) hour,minute,x(j),y(j),sind(j),cosd(j),mu(j),
     &      l1(j),l2(j)
      enddo
 5    format(a,i4)
 10   format(a,f24.8)
 20   format(a,f24.8)
 30   format(2i3.2,4f10.6,f11.5,2f10.6)
!
      end