      subroutine uvdist(modplot)
C
C Plots observed visibilities versus baseline length
C
C     C. A. Hummel          18 June 1993
C
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'
C
      integer*4 nsteps
      parameter (nsteps=40)
      real*4 blength(mxrec2),vis(mxrec2),modvis(nsteps)
      real*8 flux,arg,bessj1,diammas,dvis,lambda,ccolor
      integer*4 color(mxrec2),symbol(mxrec2)
      character header*80,device*4,plotfile*15
      logical modplot

      write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?'
      read(inc,2)device
    2 format(a4)
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotdist.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotdist.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotdist.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
c
      bmax=0.
      i=0
      do id=1,ndata
       do ig=1,nfilt
        if(viserr(id,ig).ne.2)then
         i=i+1
         blength(i)=real(sqrt(udata(id,ig)**2+vdata(id,ig)**2)/1.d6)
         bmax=max(blength(i),bmax)
         vis(i)=real(obsvis2(id,ig))
         if(ig.eq.1)then
          color(i)=2
          symbol(i)=4
         else if(ig.eq.2)then
          color(i)=3
          symbol(i)=6
         else if(ig.eq.3)then
          color(i)=11
          symbol(i)=7
         endif
        endif
       enddo
      enddo
c
c     Initialize plotting
c
      if(device.eq.'v'.or.device.eq.'V')then
       call pgbegin(0,plotfile,1,2)
      else
       call pgbegin(0,plotfile,1,1)
      endif
      call pgsch(1.5)
      call pgenv(0.,bmax,0.,1.,0,0)
      write(header,'(1x,a,i4,a)')'FK5 ',starno,' (all data)'
      if(device.eq.'v'.or.device.eq.'V')then
       call pglabel('Baseline length [millions of wavelengths]',
     &              'Visibility squared',' ')
      else
       call pglabel('Baseline length [millions of wavelengths]',
     &              'Visibility squared',header)
      endif
c
      call pgsch(1.0)
      do n=1,i
       call pgsci(color(n))
       call pgpoint(1,blength(n),vis(n),symbol(n))
      enddo
c
      if(.not.modplot)goto 111
      bstep=bmax/nsteps
      lambda=500.d0       ! for the time being, use this wavelength
      ccolor=ccoeff(3)*lambda*lambda+
     &       ccoeff(2)*lambda       +
     &       ccoeff(1)
     
      do i=1,nsteps
       modvis(i)=0.
       do ic=1,ncomp
        diammas=dcoeff(ic,3)*lambda*lambda+
     &          dcoeff(ic,2)*lambda       +
     &          dcoeff(ic,1)
        blength(i)=(i-1)*bstep
        arg=0.015231d0*dble(blength(i))*diammas
        if(arg.ne.0)then
         dvis=2.D0*abs(bessj1(arg))/arg
        else
         dvis=1.D0
        endif
        if(ic.ne.1)then
         flux=1.D0/(1.D0+10.D0**(-0.4D0*(delmag(ic,1)+ccolor)))
     &                 *(10.D0**(-0.4D0*(delmag(ic,1)+ccolor)))
        else
         flux=1.D0/(1.D0+10.D0**(-0.4D0*(delmag(2,1)+ccolor)))
        endif
        modvis(i)=modvis(i)+real(flux*dvis)
       enddo
       modvis(i)=modvis(i)**2
      enddo
      call pgsci(1)
      call pgline(nsteps,blength,modvis)
      write(header,'(a)')'(Uniform disc model calculated at 600 nm)'
      call pgsch(1.5)
      if(device.ne.'v'.and.device.eq.'V')
     &   call pgmtext('T',1.0,0.5,0.5,header)
c
  111 call pgend
c    
      return
      end

