c----------------------------------------------------------------------
c   Programme to plot fits to spectrum of the X-ray and Gamma-ray
c   backgrounds
c----------------------------------------------------------------------
      common/ c_plot/ xmin, xmax, ymin, ymax
      dimension xplot(400), yplot1(400), yplot2(400),
     > yplot3(400), yplot4(400), yplot5(400)
c
      xmin = 0.01
      xmax = 25.
      ymin = 0.7
      ymax = 1.2
      call plot_start
c
      nbin = 200
      dx = (xmax - xmin)/float(nbin-1)
c
      do i = 1, nbin
         x = xmin + float(i-1)*dx
         xplot(i) = x
         yplot1(i) = V(x, 0.)
         yplot2(i) = V(x, 2.*sqrt(3.))
         yplot3(i) = V(x, 3.75)
         yplot4(i) = V(x, 4.0)
         yplot5(i) = V(x, 4.5)
      end do
      call pgsch(1.4)
        
        h = 2.*sqrt(3.)
c       x = h**2/2.*(1. - sqrt(h*h-12))
        x = 6.
        call pgpoint(1, x, V(x, h), 17)
        h = 3.75
        x = 0.5*(h*h + h*sqrt(h*h-12.))
        call pgpoint(1, x, V(x, h), 17)
        h = 4.0
        x = 0.5*(h*h + h*sqrt(h*h-12.))
        call pgpoint(1, x, V(x, h), 17)
        h = 4.5
        x = 0.5*(h*h + h*sqrt(h*h-12.))
        call pgpoint(1, x, V(x, h), 17)
c
      call pgline(nbin, xplot, yplot1)
      call pgline(nbin, xplot, yplot2)
      call pgline(nbin, xplot, yplot3)
      call pgline(nbin, xplot, yplot4)
      call pgline(nbin, xplot, yplot5)
      call pgsch(0.8)
      call pgtext(10., 0.81, 'h\\u''\\d=0')
      call pgtext(6., 0.86, '2*(3)\\u1/2\\d')
      call pgtext(6., 0.91, '3.75')
      call pgtext(6., 0.95, '4.0')
      call pgtext(6., 1.06, '4.5')
      call pgend
      stop      
      end
c------------------------------------------------------------------
      function V(x, h)
      V = (1. - 2./x)*(1 + h*h/x/x)
      return
      end
c------------------------------------------------------------------
      subroutine plot_start
      common/ c_plot/ xmin, xmax, ymin, ymax
      call pgbegin(0, '?', 1, 1)
      call pgscf(2)
      print*, ' enter line weight'
      read*, ilw
      call pgslw(ilw)
c
      bord = 0.15
      dy = (1. - 2.*bord)
      xlo = bord
      xhi = 0.7*(1. - bord)
      ylo = bord
      yhi = (1. - bord)
c
      yhi = ylo + 0.8*dy
c
      xlo1 = xlo - 0.2*bord
      call pgsch(1.3)
      call pgvport(xlo1, xhi, ylo, yhi)
      call pglabel('x = 2r/r\\ds\\u', 'V(x)', ' ')
      call pgvport(xlo, xhi, ylo, yhi)
      call pgwindow(xmin, xmax, ymin, ymax)
      call pgbox('bncts', 5., 5, 'bcnvts', 0.1, 5)
      call pgsch(1.0)
      return
      end