      SUBROUTINE READFILE
C
C  Adapted from SOLVE (David Mozurkewich) by
C
C  J.T. Armstrong     2 Jan 91
C
C  This subroutine takes a list of .CAL files and reads data for 
C  specified stars and filters.  
C-----------------------------------------------------------------------

      INCLUDE 'BINFIT.INC'

      INTEGER*4    INCAT, INDAT, IDAY, IMON, IYEAR
      INTEGER*4    NBODY, ILEN, IERR, IDATE(MXREC)
      INTEGER*4    NDATREAD, LENGTH, LEN1,LEN2
      REAL*8       B(3), U(3), V(3), W(3)
      REAL*8       BX(MXBAS), BY(MXBAS), BZ(MXBAS)
      REAL*8       RAM, DECM, PMRA, PMDEC, PAR, RV, RAA, DECA
      REAL*8       TDTOFF, DJUL, UT1, GAST, RA, DEC, HA
      CHARACTER*64 FILELIST, BFILENAME
      LOGICAL*1    ONEFILE
      REAL*8       JD, DOT
      EXTERNAL     JD, DOT, LEN1

      PARAMETER (INCAT=11)        ! Star catalog input unit
      PARAMETER (INDAT=12)        ! Data input unit
      PARAMETER (TDTOFF=55.D0/86400.D0)   ! Time offset late 1990
      real*8 mepoch(3),errsum(3),weight(3)

C-----------------------------------------------------------------------
C  Open baseline data file
C
D     WRITE(OUTC,'(A)') ' Opening BASELINE.DAT data file: '
      OPEN(UNIT=INDAT,FILE='e:\\catalogs\\fk5\\BASELINE.DAT',
     $               STATUS='OLD',IOSTAT=IERR)
      IF ( IERR .NE. 0 ) THEN
          WRITE (OUTC,'(A)')
     $           ' File BASELINE.DAT not in current directory'
          WRITE (OUTC,'(A,$)') ' Provide name of baseline file:'
          READ (INC,'(A)',IOSTAT=IERR) BFILENAME
          OPEN(UNIT=INDAT,FILE=BFILENAME,STATUS='OLD',IOSTAT=IERR)
          IF (IERR .NE. 0) THEN
              WRITE(OUTC,'(/,A)') ' Baseline file not found'
              STOP
         END IF
      END IF
C
      READ(INDAT,'(1(/))') ! skip header lines in file
      DO I = 1, MXBAS
          READ(INDAT,'(I5,3(F15.6))',IOSTAT=IERR)
     $                  J, BX(I), BY(I), BZ(I)
         IF ( (IERR .NE. 0) .OR. (J .NE. I) ) THEN
            WRITE(OUTC,'(A,$)') ' BASELINE FILE IS CORRUPT' 
            STOP
         END IF
      END DO
      CLOSE(INDAT)
D     WRITE(OUTC,'(A)') ' Baseline file successfully read:'
C-----------------------------------------------------------------------
C     Read the list of input files.
C
      IERR = -1
      ONEFILE = .FALSE.
      DO WHILE ( IERR .NE. 0 ) 
          WRITE (OUTC,'(A,$)') 
     1      ' Input file with list of .CAL files, or single .CAL file: '
          READ  (INC,'(A)',IOSTAT=IERR) FILELIST
          LENGTH = LEN1(FILELIST)
          LEN2 = LENGTH-2
          IF(FILELIST(LEN2:LENGTH).EQ.'CAL' .OR.
     1       FILELIST(LEN2:LENGTH).EQ.'cal'      ) THEN
              ONEFILE = .TRUE.
              INFILE(1) = FILELIST
              IERR = 0
          END IF
          IF ( IERR .EQ. 0 .AND. .NOT.ONEFILE) THEN
              OPEN ( UNIT=3,FILE=FILELIST, STATUS='OLD', IOSTAT=IERR )
              IF ( IERR .NE. 0 ) THEN 
                  WRITE (OUTC,'(A,A,A)')  ' Cannot open FILE LIST ',
     1                    FILELIST,'.  Try again.'
D             ELSE
D                 WRITE (OUTC,*) ' .CAL file list read successfully'
              END IF
              ILEN = INDEX(FILELIST, ' ' ) - 1
