      subroutine plotvis(modplot,modvis2,modphas,mode,device)
C
C Plots observed visibilities and model
C
C     C. A. Hummel          16 September 1991
C
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'
C
      real*8 modvis2(mxrec,mxfilt),modphas(mxrec,mxfilt)
      real*8 modv2(mxpl,mxfilt),modph(mxpl,mxfilt),mhours(mxpl)
c     real*8 modv3(mxpl,mxfilt),modph3(mxpl,mxfilt)
      real*8 lambda(mxpl,mxfilt),dlambda(mxpl,mxfilt)
      real*8 u(mxpl,mxfilt),v(mxpl,mxfilt)
      real*8 timstep,chisq,mseprad,mparad
      real*4 y(mxpl),t(mxpl),vis2max(mxfilt),vis2min(mxfilt)
      real*4 timmax,timmin
      integer*4 firstid,lastid,mode
      integer*4 mdate(mxpl),mbaseline(mxpl),mstars(mxpl)
      logical modplot
      character header*80,device*4,plotfile*15

      if(mode.eq.0)then
       write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
       read(inc,2)device
      else
       modplot=.true.
      endif
    2 format(a4)
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotvis.lpm'//'/lpm'
       iy=4
       csize=2.5
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotvis.ps'//'/ps'
       iy=4
       csize=2.5
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotvis.ps'//'/vps'
       iy=5
       csize=2.0
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
       iy=2
       csize=2.5
      else
       return
      endif
c
c     Initialize plotting
c
      if(mode.eq.0)call pgbegin(0,plotfile,3,iy)
      call pgsch(csize)
      firstid=1
      lastid=1
      idate=1
      do ig=1,nfilt
       vis2max(ig)=0
       vis2min(ig)=0
      enddo
      timmin=real(hours(1))
      timmax=timmin
c........Interim: plot second set of magnitude differences
c        dmag1=delmag(2,1)
c        dmag2=delmag(2,2)
c        dmag3=delmag(2,3)
c        read(16,*)dnag1,dnag2,dnag3
c
c     Go through the whole data file and ...
c 
      DO ID=1,NDATA
       do ig=1,nfilt
        if(viserr(id,ig).ne.2)then
        vis2max(ig)=max(vis2max(ig),real(obsvis2(id,ig)))
        vis2min(ig)=min(vis2min(ig),real(obsvis2(id,ig)))
        endif
       enddo
       timmax=max(timmax,real(hours(id)))
       timmin=min(timmin,real(hours(id)))
c
c      ... output plot, when date or baseline changes
c
       if(date(id+1).ne.date(id).or.baseline(id+1).ne.baseline(id)
     &                          .or.id+1.gt.ndata)then
        timmin=timmin-0.5
        timmax=timmax+0.5

        if(modplot)then  ! Calculate model if requested
         timstep=0.2D0
         do k=1,mxpl
          mstars(k)=stars(id)
          mdate(k)=date(id)
          mbaseline(k)=baseline(id)
          mhours(k)=dble(timmin)+timstep*(k-1)
          do ig=1,nfilt
           lambda(k,ig)=lambda0(id,ig)
           dlambda(k,ig)=dlambda0(id,ig)
          enddo
          if(mhours(k).gt.dble(timmax))goto 40
         enddo
   40    ndat=min(k,mxpl)
         call calcuv(ndat,nfilt,mxpl,mxfilt,mstars,mdate,mbaseline,
     &               mhours,lambda,u,v)
         if(.not.calib)then
          imode=0
          call binavis(ndat,nfilt,mxpl,mxfilt,mdate,mbaseline,mhours,
     &                 lambda,dlambda,modv2,modph,u,v,
     &                 mseprad,mparad,imode,ftcurve)
c         delmag(2,1)=dnag1
c         delmag(2,2)=dnag2
c         delmag(2,3)=dnag3
c         call cadjmod
c         call binavis(ndat,nfilt,mxpl,mxfilt,mdate,mbaseline,mhours,
c    &                 lambda,dlambda,modv3,modph3,u,v,
c    &                 mseprad,mparad,imode,ftcurve)
         else
          call calvis(ndat,nfilt,mxpl,mxfilt,mdate,mbaseline,mhours,
     &                lambda,modv2,u,v,mstars,calstars,caldiams)
         endif
