c
c
c------------------------------------------------------------------------
c
      subroutine pplot3 (a,ndim,mstr,mpts,frac,nplt,ltit,xlab,isym,ntou)
c
c------------------------------------------------------------------------
c
c     new plot for ferd, giving window and synthetic windows at once .ge
c     a is matrix to be plotted, dimension of first index  is ndim
c     mpts is maximum number of lines to plot., nplt the number of col
c     lines at beginning and end are omitted if no values (abs).gt.frac
c     isym is vector of plotting symbols, higher indices override.
c     ntou is logical number of output tape.
c     data are printed after scaling. only four columns can be printed.
c     xlab is the abscissa labeling array
c     ltit is the title for xlab (1a6)
c     91 plotting locations, zero can be anywhere.     2-1-69   rwp.
c
c     modified to plot only 3 curves.   9-22-80  dti.
c     mstr  is first line to be plotted.   9-22-80  dti.
c
c------------------------------------------------------------------------
c
      dimension a(ndim,nplt),jj(10)   ,xlab(ndim),gridp(4),buf(10)
c
      character*2  ibuf(92,2),igrid(2,2), isym(3), isymb(4) 
c
c------------       equivalence (buf(1),ibuf(1))
c
        data isymb/ ' ','-','|','+' /
c
c     col 1 of ibuf and igrid for lines not mod 10,row 1 of igrid for 
c     buffer location not mod 10.
c
c------------------------------------------------------------------------
c
      igrid(1,1) = isymb(1)
      igrid(1,2) =isymb(2)
      igrid(2,1) = isymb(3)
      igrid(2,2) =isymb(4)
c *** calculate minimum and maximum value to be plotted
      vmax=-1.0e+35 
      vmin=1.0e+35
      do 30 i=mstr,mpts
         do 30 j=1,nplt
            ael=a(i,j)
         if(ael-vmax) 25,25,24
   24    vmax=ael
   25    if(vmin-ael) 30,30,26
   26    vmin=ael
   30    continue
c *** calculate scale factor
      range=vmax-vmin
      if(range .eq.0.)  range=1.0e-20
      ex =0.43429448*alog(range/0.8)
      if(ex ) 28,29,29
   28    nexp=ex -1.
         go to 35
   29    nexp=ex
   35 fexp =10.0**nexp
      scale = fexp
c     range/scale between 0.8 and 8.
   33 if(range/fexp-4.) 31,32,32
   31 fexp=fexp/2.
         go to 33
c     range/fexp between 4 and 8.graph 9 units wide. now set graph limit
   32 svmax=vmax/fexp
      svmin=vmin/fexp
      delta=fexp/10.
c     now compute ratio,the unit coordinate at left edge of graph.
      nratio=svmin
      if(svmin)34,52,52
   34 if(float(nratio)-svmin-.049    ) 52,52,50
   50 ratio=nratio-1
      go to 55
   52 ratio=nratio
   55 gridlo=ratio*fexp
c     now scan to determine line1 and linel. points will be plotted betw
      flim=frac*range
      do 40 line1=mstr,mpts
         do 40 j=1,nplt
            if(abs( a(line1,j) ) -flim) 40,40,39
   40       continue
   39 line1 = max0(mstr,line1-2)
      do 43 i=mstr,mpts
         linel=mpts + mstr - i
         do 43 j=1,nplt
            if(abs(a(linel,j) ) -flim) 43,43,42
   43        continue
   42 if(linel+2-mpts)48,48,44
   48 linel=linel+2 
c     now set up grid buffers.  will correct changes after each line. 
   44 do 45 l=1,92
         md5=1
         if(mod(l-1,10) .ne. 0) go to 47
   46    md5=2
   47   ibuf(l,1) =igrid(md5,1)
        ibuf(l,2)= igrid(md5,2)
   45    continue
      write (ntou,900) scale,flim,ltit,isym
  900 format(//37x,'scale factor =',1pe9.1,'  minimum value plotted=',
     1       e9.1 / ' ' , 1x, a4, 3x, (5x,a1,2x))
c
c *** graph is now printed, line by line
c
      mb=0
      do 200 ln=line1,linel
      jmax=max0(3,nplt)