D             WRITE(OUTC,*)  ' File ', FILELIST(1:ILEN), ' opened.'
          ELSE IF ( IERR .NE. 0 .AND. .NOT.ONEFILE ) THEN
              WRITE(OUTC,*) ' Error reading name of .CAL file list.'
          END IF
      END DO
      NFILES = 1
      IERR = 0
      IF(.NOT.ONEFILE) THEN
          DO WHILE ( (NFILES .LT. MXFILES) .AND. (IERR .EQ. 0) )
              READ(3,'(A)',IOSTAT=IERR) INFILE(NFILES)
              IF ( IERR .EQ. 0 ) THEN
                  OPEN(UNIT=4,FILE=INFILE(NFILES),STATUS='OLD',
     1                IOSTAT=IERR)
                  IF ( IERR .EQ. 0 ) THEN
                      K = INDEX ( INFILE(NFILES), ' ' )
                      WRITE(OUTC,'(A,I5,2X,A,A)') ' File ',NFILES,
     1                    INFILE(NFILES)(1:K),' added to list.'
                      NFILES = NFILES + 1
                      CLOSE(4)
                  ELSE
                      K = INDEX ( INFILE(NFILES), ' ' )
                      WRITE(OUTC,'(A,7X,A,A)') ' File ',
     1                    INFILE(NFILES)(1:K),' not found.'
                      IERR = 0
                  END IF
D             ELSE
D                 WRITE(OUTC,*) ' End of .CAL list encountered.' 
              END  IF
          END DO
          CLOSE (3)
          NFILES = NFILES - 1
      END IF
C-----------------------------------------------------------------------
C         Get star number to use: 
C
      WRITE (OUTC,'(A)') '$Enter the star number to use: '
      READ (INC,*) STARNO
D     WRITE(OUTC,'(A,I5,2X,A)') ' Looking for star number ',STARNO,
D    1    STARNAME
      CALL READFK5( STARNO, STARNAME, RAM, DECM, PMRA, PMDEC,
     1                    PAR, RV )
D     WRITE(OUTC,*) ' RA, Dec, PMRA, PMDec, Parallax, RV:'
D     WRITE(OUTC,'(1X,F5.2,2X,F6.2,2X,F4.2,2X,F6.2,2X,F6.4,2X,F5.1)')
D    1    RAM, DECM, PMRA, PMDEC, PAR, RV

C----------------------------------------------------------------------
C  Get all the data from all of the files.
C
      K = 1
cc    NFILT = 0
      NDATA = 0
      DO J = 1, NFILES
          nfilt=0
          CALL READ_VIS (INFILE(J),STARNO,IDATE(K),
     1                   NDATA, NDATREAD)
D         WRITE(OUTC,'(A,I2,A,I7)') ' IDATE(',K,') = ',IDATE(K)
          IF ( NDATREAD .GT. 0 ) THEN
              NDATA = NDATA + NDATREAD
              IDAY  = MOD(IDATE(NDATA-NDATREAD+1),       100)
              IMON  = MOD(IDATE(NDATA-NDATREAD+1)/100,   100)
              IYEAR = MOD(IDATE(NDATA-NDATREAD+1)/10000, 100) + 1900
              WRITE(OUTC,'(A,I4,A,I10,3I5,A,I3)') ' Scan',
     1           NDATA-NDATREAD+1,' Date(YYMMDD) = ',
     2           IDATE(NDATA-NDATREAD+1),IYEAR, IMON, IDAY, 
     3           '; baseline = ', BASELINE(NDATA-NDATREAD+1)
          END IF
      END DO
      WRITE(OUTC,'(A,I5,A,I5)')  ' Star = ', STARNO,
     1                        '; number of data points = ', NDATA
C      IF ( NDATA .LT. 3 ) GO TO 200
C-------------------------------------------------------------------
C  Calculate frequency band parameters:
C     WRITE(OUTC,*) ' Enter calibration uncertainties for:'
      DO IG = 1, NFILT
          LAM0(IG) = LAMBDA0(IG) * 1.0D-9       ! Convert to meters
          DLAM(IG) = DLAMBDA(IG) * 1.0D-9
          NU0(IG)  = C/LAM0(IG)
          NULO(IG) = C / ( LAM0(IG) + DLAM(IG)/2.0D0 )
          NUHI(IG) = C / ( LAM0(IG) - DLAM(IG)/2.0D0 )
          DNU(IG)  = ( NUHI(IG) - NULO(IG) ) / NUSTEPS
C         ----------------------------------------------------------
C         Needed for maximum entropy calculation:
C         WRITE(OUTC,'(1X,F5.1,A)') LAMBDA0(IG),' nm: '
C         ----------------------------------------------------------
C         READ(INC,*) SIGCAL(IG)
D         WRITE(OUTC,'(A,I2,A,3(1X,G9.3))') 
D    1              ' Filter ',IG,' freq ctr, hi, lo: ',
D    1              NU0(IG),NULO(IG),NUHI(IG)
      END DO
