subroutine readvis
c
c     Reads concatenated .CAL file into common block /data/
c     Reads filter information or asks user for filter widths
c     Also multiplies visibilities and uncertainties with
c     values optionally indicated in the .CAL file by letter 'M'.
c
      save nnfilt,filter,fwidth
c
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'

      integer mxfk5
      parameter (mxfk5=2500)  ! Maximum number of stars in diameter.dat

      REAL*8 RAM, DECM, PMRA, PMDEC, PAR, RV
      real*8 filter(mmxfilt),fwidth(mmxfilt),exfilter(mmxfilt)
      real*8 time,lambda1,lambda2,lambda3,lamhi,lamlo,lamst
      real*8 mv(3),me(3)
      integer*4 exstar(mxstars)
      character data*256,fknam*6,digits(10)
      integer*4 basel,datum,fk
      real*4 rdatum
      logical foundnum,calibs,query
c
      data digits/'0','1','2','3','4','5','6','7','8','9'/
      data nnfilt/0/
c
c ... Set up filter information
      if(nnfilt.ne.0)goto 63 ! there was a call to READVIS before
      do j=1,mmxfilt         ! initialize filter information array
       do i=1,50             ! number of transmission measurements
        do n=1,2             ! 1=wavelength [nm], 2=%Transmission
         ftcurve(j,i,n)=0.d0
        enddo
       enddo
      enddo
      open(1,file='/home/chummel/mark3/binafit/ftcurve.dat',
     $     status='old',err=61)
      do 60 i=1,1000
       read(1,1,end=62)data
       ipos=index(data,'filter')
       if(ipos.ne.0)then
        nnfilt=nnfilt+1
        read(data(ipos+6:256),*)filter(nnfilt)
        fwidth(nnfilt)=nnfilt  ! serves as index to array ftcurve
        if(nnfilt.gt.mmxfilt)then
         write(outc,*)
     &   '*** ERROR: Maximum number of tr. curves exceeded!'
         stop
        endif
        j=0
        goto 60
       endif
       j=j+1
       if(j.gt.50)then
        write(outc,*)'Maximum number of tr.curve entries exceeded!'
        stop
       endif
       read(data,*)ftcurve(nnfilt,j,1),ftcurve(nnfilt,j,2)
   60 continue

   62 write(outc,*)'Read filter information successfully!'
      close(1)
      goto 63
   61 write(outc,*)'Proceeding without file FTCURVE.DAT!'
   63 continue
c
c ... Now enter file name for calibrated visibilities
      write(outc,*)
   10 write(outc,'(a,$)')' Enter name of .CAL file: '
      read(inc,1,end=10)data
    1 format(a256)
      open(1,file=data,err=10,status='old')
      write(outdat,'(a,a)')'! Data file name: ',data(1:60)
c
c ... Enter list of stars
      if(starno.eq.0)then
   20  write(outc,'(a)')       ' FK5 number of binary, '
       write(outc,'(a,$)')'  or list of calibrators: '
       read(inc,1)data
       foundnum=.false.
       innums=0
       do 40 i=1,80
       do j=1,10
        if(data(i:i).eq.digits(j))then
         foundnum=.true.
         goto 40
        endif
       enddo
       if(foundnum)then
        innums=innums+1
        foundnum=.false.
       endif
   40  continue
       if(innums.gt.1)then
        calibs=.true.
       else
        calibs=.false.
       endif
       do k=1,mxstars
        calstars(k)=0
       enddo
       read(data,*,err=20)(calstars(k),k=1,innums)
       if(.not.calibs)then
        if(starno.ne.0.and.starno.ne.calstars(1))then
         write(outc,*)'*** ERROR: selected stars in data and model ',
     &                'not the same!'
         goto 20
        endif
        starno=calstars(1)
       endif
      else
       innums=1
       calstars(1)=starno
       calibs=.false.
      endif
c
c ... Some stars are not calibrators in every filter (e.g. Halpha)
      if(calibs)then
       if(query('Exclude data for a specific star and filter ? (y/n) '))
     &                                     then
        write(outc,*)'Enter star and filter, end with ctrl z! '
        ierr=0
        iex=0
        do while(ierr.eq.0)
         iex=iex+1
         read(inc,*,iostat=ierr)exstar(iex),exfilter(iex)
        enddo
        iex=iex-1
c       rewind inc
       endif
      endif
c
c ... Get stellar diameter from diameter.dat file
      open(3,file='/home/chummel/mark3/binafit/diameter.dat',
     $     status='old',err=50)
      do j=1,mxstars
       caldiams(j)=0.D0
      enddo
      ndiams=0
      do i=1,mxfk5
       read(3,51,end=52,err=52)fk,diameter
   51  format(2x,i4,2x,f5.2)
       do j=1,innums
        if(calstars(j).eq.fk)then
         caldiams(j)=dble(diameter)
         ndiams=ndiams+1
        endif
       enddo
      enddo
   52 close(3)
   50 fknam='      '
      write(outc,*)
      if(ndiams.gt.1)then
       write(outc,'(a,i2,a)')' Read ',ndiams,
     &         ' diameters from file: diameter.dat'
      else
       write(outc,'(a,i2,a)')' Read ',ndiams,
     &         ' diameter from file: diameter.dat'
      endif
      write(fknam(2:5),'(i4)')starno
      if(.not.calibs)then
       CALL READFK5( STARNO, STAR, RAM, DECM, PMRA, PMDEC, PAR, RV )
       write(outdat,'(a,i4,a,a)')'! Star is FK5 ',starno,', ',star
      endif
