c
c
C       NOTE;  COMPILE THIS WITH NORMAL F77 COMPILER!
C              DO NOT USE THE R8 OPTION!
C
c
c*************************************************************************
c
      subroutine SCPLT1 ( x, y, m, isymbl, istyle, 
c
     *                                   title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
c  subroutine to plot  y  vs.  x
c
c       m = no. data pts.
c
c  isymbl = symbol code   isymbl =  0 : no symbols
c                                =  1 : circles
c                                =  2 : squares
c                                =  3 : triangles
c                                =  4 : deltas
c                                =  5 : diamonds
c                                ..................
c                                = 23 : (see vg manual)
c
c  istyle = line code     istyle =  o : no line
c                                =  1 : solid line
c                                =  2 : dashed line
c                                ...................
c                                = 11 : (see vg manual)
c
c   title = character string containing title for plot
c
c  xlabel = character string containing label for x-axis
c
c  ylabel = character string containing label for y-axis
c
c------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  x(m), y(m)
c
      real  xpl(npl), ypl(npl)
c
c-------------------------------------------------------------------------
c
      if ( m .gt. npl ) then
c
                  print *, ' in SCPLT1, npl =', npl, ' is too small'
c                ++++++++
                  return
c                ++++++++
      endif
c
      do 10 i=1,m
         xpl(i) = x(i)
         ypl(i) = y(i)
   10 continue
c
 100  continue
c
         call axcode ( 0, 0, 2 )
c
         call howplt ( isymbl, istyle, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call curv ( m, xpl, ypl )
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c***************************************************************************
c
      subroutine SCPLT2 ( xp, yp, np, xc, yc, nc, isymbl, istyle, 
c
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
c  subroutine to plot  ( yp vs. xp )  and  ( yc vs. xc )  on same graph
c
c      xp = independent variable to be plotted as discrete points
c
c      yp = dependent variable to be plotted as discrete points
c
c      np = no. ( xp , yp ) data points
c
c      xc = indepencent variable to be plotted as curve
c
c      yc = dependent variable to be plotted as curve
c
c      nc = no. data pts.
c
c  isymbl = symbol code   isymbl =  0 : no symbols (print as solid line)
c           ( for  yp )          =  1 : circles
c                                =  2 : squares
c                                =  3 : triangles
c                                =  4 : deltas
c                                =  5 : diamonds
c                                ..................
c                                = 23 : (see vg manual)
c
c  istyle = line code     istyle =  o : no line
c           ( for  yc )          =  1 : solid line
c                                =  2 : dashed line
c                                ...................
c                                = 11 : (see vg manual)
c
c   title = character string containing title for plot
c
c  xlabel = character string containing label for x-axis
c
c  ylabel = character string containing label for y-axis
c
c------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np), yp(np), xc(nc), yc(nc)
c
      dimension  xppl(npl), yppl(npl), xcpl(npl), ycpl(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
c
                  print *, ' in SCPLT2, npl =', npl, ' is too small'
c                ++++++++
                  return
c                ++++++++
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
 100  continue
c
         call axcode ( 0, 0, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         if ( isymbl .ne. 0 ) then
                              call howplt ( isymbl, 0, 2 )
                      else 
                           call howplt ( 0, 1, 2 )
         end if
c
         call curv ( np, xppl, yppl )
c
         if ( istyle .eq. 0 ) then
                              call howplt ( isymbl+1 , 0, 2 )
                     else
                          call howplt ( 0, istyle, 2 )
         end if
c
         call curv ( nc, xcpl, ycpl )
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c*************************************************************************
c
c
      subroutine PLOT3 ( xp,  yp,  np,  ipsym,  ipsty,
     *                   xc,  yc,  nc,  icsym,  icsty, 
     *                   xc2, yc2, nc2, ic2sym, ic2sty,
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np), yp(np), xc(nc), yc(nc), 
     1                                xc2(nc2), yc2(nc2)
c
      dimension  xppl(npl), yppl(npl), xcpl(npl), ycpl(npl),
     1                                xcpl2(npl), ycpl2(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
                  print *, ' in PLOT3, npl =', npl, ' is too small'
c         ********
           return
c         ********
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
      do 30 i=1,nc2
          xcpl2(i) = xc2(i)
          ycpl2(i) = yc2(i)
   30 continue
c
 100  continue
c
         call axcode ( 0, 0, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call howplt ( ipsym, ipsty, 3 )
         call curv ( np, xppl, yppl )
c
         call howplt ( icsym, icsty, 3 )
         call curv ( nc, xcpl, ycpl )
c
         call howplt ( ic2sym, ic2sty, 3 )
         call curv ( nc2, xcpl2, ycpl2 )         
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c*************************************************************************
c
c
      subroutine PLOT4 ( xp,  yp,  np,  ipsym,  ipsty,
     *                   xc,  yc,  nc,  icsym,  icsty, 
     *                   xc2, yc2, nc2, ic2sym, ic2sty,
     *                   xc3, yc3, nc3, ic3sym, ic3sty,
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np), yp(np), xc(nc), yc(nc), 
     1                                xc2(nc2), yc2(nc2),
     2                                xc3(nc3), yc3(nc3)
c
      dimension  xppl(npl), yppl(npl), xcpl(npl), ycpl(npl),
     1                                xcpl2(npl), ycpl2(npl),
     2                                xcpl3(npl), ycpl3(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
                  print *, ' in PLOT3, npl =', npl, ' is too small'
c         ********
           return
c         ********
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
      do 30 i=1,nc2
          xcpl2(i) = xc2(i)
          ycpl2(i) = yc2(i)
   30 continue
c
      do 40 i=1,nc3
          xcpl3(i) = xc3(i)
          ycpl3(i) = yc3(i)
   40 continue
c
 100  continue
c
         call axcode ( 0, 0, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call howplt ( ipsym, ipsty, 3 )
         call curv ( np, xppl, yppl )
c
         call howplt ( icsym, icsty, 3 )
         call curv ( nc, xcpl, ycpl )
c
         call howplt ( ic2sym, ic2sty, 3 )
         call curv ( nc2, xcpl2, ycpl2 )         
c
         call howplt ( ic3sym, ic3sty, 3 )
         call curv ( nc3, xcpl3, ycpl3 )         
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c*************************************************************************
c
c
      subroutine PLOT5 ( xp,  yp,  np,  ipsym,  ipsty,
     *                   xc,  yc,  nc,  icsym,  icsty, 
     *                   xc2, yc2, nc2, ic2sym, ic2sty,
     *                   xc3, yc3, nc3, ic3sym, ic3sty,
     *                   xc4, yc4, nc4, ic4sym, ic4sty,
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np), yp(np), xc(nc), yc(nc), 
     1                                xc2(nc2), yc2(nc2),
     2                                xc3(nc3), yc3(nc3),
     3                                xc4(nc4), yc4(nc4)
c
      dimension  xppl(npl), yppl(npl), xcpl(npl), ycpl(npl),
     1                                xcpl2(npl), ycpl2(npl),
     2                                xcpl3(npl), ycpl3(npl),
     3                                xcpl4(npl), ycpl4(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
                  print *, ' in PLOT3, npl =', npl, ' is too small'
c         ********
           return
c         ********
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
      do 30 i=1,nc2
          xcpl2(i) = xc2(i)
          ycpl2(i) = yc2(i)
   30 continue
c
      do 40 i=1,nc3
          xcpl3(i) = xc3(i)
          ycpl3(i) = yc3(i)
   40 continue
c
      do 50 i=1,nc4
          xcpl4(i) = xc4(i)
          ycpl4(i) = yc4(i)
   50 continue
c
 100  continue
c
         call axcode ( 0, 0, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call howplt ( ipsym, ipsty, 3 )
         call curv ( np, xppl, yppl )
c
         call howplt ( icsym, icsty, 3 )
         call curv ( nc, xcpl, ycpl )
c
         call howplt ( ic2sym, ic2sty, 3 )
         call curv ( nc2, xcpl2, ycpl2 )         
c
         call howplt ( ic3sym, ic3sty, 3 )
         call curv ( nc3, xcpl3, ycpl3 )         
c
         call howplt ( ic4sym, ic4sty, 3 )
         call curv ( nc4, xcpl4, ycpl4 )         
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c****************************************************************************
c
c
      subroutine HSTPLT ( vgx, vgy, nbox, barw,
c
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  barw, vgx(nbox), vgy(nbox)
c
      real  barwid, vgxhst(npl), vgyhst(npl)
c
c-----------------------------------------------------------------------
c
      barwid = barw
      do 10 i=1,nbox
         vgxhst(i) = vgx(i)
         vgyhst(i) = vgy(i)
   10 continue
c
  100 continue
c
         call axcode ( 0, 0, 3 )
c
         call howplt ( 0, 1, 3 )
c
         call sidtex ( title, 3, xlabel, 3, ylabel, 3, ' ', 3 )
c
         call hist ( 30, vgxhst, vgyhst, barwid )
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c*************************************************************************
c
c
      subroutine LNLGPL2 ( xp,  yp,  np,  ipsym,  ipsty,
     *                     xc,  yc,  nc,  icsym,  icsty, 
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np), yp(np), xc(nc), yc(nc) 
c
      dimension  xppl(npl), yppl(npl), xcpl(npl), ycpl(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
                  print *, ' in LNLGPL, npl =', npl, ' is too small'
c         ********
           return
c         ********
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
 100  continue
c
         call linlog ( 0, 1, 0 )
         call axcode ( 0, 0, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call howplt ( ipsym, ipsty, 3 )
         call curv ( np, xppl, yppl )
c
         call howplt ( icsym, icsty, 3 )
         call curv ( nc, xcpl, ycpl )
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c*************************************************************************
c
c
      subroutine LNLGPL3 ( xp,  yp,  np,  ipsym,  ipsty,
     *                     xc,  yc,  nc,  icsym,  icsty, 
     *                     xc2, yc2, nc2, ic2sym, ic2sty,
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np), yp(np), xc(nc), yc(nc), 
     1                                xc2(nc2), yc2(nc2)
c
      dimension  xppl(npl), yppl(npl), xcpl(npl), ycpl(npl),
     1                                xcpl2(npl), ycpl2(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
                  print *, ' in LNLGPL, npl =', npl, ' is too small'
c         ********
           return
c         ********
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
      do 30 i=1,nc2
          xcpl2(i) = xc2(i)
          ycpl2(i) = yc2(i)
   30 continue
c
 100  continue
c
         call linlog ( 0, 1, 0 )
         call axcode ( 0, 0, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call legpos ( 0.8, 0.8 )
c
         call howplt ( ipsym, ipsty, 3 )
         call legnd ( '|s_j(A)' )
         call curv ( np, xppl, yppl )
c
         call howplt ( icsym, icsty, 3 )
         call legnd ( '||U^Tb||_j' )
         call curv ( nc, xcpl, ycpl )
c
         call howplt ( ic2sym, ic2sty, 3 )
         call curv ( nc2, xcpl2, ycpl2 )         
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c
c*************************************************************************
c
c
      subroutine LNLGPL4 ( xp,  yp,  np,  ipsym,  ipsty,
     *                     xc,  yc,  nc,  icsym,  icsty, 
     *                     xc2, yc2, nc2, ic2sym, ic2sty,
     *                     xc3, yc3, nc3, ic3sym, ic3sty,
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np),   yp(np),   xc(nc),   yc(nc), 
     1                  xc2(nc2), yc2(nc2), xc3(nc3), yc3(nc3)
c
      dimension  xppl(npl),  yppl(npl),  xcpl(npl),  ycpl(npl),
     1           xcpl2(npl), ycpl2(npl), xcpl3(npl), ycpl3(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
                  print *, ' in LNLGPL, npl =', npl, ' is too small'
c         ********
           return
c         ********
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
      do 30 i=1,nc2
          xcpl2(i) = xc2(i)
          ycpl2(i) = yc2(i)
   30 continue
c
      do 40 i=1,nc3
          xcpl3(i) = xc3(i)
          ycpl3(i) = yc3(i)
   40 continue
c
 100  continue
c
         call linlog ( 0, 1, 0 )
         call axcode ( 0, 0, 2 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call legpos ( 0.8, 0.8 )
c
         call howplt ( ipsym, ipsty, 3 )
         call legnd ( '|s_j(A)' )
         call curv ( np, xppl, yppl )
c
         call howplt ( icsym, icsty, 3 )
         call legnd ( '||U^Tb||_j' )
         call curv ( nc, xcpl, ycpl )
c
         call howplt ( ic2sym, ic2sty, 3 )
         call curv ( nc2, xcpl2, ycpl2 ) 
c
         call howplt ( ic3sym, ic3sty, 3 )
         call curv ( nc3, xcpl3, ycpl3 )
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c***************************************************************************
c
      subroutine LGLGPL2 ( xp, yp, np, xc, yc, nc, isymbl, istyle, 
c
     *                                      title, xlabel, ylabel )
c
c-------------------------------------------------------------------------
c
c  subroutine to plot  ( yp vs. xp )  and  ( yc vs. xc )  on same graph
c
c      xp = independent variable to be plotted as discrete points
c
c      yp = dependent variable to be plotted as discrete points
c
c      np = no. ( xp , yp ) data points
c
c      xc = indepencent variable to be plotted as curve
c
c      yc = dependent variable to be plotted as curve
c
c      nc = no. data pts.
c
c  isymbl = symbol code   isymbl =  0 : no symbols (print as solid line)
c           ( for  yp )          =  1 : circles
c                                =  2 : squares
c                                =  3 : triangles
c                                =  4 : deltas
c                                =  5 : diamonds
c                                ..................
c                                = 23 : (see vg manual)
c
c  istyle = line code     istyle =  o : no line
c           ( for  yc )          =  1 : solid line
c                                =  2 : dashed line
c                                ...................
c                                = 11 : (see vg manual)
c
c   title = character string containing title for plot
c
c  xlabel = character string containing label for x-axis
c
c  ylabel = character string containing label for y-axis
c
c------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title, xlabel, ylabel
c
      double precision  xp(np), yp(np), xc(nc), yc(nc)
c
      dimension  xppl(npl), yppl(npl), xcpl(npl), ycpl(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np .gt. npl) .or. (nc .gt. npl) ) then
c
                  print *, ' in SCPLT2, npl =', npl, ' is too small'
c                ++++++++
                  return
c                ++++++++
      endif
c
      do 10 i=1,np
          xppl(i) = xp(i)
          yppl(i) = yp(i)
   10 continue
c
      do 20 i=1,nc
          xcpl(i) = xc(i)
          ycpl(i) = yc(i)
   20 continue
c
 100  continue
c
         call axcode ( 0, 0, 2 )
c
         call linlog ( 1, 1, 0 )
c
         call sidtex ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         if ( isymbl .ne. 0 ) then
                              call howplt ( isymbl, 0, 2 )
                      else 
                           call howplt ( 0, 1, 2 )
         end if
c
         call curv ( np, xppl, yppl )
c
         if ( istyle .eq. 0 ) then
                              call howplt ( isymbl+1 , 0, 2 )
                     else
                          call howplt ( 0, istyle, 2 )
         end if
c
         call curv ( nc, xcpl, ycpl )
c
         call vg
c
      if ( loopin() .eq. 1 ) go to 100
c
      call newpag
c
      return
c
      end
c
c
c*******************************************************************************
c
c
      subroutine TWOPLT ( xp1, yp1, np1, xc1, yc1, nc1, isym1, isty1,
c
     1                                   title1, xlabl1, ylabl1,
c
     2                    xp2, yp2, np2, xc2, yc2, nc2, isym2, isty2,
c
     3                                   title2, xlabl2, ylabl2, ts )
c
c----  Subroutine to make two VG plots (of two curves each) on same page
c
c-------------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      logical  ts
c
      character*(*)  title1, xlabl1, ylabl1,
     1               title2, xlabl2, ylabl2
c
      double precision  xp1(np1), yp1(np1), xc1(nc1), yc1(nc1),
     1                  xp2(np2), yp2(np2), xc2(nc2), yc2(nc2)
c
      dimension  xppl1(npl), yppl1(npl), xcpl1(npl), ycpl1(npl),
     1           xppl2(npl), yppl2(npl), xcpl2(npl), ycpl2(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np1 .gt. npl) .or. (nc1 .gt. npl)
     1                    .or.
     2     (np2 .gt. npl) .or. (nc2 .gt. npl) ) then
c
                  print *, ' in TWOPLT, npl =', npl, ' is too small'
c                ++++++++
                  return
c                ++++++++
      endif
c
      do 10 i=1,np1
          xppl1(i) = xp1(i)
          yppl1(i) = yp1(i)
   10 continue
c
      do 20 i=1,nc1
          xcpl1(i) = xc1(i)
          ycpl1(i) = yc1(i)
   20 continue
c
      do 11 i=1,np2
          xppl2(i) = xp2(i)
          yppl2(i) = yp2(i)
   11 continue
c
      do 21 i=1,nc2
          xcpl2(i) = xc2(i)
          ycpl2(i) = yc2(i)
   21 continue
c
c-------------------------------------------------------------------------
c
 100  continue
c
      if ( ts ) then
                     call PLAC ( 1.0, 0.52, 0.0, 1.0 )
                else 
                     call PLAC ( 0.8, 0.15 , 0.0, 0.48 )
      end if
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title1, 2, xlabl1, 2, ylabl1, 2, ' ', 2 )
c
         if ( isym1 .ne. 0 ) then
                              call HOWPLT ( isym1, 0, 2 )
                      else 
                           call HOWPLT ( 0, 1, 2 )
         end if
c
         call CURV ( np1, xppl1, yppl1 )
c
         if ( isty1 .eq. 0 ) then
                              call HOWPLT ( isym1+1 , 0, 2 )
                     else
                          call HOWPLT ( 0, isty1, 2 )
         end if
c
         call CURV ( nc1, xcpl1, ycpl1 )
c
         call VG
c
      if ( ts ) then
                     call PLAC ( 0.48, 0.0, 0.0, 1.0 )
                else
                     call PLAC ( 0.8, 0.15, 0.52, 1.0)
      end if
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title2, 2, xlabl2, 2, ylabl2, 2, ' ', 2 )
c
         if ( isym2 .ne. 0 ) then
                              call HOWPLT ( isym2, 0, 2 )
                      else 
                           call HOWPLT ( 0, 1, 2 )
         end if
c
         call CURV ( np2, xppl2, yppl2 )
c
         if ( isty2 .eq. 0 ) then
                              call HOWPLT ( isym2+1 , 0, 2 )
                     else
                          call HOWPLT ( 0, isty2, 2 )
         end if
c
         call CURV ( nc2, xcpl2, ycpl2 )
c
         call VG
c
      if ( loopin() .eq. 1 ) go to 100
c
      call NEWPAG
c
      return
c
      end
c
c*******************************************************************************
c
      subroutine THRPLT ( xp1, yp1, np1, xc1, yc1, nc1, isym1, isty1,
c
     1                                   title1, xlabl1, ylabl1,
c
     2                    xp2, yp2, np2, xc2, yc2, nc2, isym2, isty2,
c
     3                                   title2, xlabl2, ylabl2,   
c
     4                     x, y, m, isymbl, istyle, 
c
     5                                   title, xlabel, ylabel,  ts )
c
c  Subroutine to make 3 VG plots (2 w/ 2 curves & 1 w/ 1 curve) on same page
c
c-------------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      logical  ts
c
      character*(*)  title1, xlabl1, ylabl1,
     1               title2, xlabl2, ylabl2,
     2               title,  xlabel, ylabel
c
      double precision  xp1(np1), yp1(np1), xc1(nc1), yc1(nc1),
     1                  xp2(np2), yp2(np2), xc2(nc2), yc2(nc2),
     2                  x(m),     y(m)
c
      dimension  xppl1(npl), yppl1(npl), xcpl1(npl), ycpl1(npl),
     1           xppl2(npl), yppl2(npl), xcpl2(npl), ycpl2(npl),
     2           xpl(npl),     ypl(npl)
c
c-------------------------------------------------------------------------
c
      if ( (np1 .gt. npl) .or. (nc1 .gt. npl)
     1                    .or.
     2     (np2 .gt. npl) .or. (nc2 .gt. npl) 
     3                    .or. 
     4                         ( m  .gt. npl) ) then
c
                  print *, ' in TWOPLT, npl =', npl, ' is too small'
c                ++++++++
                  return
c                ++++++++
      endif
c
      do 10 i=1,np1
          xppl1(i) = xp1(i)
          yppl1(i) = yp1(i)
   10 continue
c
      do 20 i=1,nc1
          xcpl1(i) = xc1(i)
          ycpl1(i) = yc1(i)
   20 continue
c
      do 11 i=1,np2
          xppl2(i) = xp2(i)
          yppl2(i) = yp2(i)
   11 continue
c
      do 21 i=1,nc2
          xcpl2(i) = xc2(i)
          ycpl2(i) = yc2(i)
   21 continue
c
      do 31 i=1,m
         xpl(i) = x(i)
         ypl(i) = y(i)
   31 continue
c
c-------------------------------------------------------------------------
c
 100  continue
c
      if ( ts ) then
                     call PLAC ( 1.0, 0.52, 0.0, 1.0 )
                else 
                     call PLAC ( 1.0, 0.52 , 0.0, 0.48 )
      end if
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title1, 2, xlabl1, 2, ylabl1, 2, ' ', 2 )
c
         if ( isym1 .ne. 0 ) then
                              call HOWPLT ( isym1, 0, 2 )
                      else 
                           call HOWPLT ( 0, 1, 2 )
         end if
c
         call CURV ( np1, xppl1, yppl1 )
c
         if ( isty1 .eq. 0 ) then
                              call HOWPLT ( isym1+1 , 0, 2 )
                     else
                          call HOWPLT ( 0, isty1, 2 )
         end if
c
         call CURV ( nc1, xcpl1, ycpl1 )
c
         call VG
c
      call PLAC ( 0.48, 0.0, 0.0, 0.48 )
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title2, 2, xlabl2, 2, ylabl2, 2, ' ', 2 )
c
         if ( isym2 .ne. 0 ) then
                              call HOWPLT ( isym2, 0, 2 )
                      else 
                           call HOWPLT ( 0, 1, 2 )
         end if
c
         call CURV ( np2, xppl2, yppl2 )
c
         if ( isty2 .eq. 0 ) then
                              call HOWPLT ( isym2+1 , 0, 2 )
                     else
                          call HOWPLT ( 0, isty2, 2 )
         end if
c
         call CURV ( nc2, xcpl2, ycpl2 )
c
         call VG
c
      call PLAC ( 0.48, 0.0, 0.52, 1.0 )
c
         call AXCODE ( 0, 0, 2 )
c
         call HOWPLT ( isymbl, istyle, 2 )
c
         call SIDTEX ( title, 2, xlabel, 2, ylabel, 2, ' ', 2 )
c
         call CURV ( m, xpl, ypl )
c
         call VG
c
      if ( loopin() .eq. 1 ) go to 100
c
      call NEWPAG
c
      return
c
      end
c
c*******************************************************************************
c
      subroutine FOURPL ( xp1, yp1, np1, isym1, isty1,
c
     1                    sumytr, yp1tr, isym12, isty12,
c
     1                                   title1, xlabl1, ylabl1,
c
     2                    xp2,   yp2, np2,  xc2, yc2, nc2, 
c
     2                    isymp2, istyp2,  isymc2, istyc2,
c
     2                                   title2, xlabl2, ylabl2,   
c
     3                    xp3, yp3, np3, isym3, isty3, 
c
     3                                   title3, xlabl3, ylabl3,
c
     4                    xp4, yp4, np4,  xc41, yc41, nc41, 
c
     4                                    xc42, yc42, nc42,
c
     4                 isymp4, istyp4, isym41, isty41, isym42, isty42,
c
     4                                   title4, xlabl4, ylabl4,
c
     4                                   legnd4  )
c
c-------------------------------------------------------------------------------
c
c  Subroutine to make 34 VG plots on same page
c        (#1 w/ 1 curve, #2 w/ 2 curves, #3 w/ 1 curve, &  #4 w/ 3 curves) 
c
c-------------------------------------------------------------------------------
c
      parameter ( npl = 40000 )
c
      character*(*)  title1, xlabl1, ylabl1,
     1               title2, xlabl2, ylabl2,
     2               title3, xlabl3, ylabl3,
     3               title4, xlabl4, ylabl4, legnd4
c
      double precision  xp1(np1), yp1(np1), yp1tr(np1),
     1                  xp2(np2), yp2(np2), xc2(nc2),   yc2(nc2),
     2                  xp3(np3), yp3(np3),
     3                  xp4(np4), yp4(np4), xc41(nc41), yc41(nc41),
     4                                      xc42(nc42), yc42(nc42)
c
      dimension  xppl1(npl), yppl1(npl), ytrpl1(npl),
     1           xppl2(npl), yppl2(npl), xcpl2(npl), ycpl2(npl),
     2                  xppl3(npl), yppl3(npl),
     3                  xppl4(npl), yppl4(npl),
     4                                      xcpl41(npl), ycpl41(npl),  
     4                                      xcpl42(npl), ycpl42(npl)          
c
c-------------------------------------------------------------------------------
c
      do i=1,np1
         xppl1(i) = xp1(i)
         yppl1(i) = yp1(i)
         ytrpl1(i) = yp1tr(i)
      end do
c
      do i=1,np2
         xppl2(i) = xp2(i)
         yppl2(i) = yp2(i)
      end do
c
      do i=1,nc2
         xcpl2(i) = xc2(i)
         ycpl2(i) = yc2(i)
      end do
c
      do i=1,np3
         xppl3(i) = xp3(i)
         yppl3(i) = yp3(i)
      end do
c
      do i=1,np4
         xppl4(i) = xp4(i)
         yppl4(i) = yp4(i)
      end do
c
      do i=1,nc41
         xcpl41(i) = xc41(i)
         ycpl41(i) = yc41(i)
      end do
c
      do i=1,nc42
         xcpl42(i) = xc42(i)
         ycpl42(i) = yc42(i)
      end do
c
c-------------------------------------------------------------------------
c
c--------------
 100  continue
c--------------
c
      call PLAC ( 1.0, 0.52 , 0.0, 0.48 )
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title1, 2, xlabl1, 2, ylabl1, 2, ' ', 2 )
c
         call HOWPLT ( isym1, isty1, 2 )
c
         call CURV ( np1, xppl1, yppl1 )
c
         if ( sumytr .ne. 0.0 ) then
c
                                call HOWPLT ( isym12, isty12, 2 )
c
                                call CURV ( np1, xppl1, ytrpl1 )
c 
         end if
c
         call VG
c
c-------------------------------------------------------------------------------
c
      call PLAC ( 1.0, 0.52, 0.52, 1.0 )
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title2, 2, xlabl2, 2, ylabl2, 2, ' ', 2 )
c
         call HOWPLT ( isymp2, istyp2, 2 )
c
         call CURV ( np2, xppl2, yppl2 )
c
         call HOWPLT ( isymc2 , istyc2, 2 )
c
         call CURV ( nc2, xcpl2, ycpl2 )
c
         call VG
c
c-------------------------------------------------------------------------------
c
      call PLAC ( 0.48, 0.0, 0.0, 0.48 )
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title3, 2, xlabl3, 2, ylabl3, 2, ' ', 2 )
c
         call HOWPLT ( isym3, isty3, 2 )
c
         call CURV ( np3, xppl3, yppl3 )
c
         call VG
c
c-------------------------------------------------------------------------------
c
      call PLAC ( 0.48, 0.0, 0.52, 1.0 )
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( title4, 2, xlabl4, 2, ylabl4, 2, ' ', 2 )
c
         call HOWPLT ( isymp4, istyp4, 2 )
c
         call LEGND ( legnd4 )
c
         call CURV ( np4, xppl4, yppl4 )
c
         call HOWPLT ( isym41, isty41, 2 )
c
         call CURV ( nc41, xcpl41, ycpl41 )
c
         call HOWPLT ( isym42, isty42, 2 )
c
         call CURV ( nc42, xcpl42, ycpl42 )
c
         call VG
c
c-------------------------------------------------------------------------------
c
c                           -----------
      if ( loopin() .eq. 1 ) go to 100
c                           -----------
c
      call NEWPAG
c
      return
c
      end
c
c
c*******************************************************************************
c
c
      subroutine LCURVPL ( nspl,         rnspldp, xnspldp, 
c
     1                     ntau, taudp,  rnormdp, xnormdp, method,
c
     2                     nrad, xraddp, yraddp )
c
c-------------------------------------------------------------------------------
c
      parameter ( npl = 4000 )
c
      character*(*) method
c
      character*10 tauch
c
      double precision  xnspldp(nspl), rnspldp(nspl)
c
      dimension         xnspl(npl),    rnspl(npl)
c
      double precision  taudp(ntau), rnormdp(ntau), xnormdp(ntau)
c
      dimension         tau(npl),    rnorm(npl),    xnorm(npl)
c
      double precision  xraddp(nrad), yraddp(nrad)
c
      dimension         xrad(npl),    yrad(npl)
c
      dimension isymb(14)
c
      data  isymb / 1, 2, 3, 4, 5, 6, 10, 11, 22, 21, 20, 19, 18, 17 /
c
c-------------------------------------------------------------------------------
c
      do k = 1, nspl
         xnspl(k) = xnspldp(k)
         rnspl(k) = rnspldp(k)
      end do
c
      do k = 1, ntau
         tau(k)   = taudp(k)
         rnorm(k) = rnormdp(k)
         xnorm(k) = xnormdp(k)
      end do
c
      do k = 1, nrad
         xrad(k) = xraddp(k)
         yrad(k) = yraddp(k)
      end do
c
  300 continue
c
         call AXCODE ( 0, 0, 2 )
c
         call SIDTEX ( 'L-curve for ' // method, 2, 
c
     1                 '||||b - Ax(|t)||||',   2, 
c
     2                 '||||x(|t)||||',        2, ' ', 2 )
c
         call LINLOG ( 1, 1, 0 )
c
         call HOWPLT ( 0, 1, 2 )
c
         call CURV ( nspl, rnspl, xnspl )
c
         do k = 1, ntau
c
            call HOWPLT ( isymb(k), 0, 2 )
c
            write ( tauch, '(1pe10.3)' ) tau(k)
c
            call LEGND ( '|t =' // tauch )
c
            call CURV ( 1, rnorm(k), xnorm(k) )
c
         end do
c
         call HOWPLT ( 0, 1, 2 )
c
         call CURV ( nrad, xrad, yrad )
c
         call VG
c
      if ( loopin() .eq. 1 ) go to 300
c
      call NEWPAG
c
c
      return
c
      end