D     WRITE(OUTC,*) ' Wavelengths, frequencies calculated'
C-------------------------------------------------------------------
C  Calculate the (u,v) position for each observation.
C
      DO J = 1, NDATA
          IF (J.EQ.1) THEN
              IDAY  = MOD(IDATE(J),       100)
              IMON  = MOD(IDATE(J)/100,   100)
              IYEAR = MOD(IDATE(J)/10000, 100) + 1900
D             WRITE(OUTC,1300) J, IDATE(J), IDAY, IMON, IYEAR
 1300         FORMAT( ' SCAN = ', I5, ' DATE = ', I10, 3I5,/ )
              DJUL  = JD ( IDAY, IMON, IYEAR )
D             WRITE(OUTC,'(A,F20.5)') ' DJUL = ',DJUL
              B(1)  = BX(BASELINE(J))     ! Baseline components in meters
              B(2)  = BY(BASELINE(J))
              B(3)  = BZ(BASELINE(J))
              IF(B(1).EQ.0..AND.B(2).EQ.0..AND.B(3).EQ.0.) THEN
                  WRITE(OUTC,'(A,I3,A)')
     1                ' Baseline ',BASELINE(J),' has zero length!  Bye.'
                  CALL EXIT
              END IF
D             WRITE(OUTC,'(A,3(1X,G12.5))') ' B(1,2,3) (m): ',
D    1                                       B(1),B(2),B(3)
          ELSE IF ( (IDATE(J).NE.IDATE(J-1))
     $                  .OR. (BASELINE(J).NE.BASELINE(J-1)) ) THEN
              IDAY  = MOD(IDATE(J),       100)
              IMON  = MOD(IDATE(J)/100,   100)
              IYEAR = MOD(IDATE(J)/10000, 100) + 1900
              WRITE(OUTC,1300) J, IDATE(J), IDAY, IMON, IYEAR
              DJUL  = JD ( IDAY, IMON, IYEAR )
D             WRITE(OUTC,'(A,F20.5)') ' DJUL = ',DJUL
              B(1)  = BX(BASELINE(J))     ! Baseline components in meters
              B(2)  = BY(BASELINE(J))
              B(3)  = BZ(BASELINE(J))
              IF(B(1).EQ.0..AND.B(2).EQ.0..AND.B(3).EQ.0.) THEN
                  WRITE(OUTC,'(A,I3,A)')
     1                ' Baseline ',BASELINE(J),' has zero length!  Bye.'
                  CALL EXIT
              END IF
D             WRITE(OUTC,'(A,3(1X,G12.5))') ' B(1,2,3) (m): ',
D    1                                       B(1),B(2),B(3)
          END IF
          UT1  = HOURS(J) / 24.D0
D         WRITE(OUTC,'(A,I2,A,G20.10)') ' HOURS(',J,') = ',HOURS(J)
C         WRITE(OUTC,'(A,I4,A)') ' Calculating TJD(',J,')'
          TJD(J)  = DJUL + UT1 + TDTOFF
D         WRITE(OUTC,'(A,I2,A,F20.5)') ' TJD(',J,') = ',TJD(J)
          NBODY = 3
D         WRITE(OUTC,*) ' Going to APSTAR'
D         WRITE(OUTC,*) ' TJD, NBODY, RAM, DECM:'
D         WRITE(OUTC,998) TJD(J), NBODY, RAM, DECM
D         WRITE(OUTC,*) ' PMRA, PMDEC, PAR, RV:'
D         WRITE(OUTC,999) PMRA, PMDEC, PAR, RV
  998     FORMAT(1X,F10.1,2X,I2,2X,F5.2,2X,F6.2)
  999     FORMAT(1X,F4.2,2X,F6.2,2X,F6.4,2X,F4.0)
          CALL APSTAR (TJD(J), NBODY, RAM, DECM, PMRA, PMDEC, PAR,
     $                         RV, RAA, DECA)
D         WRITE(OUTC,*) ' Going to SIDTIM'
          CALL SIDTIM ( DJUL, UT1, 1, GAST )
D         WRITE(OUTC,'(A,G20.10)') ' GAST = ',GAST
D         WRITE(OUTC,*) ' Going to DIURN'
          CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC )
D         WRITE(OUTC,'(A,G20.10,1X,G20.10,1X,G20.10)') 
D    1                             ' GAST, RA, LONG = ',GAST,RA,LONG
          HA = 15.D0 *( GAST - RA ) - LONG
          IF ( HA .LT. -180 ) HA = HA + 360.D0