c
c ... Now read data and fill arrays obsvis2, viserr, etc...
      id=1
      ied=0
      nfilt=3
      idateold=0
      read(1,1)data
   30 read(data,*,end=112)rdatum,basel,lambda1,lambda3,lambda2
      datum=int(rdatum)
      ipos=index(data,'M')  ! read visibility multipliers
      if(ipos.ne.0.and..not.calibs)then
       write(outc,'(a,i6)')' *** WARNING: multipliers defined: ',datum
       read(data(ipos+1:256),*,end=112,err=112)
     &      mv(1),me(1),mv(3),me(3),mv(2),me(2)
      else
       do ig=1,nfilt
        mv(ig)=1.D0
        me(ig)=1.D0
       enddo
      endif

      do 100 n=1,100000
      read(1,1,end=111)data
      if(data.eq.'                                            ')goto 100
      read(data,*)time
      if(time.gt.24)goto 30 
      IF(INDEX(DATA,FKNAM).NE.0.or.CALIBS)THEN  ! to accelerate reading
      read(data,*,end=100)time,fk,d1,obsvis2(id,1),viserr(id,1),
     &d2,obsvis2(id,3),viserr(id,3),d3,obsvis2(id,2),viserr(id,2)
      do i=1,innums
      if(fk.eq.calstars(i))then
       stars(id)=fk
       baseline(id)=basel
       lambda0(id,1)=lambda1
       lambda0(id,2)=lambda2
       lambda0(id,3)=lambda3
       date(id)=datum
       hours(id)=time
c
c      We should flag bad data by setting vis to 1.0 and viserr to 2.0,
c      we also make sure that an error=0 causes no trouble!
       do j=1,mxfilt
        obsvis2(id,j)=obsvis2(id,j)*mv(j)
        viserr(id,j)=viserr(id,j)*me(j)
        if(obsvis2(id,j).gt.1.2D0.or.obsvis2(id,j).lt.-0.2D0)then
         obsvis2(id,j)=1.0D0
         viserr(id,j)=2.0D0
         ied=ied+1
        endif
        if(viserr(id,j).eq.0.D0)then
         viserr(id,j)=2.0D0
         obsvis2(id,j)=1.0D0
         ied=ied+1
        endif
c       if(fk.eq.193)then     ! try this for Capella
c        if(obsvis2(id,j).lt.0.05D0)viserr(id,j)=0.005D0
c        viserr(id,j)=max(viserr(id,j),0.005D0) 
c       endif
        if(calibs)then
         do l=1,iex
          if(fk.eq.exstar(l).and.lambda0(id,j).eq.exfilter(l))then
           obsvis2(id,j)=1.d0
           viserr(id,j)=2.d0
           ied=ied+1
          endif
         enddo
        endif
       enddo
c
c      Read user input for filter widths
       do j=1,mxfilt
        dlambda0(id,j)=0.D0
        do k=1,nnfilt
         if(lambda0(id,j).eq.filter(k))dlambda0(id,j)=k
        enddo
        if(dlambda0(id,j).eq.0)then
         nnfilt=nnfilt+1
         if(nnfilt.gt.mmxfilt)then
          write(outc,*)' Maximum number of filters exceeded !',
     &    '  Increase mmxfilt and recompile !'
          stop
         endif
         write(outc,'(a,f5.1,a,$)')' Found filter ',
     &   lambda0(id,j),' nm, please enter width [nm]: '
         read(inc,*)fwidth(nnfilt)
         filter(nnfilt)=lambda0(id,j)
         lamhi=filter(nnfilt)+fwidth(nnfilt)/2.d0
         lamlo=filter(nnfilt)-fwidth(nnfilt)/2.d0
         lamst=(lamhi-lamlo)/dble(nusteps-1)
         do l=1,nusteps
          ftcurve(nnfilt,l,1)=dble(l-1)*lamst+lamlo
          ftcurve(nnfilt,l,2)=100.d0 ! 100% Transmission
         enddo
         dlambda0(id,j)=nnfilt
        endif
       enddo
c
       id=id+1
       if(id.gt.mxrec)then
        write(outc,'(a,i4,a)')
     &    ' *** ERROR: more than ',mxrec,' records !'
        write(outc,*)'*** Increase MXREC and recompile!'
        stop
       endif
       if(datum.eq.idateold)then
        nscans=nscans+1
        if(nscans.gt.mxscn)then
         write(outc,'(a,i3,a,i6,a)')
     &   ' *** ERROR: more than ',mxscn,' records for night ',datum,' !'
         write(outc,*)'*** Increase MXSCN and recompile!'
         stop
        endif
       else
        idateold=datum
        nscans=1
       endif

      endif
      enddo
      ENDIF
  100 continue
  111 ndata=id-1
      write(outdat,'(a,i5)')'! Number of visibilities: ',ndata*nfilt
      write(outdat,'(a,i4,a)')'! Edited ',ied,' visibilities.'
      write(outc,'(a,i5)')' Number of visibilities read: ',ndata*nfilt
      write(outc,'(a,i4,a)')' Edited ',ied,' visibilities.'
      goto 113
  112 write(outc,*)'*** ERROR reading header line of cal file!'
      stop
  113 close(1)
      return
      end