c     if all 3 lines are below flim, nothing until values again .gt
      amax=0.
      do 180 j=1,jmax
      alnj = 0.0
      if(j.le.nplt)  alnj=a(ln,j)
      gridp(j)=alnj/scale
      if(abs(alnj) -amax) 180,180,179
  179 amax=abs(alnj)
  180 continue
      if(amax-flim)182,189,189
  182 mb=mb+1
      if(mb-4) 190,183,200
  183 write(ntou,903)
      go to 200
  189 mb=0
  190 i5=1
      if(mod(ln-1,10) .ne. 0) go to 53
   51    i5=2
c     calculate index of symbol for a(ln,.)
   53    do 60 l=1,nplt
            jamp=(a(ln,l) -gridlo)/delta  +1.5
            jj(l)=jamp
   60       ibuf(jamp,i5) = isym(l)
      write(ntou,901) xlab(ln), (gridp(j),j=1,jmax), ln, 
     1                   (ibuf(j,i5),j=1,92) 
  901 format(' ', 1pe8.2, 0p2f7.3, f6.3, i4, 1x, 92a1 ) 
c
c     now repair the print buffer.
      do 70 l=1,nplt
            j=jj(l) 
            md5=1
            if(mod(j-1,10) .ne. 0) go to 72
   71       md5=2
   72       ibuf(j,i5)=igrid(md5,i5)
   70       continue
  200 continue
c *** scale is now printed
      tensp=fexp/scale
      buf(1)=ratio*tensp
      do 220 i=2,10 
      ii=i-1
      buf(i)=buf(ii) +tensp
  220 continue
      write (ntou,902) buf
c
c    --------
  800 return
c    --------
c
  902  format(1h ,31x,9f10.2,f9.2)
  903 format(1h0)
      end 
c
c
c*********************************************************************
c
c
      subroutine aplot(a,ndc,ndr,ncr,nrc,nth,nppp,x,nt)
c
c---------------------------------------------------------------------
c
      dimension a(2*ndc),x(2*ndc),plot(101),escale(11),save(5),sym(5),
c
     1          l(5),loc(5)
c
c  define kkk (at compile time) to test for first or later entry
c
      data kkk/137/ 
c
c---------------------------------------------------------------------
c
      bakf(y,f) = 0.4342944819 * alog(y+f)
      forf(x,f) = 10.0**x-f