D         WRITE(OUTC,'(A,G20.10)') ' HA (deg) = ', HA
D         WRITE(OUTC,'(A,G20.10)') ' HA*DEGRAD = ',HA*DEGRAD
C         Components of u, v, w axes in local coordinate system
          U(1) = -SIN(HA*DEGRAD)
          U(2) = -COS(HA*DEGRAD)
          U(3) =  0.D0
D         WRITE(OUTC,'(A,3(1X,G12.5))') ' U(1,2,3) : ',U(1),U(2),U(3)
          V(1) =  SIN(DEC*DEGRAD) * COS(HA*DEGRAD)
          V(2) = -SIN(DEC*DEGRAD) * SIN(HA*DEGRAD)
          V(3) = -COS(DEC*DEGRAD)
D         WRITE(OUTC,'(A,3(1X,G12.5))') ' V(1,2,3) : ',V(1),V(2),V(3)
          W(1) = -COS(DEC*DEGRAD) * COS(HA*DEGRAD)
          W(2) =  COS(DEC*DEGRAD) * SIN(HA*DEGRAD)
          W(3) = -SIN(DEC*DEGRAD)
D         WRITE(OUTC,'(A,3(1X,G12.5))') ' W(1,2,3) : ',W(1),W(2),W(3)
          DO IG = 1, NFILT
              UDATA(J,IG) = DOT( B, U ) / (LAMBDA0(IG)*1.0E-9)
              VDATA(J,IG) = DOT( B, V ) / (LAMBDA0(IG)*1.0E-9)
              DELAY(J,IG) = DOT( B, W ) / (LAMBDA0(IG)*1.0E-9)
D         WRITE(OUTC,'(A,I4,1X,I2,3(1X,G12.5))') 
D    1               ' Scan, Filter, U, V, W (wavel): ', J, IG,
D    2               UDATA(J,IG),VDATA(J,IG),DELAY(J,IG)
          END DO
      END DO
D     WRITE(OUTC,*) ' U, V coordinates calculated'
  200 CONTINUE
C Count the valid data
      do ig=1,nfilt
      errsum(ig)=0.
      mepoch(ig)=0.
      enddo
      NZTOT = 0
      DO IG = 1, NFILT
          NZVIS(IG) = 0
          DO ID = 1, NDATA
              IF ( VISERR(ID,IG)  .NE. 0.0
     1                      .AND.
     2             OBSVIS2(ID,IG) .NE. 2.0 ) THEN
                  NZVIS(IG) = NZVIS(IG) + 1
                  weight(ig)=1./viserr(id,ig)/viserr(id,ig)
                  mepoch(ig)=mepoch(ig)+tjd(id)*weight(ig)
                  errsum(ig)=errsum(ig)+weight(ig)
              END IF
          END DO
          NZTOT = NZTOT + NZVIS(IG)
      END DO
      mepoch3=(mepoch(1)+mepoch(2)+mepoch(3))/
     1        (errsum(1)+errsum(2)+errsum(3))
      write(outc,*)' Mean epoch= ',mepoch3
D     WRITE(OUTC,*) ' Number of points read = ', NDATREAD
D     WRITE(OUTC,*) ' Number of valid points = ', 
D    1              (NZVIS(I),I=1,NFILT),  NZTOT

      RETURN
      END
C======================================================================
      REAL *8 FUNCTION DOT( A, B, )
C
C  Do a vector dot product.
C
      IMPLICIT UNDEFINED ( A-Z )
      REAL *8 A(3), B(3)
      DOT = A(1)*B(1) + A(2)*B(2) + A(3)*B(3)  
      RETURN
      END
C=====================================================================
      SUBROUTINE SORT ( LIST, N )
C
C   Sort a list of integers and remove duplicate entries.
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 LIST(*), N, I, TEMP, I0, I1
      IF ( N .LE. 1 ) RETURN
      I0 = 1
  100 CONTINUE
          I = I0+1
  200     CONTINUE
              IF ( LIST(I) .EQ. LIST(I0) ) THEN
                  LIST(I) = LIST(N)
                  N = N - 1
              ELSE IF ( LIST(I) .LT. LIST(I0) ) THEN
                  TEMP = LIST(I0)
                  LIST(I0) = LIST(I)
                  LIST(I) = TEMP
                  I = I + 1
              ELSE
                  I = I + 1
              END IF
          IF ( I .LE. N ) GO TO 200
          I0 = I0 + 1
          IF ( I0 .LT. N ) GO TO 100
      RETURN
      END
C=====================================================================
      REAL *8 FUNCTION JD ( IDAY, IMON, IYEAR )