c        delmag(2,1)=dmag1
c        delmag(2,2)=dmag2
c        delmag(2,3)=dmag3
c        call cadjmod
        endif

       if(starno.eq.193.and.(date(id).eq.920929.or.date(id).eq.920930))
     &  then
        vis2max(1)=vis2max(3)
        vis2max(2)=vis2max(3)
        if(date(id).eq.920929)vis2max(2)=vis2max(2)/2.
       endif
c
c      loop over nfilt plots
c
       do ig=1,nfilt
        call pgsci(5)
        call pgsls(1)
        vr=vis2max(ig)-vis2min(ig)
        vis2max(ig)=vis2max(ig)+vr/8.
        call pgenv(timmin,timmax,vis2min(ig),vis2max(ig),0,1)
        if(ig.eq.1)then
         write(header,3)starno,int(lambda0(id,ig))
    3    format('FK5 ',i4,3x,i3,'nm')
         if(calib)write(header(1:8),'(a)')'calibr. '
        elseif(ig.eq.2)then
         write(header,4)date(id),int(lambda0(id,ig))
    4    format(2x,i6,3x,i3,'nm')
        elseif(ig.eq.3)then
         write(header,5)baseline(id),int(lambda0(id,ig))
    5    format(3x,'Bl ',i2,3x,i3,'nm')
        endif
        csize=csize+1
        call pgsch(csize)
        if(device.eq.'v'.or.device.eq.'V')then
c        write(header,8)int(lambda0(id,ig))
c   8    format(i3,'nm')
         call pgmtext('T',1.0,0.5,0.5,header)
         call pglabel('UT','squared visibility',' ')
        else
         call pglabel('','',header)
        endif
        csize=csize-1
        call pgsch(csize)
        if(ig.eq.1)call pgsci(2)
        if(ig.eq.2)call pgsci(3)
        if(ig.eq.3)call pgsci(11)
        chisq=0.D0
        nchisq=0
        do n=1,lastid-firstid+1
         yp=real(obsvis2(n+firstid-1,ig))
         tp=real(hours(n+firstid-1))
         if(viserr(n+firstid-1,ig).ne.2)then
          call pgpoint(1,tp,yp,1)
          obsvis2lo=yp-real(viserr(n+firstid-1,ig))
          obsvis2hi=yp+real(viserr(n+firstid-1,ig))
          call pgerry(1,tp,obsvis2lo,obsvis2hi,1.0)
          chisq=chisq
     &         +((obsvis2(n+firstid-1,ig)-modvis2(n+firstid-1,ig))
     &          /viserr(n+firstid-1,ig))**2
          nchisq=nchisq+1
         endif
        enddo
        header='                                        '//
     &         '                                        '
        if(nchisq.ne.0.and.modplot)then
         call pgsci(6)
         if(vis2rat(idate,ig).eq.0)then
          write(header,6)real(chisq)/real(nchisq)
    6     format('\gx\u2\d/\gn = ',f6.2) 
         else
          write(header,7)real(chisq)/real(nchisq),vis2rat(idate,ig)
    7     format('\gx\u2\d/\gn = ',f6.2,', O/C = ',f4.2)
         endif
         call pgmtext('T',1.0,0.5,0.5,header)
        endif
c
c       Plot model visibilities here and now
        if(modplot)then
         do k=1,ndat
          y(k)=real(modv2(k,ig))
          t(k)=real(mhours(k))
         enddo
         call pgsci(6)
         call pgsls(1)
         call pgline(ndat,t,y)
c        do k=1,ndat
c         y(k)=real(modv3(k,ig))
c         t(k)=real(mhours(k))
c        enddo
c        call pgsls(4)
c        call pgline(ndat,t,y)
        endif 
c
       enddo      ! End of ig-loop
c
c      Now prepare for the next date
       idate=idate+1
       do ig=1,nfilt
        vis2max(ig)=0
        vis2min(ig)=0
       enddo
       timmin=real(hours(id+1))
       timmax=timmin
       firstid=lastid+1
       lastid=firstid-1
c
c     continue here, if date did not change
      endif
      lastid=lastid+1
      ENDDO
      if(mode.eq.0)call pgend
      return
      end