c
c     this version plots 3 plus and 2 minus cycles on a p-log scale
c
c                   ndc = 1
c     row plots     ndr = first dimension of a
c                   ncr = number of rows in plot
c                   nrc = number of columns in plot
c
c     column plots     ndc = first dimension of a 
c                      ndr = 1
c                      ncr = number of columns in plot
c                      nrc = number of rows in plot
c
c     every nth row or column of matrix a is plotted
c
c     nppp specifies the number of plots per page (5 max)
c
c     x is the absicsa labeling array
c
c     nt is the number of the output tape
c     ***** (literaly .6. in this version) *****
c
      if(kkk-14232)10,30,10
   10 kkk=14232
      call equiv(sym(1),1h*)
      call equiv(sym(2),1h$)
      call equiv(sym(3),1h=)
      call equiv(sym(4),1h+)
      call equiv(sym(5),1h')
      do 15 i=1,101 
   15 call equiv(plot(i),1h ) 
      do 20 i=1,101,10
   20 call equiv(plot(i),'|') 
   30 index=1
      l1=1
      nplot=1+(ncr-1)/nth
   40 if(nplot-nppp)42,44,44
   42 jp=nplot
      go to 46
   44 jp=nppp
   46 nplot=nplot-nppp
      jp1=jp+1
      amax=1.0e-30
      do 60 j=1,jp
      l(j)=l1+(j-1)*nth
      do 60 i=1,nrc 
      ij=index+ndr*(i-1)+ndc*(j-1)*nth
      if(amax1(a(ij),-10.0*a(ij))-amax)60,60,50                         08/07/67
   50 amax=amax1(a(ij),-10.0*a(ij))                                     08/07/67
   60 continue
      l1=l(nppp)+nth
      nscale = bakf(amax,0)+40.95                                       08/07/67
      nscale=nscale-40
      fnsm3 = nscale-3                                                  08/09/67
      f0 = forf(fnsm3,0)                                                08/09/67
      fix=-20*nscale
      fix=fix+101.5 
      fixn=83.0-fix 
      scale = nscale
      do   70   i=1,11
      dif = i-5
      escale(i) = forf(scale-3.+.5*abs(dif ),f0)                        08/07/67
   70 escale(i) = sign(escale(i),dif)                                   08/09/67
      if(ndr-1)90,90,100
   90 write (6,905) (sym(j),l(j),j=1,jp)                                08/09/67
      write (6,902) escale                                              08/09/67
      go to 110
  100 write (6,906) (sym(j),l(j),j=1,jp)                                08/09/67
      write (6,903) escale                                              08/09/67
  110 do 130 i=1,nrc
      do 121 j=1,jp 
      ij=index+ndr*(i-1)+ndc*(j-1)*nth
      if(a(ij))112,114,114
  112 ij = -20.0*bakf(f0,-a(ij))+fixn                                   08/07/67
      go to 116
  114 ij = 20.0*bakf(a(ij),f0)+fix                                      08/07/67
  116 if(ij)117,117,118
  117 ij=1
  118 if(ij-101)120,120,119
  119 ij=101
  120 loc(j)=ij
      save(j)=plot(ij)
  121 plot(ij)=sym(j)
      write (6,904) i,x(i),plot                                         08/09/6 
      do 130 j=1,jp 
      jj=jp1-j
      ij=loc(jj)
  130 plot(ij)=save(jj)
      index=index+nppp*nth*ndc
      if(nplot)140,140,135
  135 write (6,907)                                                     08/09/67
      go to 40
  140 return
  902 format(1h0,2x,3hrow,8x,5e10.1,             3e10.1,3e9.1)
  903 format(1h0,2x,3hcol,8x,5e10.1,             3e10.1,3e9.1)
  904 format(1x,i3,2x,e10.3,3x,101a1)
  905 format(1h0,13x,5(12x,a1,4h col,i3))
  906 format(1h0,13x,5(12x,a1,4h row,i3))
  907 format(1h1)
      end 
c
c
c************************************************************************
c
c
      subroutine rsplot (a,mr,mc,elab,nr,nc,npl,ipl,ii,apl) 
c
c------------------------------------------------------------------------
c
      common /pltcon/ mi1,mi100,int100,int5,fl45,zero,
c
     1               izero,sevfiv,ten25,sixfiv,two75,fl90
c
      dimension a(mr,mc),elab(mc),ipl(7)
c
      dimension apl(4*mr)
c
ccccccc-------      common /qqextr/ xor,xstp,yor,ystp,xax,yax,idum(3) 
c
c------------------------------------------------------------------------
c
      xlong=8.0
      ylong=6.0
      call reset ('all')
      call nobrdr
      call height ( 0.2 )
      call xaxang ( 45.0 )
      call yaxang ( 45.0 )
c
       emin = r1mach(2)
       emax = -r1mach(2)
       amin = r1mach(2)
       amax = -r1mach(2)
c
       do  5 j=1,nc 
        if (elab(j).le.emin) emin = elab(j)
        if (elab(j).ge.emax) emax = elab(j)
    5  continue
c
      do 20 i=1,npl 
      ip = ipl(i)
      do 10 j=1,nc
        if ( a(ip,j).le.amin ) amin = a(ip,j)
        if ( a(ip,j).ge.amax ) amax = a(ip,j)
   10 continue
   20 continue
c
      call page (12.5,9.5)
c
      if (ii.eq.1) go to 35
c
      call area2d (xlong,ylong)
      call headin ('selected response functions $',100,2.0,1)
      call xname ('energy $',100)
      call yname ('response $',100)
      go to 38
c
   35 continue
c
      call area2d (xlong,ylong)
      call headin ('selected windows $',100,2.0,1)
      call xname ('energy $',100)
      call yname ('window response $',100)
c
   38 continue
c
      call xticks(int5)
      call yticks(int5)
      call xaxang(fl45)
      call yaxang(zero)
c
        xor = emin
        xmax = emax 
        yor = amin
        ymax = amax 
c
      call graf (xor,'scale',xmax,yor,'scale',ymax)
      do 50 ip=1,npl
      i=ipl(ip)
      do 40 j=1,nc
      apl(j)=a(i,j) 
   40 continue
      call curve (elab,apl,nc,izero)
   50 continue
      call frame
      call endpl(izero)
c
      return
      end 
c
c
c************************************************************************
c
c
      subroutine plotb ( b,s,mr,nr,babsc,xwrk,ywrk,mw,
c
     1                   workv1,workv2,workv3,workv4 )
c
c-----------------------------------------------------------------------
c
      common /pltcon/ mi1,mi100,int100,int5,fl45,zero,
c
     1                izero,sevfiv,ten25,sixfiv,two75,fl90
c
      common /plong/ xpage,ypage,xlong,ylong
c
      dimension b(mr),s(mr),babsc(mr),xwrk(mw),ywrk(mw)
      dimension workv1(2*mr),workv2(2*mr),workv3(2*mr),workv4(2*mr)
c
c-----------------------------------------------------------------------
c
      xlong = 8.0
      ylong = 6.0
      xpage = xlong + 3.0
      ypage = ylong + 3.0
      nr2 = 2*nr
      bmin = r1mach(2)
      do 10 i=1,nr
      if(b(i).le.0) go to 5
      if (b(i).lt.bmin) bmin=b(i)
    5 continue
      xwrk(i) = babsc(i)
      ywrk(i) = b(i) + s(i)
      nrpi = nr + i 
      xwrk(nrpi) = babsc(i)
      ywrk(nrpi) = b(i) - s(i)
   10 continue
c
c  linear plot of the measured counts
c
      nrp1 = nr + 1 
c
      call pltset ( babsc,nr,ywrk(1),ywrk(nrp1),1,1,nr,workv1,workv2, 
     1              workv3,npl,2,workv4,workv4 )
c
      call cplot ( workv1,workv3,workv2,npl,2,workv4,
     1             'observed counts $','pulse height $',
     2             'normalized counts +- one std. dev. $',
     3              .false. ) 
c
c
cc  logarithmic plot of measured counts 
c
      do 20 i=1,nr
      nrpi = nr + i 
      if((ywrk(i).gt.0.0).and.(ywrk(nrpi).gt.0.0)) go to 20 
      ywrk(i) = bmin
      ywrk(nrpi) = bmin
   20 continue
c
      call gphbgn (xwrk,ywrk,nr2,xorig,xcycle,yorig,ycycle) 
c
      call page (12.5,9.5)
      call area2d (xlong,ylong)
      call headin ('observed counts $',100,2.0,1) 
      call xname ('pulse height $',100) 
      call yname ('normalized counts +- 1 std. dev. $',100) 
c
      call gphend (xorig,yorig,xcycle,ycycle,xwrk,ywrk,nr)
c
      return
      end 
c
c
c**********************************************************************
c
c
      subroutine phplot ( plo, pup, n1, nw, nwrpl, kr,pabsc,
c
     1                    xwrk,ywrk,mwrk,
c
     1                    pest,workv1,workv2,workv3,workv4,label )
c
c---------------------------------------------------------------------
c
      character*80 label
c
      common /pltcon/ mi1,mi100,int100,int5,fl45,zero,
c
     1               izero,sevfiv,ten25,sixfiv,two75,fl90
c
      common /plong/ xpage,ypage,xlong,ylong
c
c
      dimension plo(kr),pup(kr),pabsc(kr),xwrk(mwrk),ywrk(mwrk)
c
      dimension pest(kr),workv1(2*kr),workv2(2*kr),workv3(2*kr),
c
     1                   workv4(2*kr)
c
c----------------------------------------------------------------------
c
      xlong = 8.0
      ylong = 6.0
      xpage = xlong + 3.0
      ypage = ylong + 3.0
c
      if(n1.gt.0) go to 4
c
c  draw linear histogram plots only
c
      npl = 2*nw
      npl2 = 2*npl
      delo2 = (pabsc(2)-pabsc(1))/2.0
      xwrk(1) = pabsc(1) - delo2
      ywrk(1) = plo(1)
      xwrk(npl+1) = xwrk(1)
      ywrk(npl+1) = pup(1)
      nwm1 = nw - 1 
c
      do 3 i=1,nwm1 
      ip = 2*i
      ip1 = ip + 1
      xwrk(ip) = (pabsc(i)+pabsc(i+1))/2.0
      ywrk(ip) = plo(i)
      xwrk(ip1) = xwrk(ip)
      ywrk(ip1) = plo(i+1)
      ipu = ip + npl
      ipu1 = ipu + 1
      xwrk(ipu) = xwrk(ip)
      ywrk(ipu) = pup(i)
      xwrk(ipu1) = xwrk(ip)
      ywrk(ipu1) = pup(i+1)
    3 continue
      delo2 = (pabsc(nw)-pabsc(nw-1))/2.0
      xwrk(npl) = pabsc(nw) + delo2
      ywrk(npl) = plo(nw)
      xwrk(2*npl) = xwrk(npl) 
      ywrk(2*npl) = pup(nw)
c
      go to 65
c
c
    4 continue
c
      npl = nwrpl
c
c
      call pltset ( pabsc(n1),npl,plo(n1),pup(n1),1,1,npl,workv1,
c
     1              workv2,workv3,nplot,3,pest(n1),workv4 ) 
c
c
      call cplot ( workv1,workv3,workv2,nplot,3,workv4,label,
c
     1             'energy $','flux $',.false. )
c
c
      npl2 = 2*npl
      isc = 0
c
      pumin = r1mach(2)
      do 2 i=1,npl
       in1 = i + n1 - 1
       if ( (pup(in1).gt.0.0).and.(pup(in1).lt.pumin) ) then
                                                        pumin = pup(in1)
       end if
    2 continue
c
      do 10 i=1,npl 
      in1 = i + n1 - 1
      if ((plo(in1).gt.0.0).and.(plo(in1).gt.pumin*0.01)) go to 5
      if(pup(in1).le.0.0) go to 10
      isc = isc + 1 
      xwrk(isc) = pabsc(in1)
      ywrk(isc) = pup(in1)
       go to 7
c
    5 continue
c
      isc = isc + 1 
      xwrk(isc) = pabsc(in1)
      ywrk(isc) = plo(in1)
c
    7 continue
c
      ipnpl = npl + isc
      xwrk(ipnpl) = pabsc(in1)
      ywrk(ipnpl) = pup(in1)
  10    continue
      nsc = isc
      nsc2 = 2*nsc
      if(nsc.ge.npl) go to 30 
      do 20 i=1,nsc 
      isc = nsc + i 
      ipup = npl + i
      xwrk(isc) = xwrk(ipup)
      ywrk(isc) = ywrk(ipup)
   20 continue
c
   30 continue
c
      call gphbgn (xwrk,ywrk,nsc2,xorig,xcycle,yorig,ycycle)
c
      call page (12.5,9.5)
      call area2d (xlong,ylong)
      call headin ('estimated spectrum $',100,2.0,1)
      call xname ('energy $',100)
      call yname ('flux $',100)
c
      do 50 i=1,npl 
      in1 = i + n1 - 1
      xwrk(i) = pabsc(in1)
      ywrk(i) = plo(in1)
      if(ywrk(i).le.yorig) ywrk(i) = yorig
      ipnpl = i + npl
      xwrk(ipnpl) = pabsc(in1)
      ywrk(ipnpl) = pup(in1)
      if(ywrk(ipnpl).le.yorig) ywrk(ipnpl) = yorig
  50    continue
c
      call gphend (xorig,yorig,xcycle,ycycle,xwrk,ywrk,npl) 
c
      do 60 i=1,npl 
      in1 = i + n1 - 1
      xwrk(i) = pabsc(in1)
      ywrk(i) = plo(in1)
      ipnpl = npl + i
      xwrk(ipnpl) = pabsc(in1)
      ywrk(ipnpl) = pup(in1)
   60 continue
c
c
   65 continue
c
c
c  do two copies of final plot
c
      do 80 iipp=1,2
c
      call reset ('all')
      call hldplt
      call linplt (xwrk,ywrk,npl2)
c
c  the following statements killed (3/20/86) for conversion 
c   to laser graphics plots
c
c-----call page (sevfiv,ten25)
c-----call physor (sixfiv,two75)
c-----call bangle (fl90)
c
      call page (12.5,9.5)
      call area2d (xlong,ylong)
      call headin ('estimated spectrum $',100,2.0,1)
      call xname ('energy $',100)
      call yname ('flux $',100)
      call xticks(int5)
      call yticks(int5)
      call xaxang(fl45)
      call yaxang(zero)
c
        xor = r1mach(2)
        xmax = -r1mach(2)
        yor = r1mach(2)
        ymax = -r1mach(2)
       do 67 i=1,npl2
        if (xwrk(i).le.xor) xor = xwrk(i)
        if (xwrk(i).ge.xmax) xmax = xwrk(i)
        if (ywrk(i).le.yor) yor = ywrk(i)
        if (ywrk(i).ge.ymax) ymax = ywrk(i)
   67  continue
c
      call graf (xor,'scale',xmax,yor,'scale',ymax)
      call curve (xwrk,ywrk,npl,izero)
c
      nplp1 = npl + 1
c
      call curve (xwrk,ywrk(nplp1),npl,izero)
c
      call frame
      call endpl (izero)
c
      do 76 i=1,npl2
        if (ywrk(i).lt.0.0) ywrk(i) = 0.0
   76 continue
c
   80 continue
c
c
      return
      end 
c
c
c***********************************************************************
c
c
      subroutine gphbgn (x,y,n,xorig,xcy,yorig,ycy)
c
c-----------------------------------------------------------------------
c
c        this subroutine stolen from b. l. kirk
c
      common /pltcon/ mi1,mi100,int100,int5,fl45,zero,
c
     1                izero,sevfiv,ten25,sixfiv,two75,fl90
c
      common /plong/ xpage,ypage,xlong,ylong
c
      dimension  x(n), y(n)
c
c----------------------------------------------------------------------
c
      call reset ('all')
      call nobrdr
      call height ( 0.2 )
      call xaxang ( 45.0 )
      call yaxang ( 45.0 )
c
      xmin = r1mach(2)
      ymin = r1mach(2)
      xmax = 0.0
      ymax = 0.0
      do 5 i=1,n
      if(x(i).lt.xmin) xmin = x(i)
      if(x(i).gt.xmax) xmax = x(i)
      if(y(i).lt.ymin) ymin = y(i)
      if(y(i).gt.ymax) ymax = y(i)
    5 continue
c
      call algplt ( xmin, xmax, xlong, xorig, xcycle )
c
      call algplt ( ymin, ymax, ylong, yorig, ycycle )
c
      xcy = xcycle
      ycy = ycycle
c
c
c
      return
      end 
c
c
c**********************************************************************
c
c
      subroutine gphend (xorig,yorig,xcycle,ycycle,x,y,n)
c
c----------------------------------------------------------------------
c
c        this subroutine stolen from b. l. kirk
c
      common /pltcon/ mi1,mi100,int100,int5,fl45,zero,
c
     1                izero,sevfiv,ten25,sixfiv,two75,fl90
c
c
      dimension  x(2*n), y(2*n)
c
c---------------------------------------------------------------------
c
      call xticks(int5)
      call yticks(int5)
      call xaxang(fl45)
      call yaxang(zero)
c
c
      call loglog (xorig,xcycle,yorig,ycycle)
      call curve (x,y,n,izero)
      do 10 i=1,n
      npi = n+i
      y(i) = y(npi) 
   10 continue
c
c
      call curve (x,y,n,izero)
      call frame
      call endpl(izero)
      return
      end 
c
c
c************************************************************************
c
      subroutine ppunch(plab,plo,pup,nw,ident,nout)
c
c-----------------------------------------------------------------------
c
c   punches plo-pup solution in user format
c
      dimension plab(nw),plo(nw),pup(nw),ident(20)
c
c------------------------------------------------------------------------
c
      write(nout,701) ident
      write(nout,702) nw
      write(nout,703) (plab(i),plo(i),pup(i),i=1,nw)
c
      return
c
  701 format(20a4)
  702 format(i10)
  703 format(6e12.5)
c
c
      end 
c
c
c***********************************************************************
c
c
      subroutine arplot (a,mr,mc,babsc,nr,nc,npl,ipl,ii,apl)
c
c-----------------------------------------------------------------------
c
      common /pltcon/ mi1,mi100,int100,int5,fl45,zero,
c
     1                izero,sevfiv,ten25,sixfiv,two75,fl90
c
      dimension  a(mr,mc),babsc(mr),ipl(7)
c
      dimension  apl(4*mr)
c
c------------------------------------------------------------------------
c
      xlong=8.0
      ylong=6.0
      xpage=xlong+3.0
      ypage=ylong+3.0
      call reset ('all')
c
      amin = 0.0
      amax=0.0
      do 20 j=1,npl 
      ip = ipl(j)
      do 10 i=1,nr
      if(a(i,ip).le.amax) go to 10
      amax = a(i,ip)
      imax=ip
   10 continue
   20 continue
      xmin = r1mach(2)
      xmax = -r1mach(2)
      do 30 i=1,nr
      apl(i)=a(i,imax)
      if (babsc(i).lt.xmin) xmin = babsc(i)
      if (babsc(i).gt.xmax) xmax = babsc(i)
   30 continue
c
      call nobrdr
      call height ( 0.2 )
      call xaxang ( 45.0 )
      call yaxang ( 45.0 )
      call page (12.5,9.5)
c
      if (ii.eq.1) go to 35
      call area2d (xlong,ylong)
      call headin ('selected response functions $',100,2.0,1)
      call xname ('pulse height $',100) 
      call yname ('response $',100)
      go to 38
c
   35 continue
c
      call area2d (xlong,ylong)
      call headin ('selected windows $',100,2.0,1)
      call xname ('energy $',100)
      call yname ('window response $',100)
c
   38 continue
c
      call xticks(int5)
      call yticks(int5)
c
      call algplt ( xmin, xmax, xlong, xorig, xcycle )
c
      call axsplt ( amin, amax, ylong, yorig, ystep, yaxis )
c
      call xlog ( xorig, xcycle, yorig, ystep )
c
      do 50 ip=1,npl
      j=ipl(ip)
      do 40 i=1,nr
      apl(i)=a(i,j) 
   40 continue
      call curve (babsc,apl,nr,izero)
   50 continue
      call frame
      call endpl(izero)
c
      return
      end 
c
c
c****************************************************************************
c
c
      subroutine pltset (smesh,nc,filo,fihi,iwlo,iwstp,iwup,
c
     1                   workv1,workv2,workv3,npl,ncv,fi3,workv4 )
c
c---------------------------------------------------------------------------
c
c
      dimension smesh(nc)
c
      dimension  filo(iwup), fihi(iwup), workv1(2*iwup),
c
     1           workv2(2*iwup),
c
     2           workv3(2*iwup), fi3(iwup), workv4(2*iwup)
c
c
c--------------------------------------------------------------------------
c
      if (iwup.lt.nc) go to 50
c
c
c  pointwise plots  (no averaging)
c
      npl = nc
      do 10 i=1,npl 
       workv1(i) = smesh(i)
       workv2(i) = fihi(i)
       workv3(i) = filo(i)
        if (ncv.gt.2) workv4(i) = fi3(i)
   10 continue
c
      go to 60
c
c
c
   50 continue
c
c
c  histogram plots (averaging)
c
      ipl = 0
      npl = 0
      do 75 iw=iwlo,iwup,iwstp
       ipl = ipl + 1
       npl = npl + 2
       workv1(npl-1) = smesh(iwstp*(ipl-1)+1)
       workv2(npl-1) = fihi(ipl)
       workv3(npl-1) = filo(ipl)
        if (ncv.gt.2) workv4(npl-1) = fi3(ipl)
       workv1(npl) = smesh (iwstp*ipl+1)
       workv2(npl) = fihi(ipl)
       workv3(npl) = filo(ipl)
        if (ncv.gt.2) workv4(npl) = fi3(ipl)
   75 continue
c
c
 60   continue
c
c
      return
c
c
      end 
c
c
c*************************************************************************
c
c
      subroutine cplot (x,ylo,yup,np,ncv,y3,label,xlab,ylab,iterm)
c
c------------------------------------------------------------------------
c
      character*40 xlab,ylab
      character*80 label
c
      logical iterm 
c
      dimension x(np),ylo(np),yup(np),y3(np)
c
c-----------------------------------------------------------------------
c
      if (ncv.eq.2) go to 25
c
c
c  first plot the point estimate only
c
c
      xlong = 7.0
      ylong = 5.0
      if(iterm) xlong = 8.0
      if(iterm) ylong = 6.0
      call reset ('all')
      call nobrdr
      call height(0.2)
      call xaxang (45.)
      call yaxang (45.)
c
c
c  for tektronix plots
c
      call page (12.5,9.5)
c
c
      call area2d (xlong,ylong)
      call headin (label,100,1.0,1)
      call xname (xlab,100)
      call yname (ylab,100)
c
      xmin = x(1)
      ymin = y3(1)
      ymax = y3(1)
      xmax = x(np)
      do 20 i=2,np
       if (x(i).lt.xmin) xmin = x(i)
       if (x(i).gt.xmax) xmax = x(i)
       if (y3(i).lt.ymin) ymin = y3(i)
       if (y3(i).gt.ymax) ymax = y3(i)
   20 continue
c
c
      call graf (xmin,'scale',xmax,ymin,'scale',ymax)
      call frame
      call curve (x,y3,np,0)
c
      call endpl(0) 
c
c  now plot the lower and upper bounds only
c
c
 25   continue
c
      xlong = 7.0
      ylong = 5.0
      if(iterm) xlong = 8.0
      if(iterm) ylong = 6.0
      call reset ('all')
      call nobrdr
      call height(0.2)
      call xaxang (45.)
      call yaxang (45.)
c
c
c  for tektronix plots
c
      call page (12.5,9.5)
c
c
      call area2d (xlong,ylong)
      call headin (label,100,1.0,1)
      call xname (xlab,100)
      call yname (ylab,100)
c
      xmin = x(1)
      ymin = ylo(1) 
      ymax = yup(1) 
      xmax = x(np)
      do 30 i=2,np
       if (x(i).lt.xmin) xmin = x(i)
       if (x(i).gt.xmax) xmax = x(i)
       if (ylo(i).lt.ymin) ymin = ylo(i)
       if (yup(i).gt.ymax) ymax = yup(i)
   30 continue
c
c
      call graf (xmin,'scale',xmax,ymin,'scale',ymax)
      call frame
      call curve (x,ylo,np,0) 
      call curve (x,yup,np,0) 
c
      call endpl(0) 
c
c
      if (ncv.eq.2) return
c
c
c
c  now plot the bounds and the point estimate
c
c
      xlong = 7.0
      ylong = 5.0
      if(iterm) xlong = 8.0
      if(iterm) ylong = 6.0
      call reset ('all')
      call nobrdr
      call height(0.2)
      call xaxang (45.)
      call yaxang (45.)
c
c
c  for tektronix plots
c
      call page (12.5,9.5)
c
c
      call area2d (xlong,ylong)
      call headin (label,100,1.0,1)
      call xname (xlab,100)
      call yname (ylab,100)
c
      xmin = x(1)
      ymin = ylo(1) 
      ymax = yup(1) 
      xmax = x(np)
      do 40 i=2,np
       if (x(i).lt.xmin) xmin = x(i)
       if (x(i).gt.xmax) xmax = x(i)
       if (ylo(i).lt.ymin) ymin = ylo(i)
       if (yup(i).gt.ymax) ymax = yup(i)
   40 continue
c
c
      call graf (xmin,'scale',xmax,ymin,'scale',ymax)
      call frame
      call curve (x,ylo,np,0) 
      call curve (x,yup,np,0) 
      if (ncv.gt.2) then
                         call dot
                         call curve (x,y3,np,0)
                         call reset ('dot')
      end if
c
      call endpl(0) 
c
c
c
      return
c
      end 
c
c
c------------------------------------------------------
c
c
      subroutine vecpr (vec,n,vecnam,iw10)
c
      logical ipr
      character*4 vecnam
      common /blokip/ ipr
      dimension vec(1)
c
c  if (iw10.lt.0) don't print,  write on unit 10 0nly
c
      if (iw10.ge.0) print 10,vecnam
   10 format (1h //1h ,a4/)
      if (ipr) write (10,10) vecnam
      do 50 jl=1,n,5
       jlo = jl
       jup = jl + 4 
       if (jup.gt.n) jup = n
       if (iw10.ge.0) print 20,jlo,(vec(j),j=jlo,jup)
   20  format (1h ,i4,1p5e14.6)
      if (ipr) write (10,20) jlo,(vec(j),j=jlo,jup)
   50 continue
c
      return
      end 