C
C           Returns the julian day for 0 hours ut.
      IMPLICIT UNDEFINED(A-Z)
      SAVE
      INTEGER *4 IDAY, IMON, IYEAR, YEAR0, YEAR1, DOY
      PARAMETER (YEAR0=1986)
      PARAMETER (YEAR1=1995)
      REAL *8 JD0(YEAR0:YEAR1)
      EXTERNAL DOY
      DATA JD0 / 2446430.5D0,  2446795.5D0, 2447160.5, 2447526.5,
     $           2447891.5D0,  2448256.5D0, 2448621.5, 2448987.5,
     $           2449352.5D0,  2449717.5D0 /
      IF (( IYEAR .LT. YEAR0 ).OR.( IYEAR .GT. YEAR1 ) ) THEN
          WRITE(0,'(A,I5)') ' INVALID YEAR FOR COMPUTING JD:', IYEAR
      ELSE
          JD = JD0(IYEAR) + DOY( IDAY, IMON, IYEAR )
      END IF
      RETURN
      END
C======================================================
      INTEGER *4 FUNCTION DOY( IDAY, IMON, IYEAR )
C
C     Returns Day of Year for date in arguements ( 1 Jan => doy=1 )
      IMPLICIT UNDEFINED(A-Z)
      SAVE
      INTEGER *4 IDAY, IMON, IYEAR, LEAP, I, NDAY(12,0:1)
      DATA NDAY / 0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30,
     $            0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30 /
C        LEAP=0 for leap lear
      LEAP = MIN( 1, MOD( IYEAR, 4 ) )
      DOY = IDAY
      DO 100 I = 1, IMON
          DOY = DOY + NDAY(I, LEAP)
  100 CONTINUE
      RETURN
      END
C=======================================================
      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
C=======================================================
      SUBROUTINE DIURN( TJD, GAST, RAA, DCA, RAT, DCT )
C
C           Calculates diurnal abberation
C    TJD   Julian date ( in Terestial dynamical time )
C    GAST  Greenwich aparent sideral time
C    RAA   Apparent geocentric right ascen.
C    DCA   Apparent geocentric decl.
C    RAT   Topocentric right ascen. ( output )
C    DCT   Topocentric decl. ( output )
C
      IMPLICIT UNDEFINED(A-Z)
      SAVE
      REAL *8 PI, TWOPI, RADCON, C, ERAD, AU, VROT, DLAT, DLON, HT
      REAL *8 SINLAT, COSLAT, RHO, ROTVEL, ABRCON
      REAL *8 TJD, RAA, DCA, RAT, DCT, GAST, HA
      REAL *8 SINHA, COSHA, SINDEC, COSDEC, ALT, SINALT
      PARAMETER ( PI = 3.141592653589793D0 )
      PARAMETER ( TWOPI = 2.D0*PI )
      PARAMETER ( RADCON = PI / 180.D0 )
      PARAMETER ( C = 2.99792538D8 )
      PARAMETER ( ERAD = 6378140.D0 )
      PARAMETER ( AU = 1.49597870D11 )
      PARAMETER ( VROT = 7.2921151467D-5 )
C           Mt Wilson geocentric coordinates.
      PARAMETER ( DLAT =  34.2166944D0 )
      PARAMETER ( DLON = 118.0591670D0 )
      PARAMETER ( HT   = 1742.0D0 )
      LOGICAL FIRST
      DATA FIRST / .TRUE. /
C
      IF ( FIRST ) THEN 
          FIRST  = .FALSE.
          SINLAT = DSIN ( DLAT * RADCON )
          COSLAT = DCOS ( DLAT * RADCON )
          RHO = ( ERAD + HT  - 21400.D0 * SINLAT*SINLAT ) / ERAD
          ROTVEL = VROT * RHO * ERAD
          ABRCON = ROTVEL * COSLAT / C
      END IF
      HA  = ( (GAST - RAA ) * 15.D0 - DLON ) * RADCON
      SINHA  = DSIN( HA )
      COSHA  = DCOS( HA )
      SINDEC = DSIN( DCA * RADCON )
      COSDEC = DCOS( DCA * RADCON )
C
C         Test for object above the hroizon
C
      SINALT = SINLAT * SINDEC + COSLAT * COSDEC * COSHA
      ALT    = DASIN ( SINALT ) / RADCON
      IF ( ALT .LE. -1.D0 ) THEN
          RAT = 99.9D0
          DCT = 99.9D0
      ELSE
C           Apply correction for dirunal abberation
          RAT = RAA + ( ABRCON * COSHA / COSDEC ) /(15.D0*RADCON)
          DCT = DCA + ( ABRCON * SINHA * SINDEC ) / RADCON 
      END IF
      RETURN
      END


