      program mcalistr
      integer day,month,year
      character text*80
      real*8 date,gmst,jdate,bepoch,lyear,jd1900
      lyear=365.242198781D0
      jd1900=2415020.31352D0
      degrad=3.141592654/180.0
      open(7,file='mca.dat')
   93 write(6,*)'Enter Besselian epoch, pa, +/- pa, separation, +/- r!'
      read(5,2,end=111)text
    2 format(a80)
      if(index(text,'n').ne.0)goto 111
      read(text,*,end=93,err=111)bepoch,pa,epa,sep,esep
      jdate=((bepoch-1900.D0)*lyear)+jd1900
  911 year=1950
      month=1
      day=1
      call julda(year,month,day,date,gmst)
      if(jdate.ge.date)then
       isign=1
      else
       isign=-1
      endif
      do i=1,100
       year=1950+(i-1)*isign
       call julda(year,month,day,date,gmst)
       ndays=int(jdate-date)
       if(abs(ndays).lt.365)goto 100
      enddo
  100 if(isign.lt.0)then
       year=year-1
       call julda(year,month,day,date,gmst)
       ndays=int(jdate-date)
      endif
      do i=1,ndays
       call nextday(day,month,year)
      enddo
      fracod=jdate-int(jdate)-0.5
      if(fracod.lt.0)fracod=fracod+1.0
      ihour=int(fracod*24.0)
      imin=int((fracod*24.0-float(ihour))*60.0)
      write(6,1)year,month,day,ihour,imin
    1 format(' Year= ',i4,' month= ',i2,' day= ',i2,', 'i2,' UT ',i2)
      utsec=ihour*3600.+imin*60.
      write(7,11)sep,pa,esep,sin(epa*degrad)*sep,pa,year-1900,month,day,
     $           utsec,zero,zero
   11 format(f5.1,1x,f5.1,1x,f4.1,1x,f4.1,1x,f5.0,1x,i2.2,i2.2,i2.2,
     $       1x,f6.0,3x,f3.1,1x,f3.1)
      goto 93
  111 end 
      
      SUBROUTINE NEXTDAY ( IDAY, IMON, IYEAR )
C
C   Increments the date by one
C
      IMPLICIT UNDEFINED(A-Z)
      INTEGER *4 IDAY, IMON, IYEAR, DAYS(12,0:1), LEAP
      DATA DAYS / 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31,
     $            31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /
      LEAP = MIN ( 1, MOD( IYEAR, 4 ) )
      IDAY = IDAY + 1
      IF ( IDAY .GT. DAYS(IMON,LEAP) ) THEN
          IMON = IMON + 1
          IDAY = 1
          IF ( IMON .EQ. 13 ) THEN
              IMON = 1
              IYEAR = IYEAR + 1
          END IF
      END IF
      RETURN
      END
        
      
        SUBROUTINE JULDA(YEAR,MONTH,DAY,DATE,GMST)
C
C               Returns DATE = Julian Date and GMST = Greenwich Mean
C               sidereal time at 0 hrs U.T. on day DAY/MONTH/YEAR.
C               TJP 1977 August 5
C               NB Accuracy not tested; GMST is correct at least to
C               nearest second
C
        INTEGER DAY,MONTH,YEAR
        INTEGER MOFF(12)
        REAL*8 DATE,GMST,DU,SMD,T,DFLOAT,DMOD
        DATA MOFF/0,31,59,90,120,151,181,212,243,273,304,334/
C
C               JD number at 12 hrs UT on Jan 0 of year YEAR:
C               (Gregorian Calendar)
C
        NYRM1=YEAR-1
        IC=NYRM1/100
        JD=1721425 + 365*NYRM1 + NYRM1/4 - IC + IC/4
C
C               Number of days from Standard Epoch 1900 Jan 0.5
C               (JD 2415020.0) to Jan 0.0 of YEAR:
C
        DU=DFLOAT(JD-2415020) - 0.5D0
C
C               Day number; is it a leap year?
C
        NDAYS=MOFF(MONTH) + DAY
        IF(MONTH.GT.2 .AND. MOD(YEAR,4).EQ.0 .AND. MOD(YEAR,100).NE.0)
     1     NDAYS=NDAYS+1
C
C               UT from Epoch to today (days), (centuries)
C
        SMD=DU+DFLOAT(NDAYS)
        T=SMD/36525.D0
C
C               Greenwich Mean Sidereal Time
C
        GMST= DFLOAT(6*3600 + 38*60) +45.836D0
     1     + 8 640 184.542D0*T  +  0.0929D0*T**2
        GMST=DMOD(GMST,86400D0)
C
C
        DATE=2415020.D0+SMD
C
        RETURN
        END


