Copyright (c) 2002 Peter Guntert. All rights reserved.
c     ==================================================================
c
c               GRAF: Make graphs and histograms
c               ________________________________
c
c
c
c               The input data file(s) consist of data and comment
c               lines. Comment lines are not required and are either
c               blank or start with '#', all other lines are data
c               lines. Optionally, a first data line contains the
c               number m of curves. If this line is missing, the
c               program sets this number by itself. The following n
c               data lines contain the m+1 real numbers x(i),y_1(i),
c               ...,y_m(i) (i=1,...,n). The format of the data lines
c               is free. An input file may have only comment lines.
c               Special comments can be used to change the default
c               values for the various display options and parameters:
c
c                 #RECTANGLE        <new default value>
c                 #VIEWPORT         <new default value> | '*'
c                 #MINMAX           <new default value> | '*'
c                 #LETTERSIZE       <new default value>
c                 #XTEXT            <new default value>
c                 #YTEXT            <new default value>
c                 #REPRESENTATIONS  <new default value>
c                 #CONNECTIONS      <new default value>
c                 #LINESTYLES       <new default value>
c
c               Here the keywords are given in the order of the inter-
c               active dialog, <new default value> is a character string
c               that is a legal answer to the corresponding question,
c               '*' indicates the device or input dependent default
c               values. The keywords must follow the comment sign '#'
c               immediately. They may be (unambiguously) abbreviated.
c
c               The special comments #NEXTCURVE and #NEXTPLOT indicate
c               that the data (and possibly special comment lines)
c               following them refer to a new set of curves (usually
c               with different x(i) values) in the same plot or to a
c               new plot on the same sheet of paper, recpectively.
c               The special comment #PLOTAGAIN works like #NEXTPLOT
c               except that the same data are used as for the previous
c               plot.
c
c               The special comments #XLABEL and #YLABEL can be used to
c               label the x and y axes in a non-default way:
c
c                 #XLABEL <x-value> [<label text>]
c                 #YLABEL <y-value> [<label text>]
c
c               If one or several #XLABEL comments are present in the
c               data file the whole default labelling of the x axis is
c               suppressed, and labels will only be placed at the
c               locations defined with <x-value>. The same applies for
c               the y axis. The special comments
c
c                 #XTICSONLY
c                 #YTICSONLY
c
c               suppress the labels but not the tics.
c
c               The special comments #LEFTTEXT, #CENTERTEXT and
c               #RIGHTTEXT can be used to place text of any size at
c               any position in the plot. The text will be left-
c               adjusted, centered or rightadjusted, respectively,
c               and the syntax is
c
c                 #LEFTTEXT   <x-value> <y-value> <rel. size> <text>
c                 #CENTERTEXT <x-value> <y-value> <rel. size> <text>
c                 #RIGHTTEXT  <x-value> <y-value> <rel. size> <text>
c
c               <rel. size> denotes the relative letter size with
c               respect to the letter size used for the labelling
c               of the axes.
c               The special comments #LEFTSTRING, #CENTERSTRING and
c               #RIGHTSTRING work in the same way except that <x-value>
c               and <y-value> are given in fractional coordinates
c               between 0 and 1.
c               <text> that starts with a backslash is printed using
c               the Symbol font.
c
c               Relative viewport specifications: Instead of defining
c               the viewport with coordinates, it may also be given
c               in one of the following symbolic formats:
c
c                  >X, >>X, >>>X, etc.   same size, right of plot X
c                  <X, <<X, <<<X, etc.   same size, left of plot X
c                  vX, vvX, vvvX, etc.   same size, below plot X
c                  ^X, ^^X, ^^^X, etc.   same size, above plot X
c
c               X denotes a previous plot (A the first, B the second
c               plot etc.). By default, the previous plot is used.
c               The number of direction signs (>, <, v, or ^)
c               determines the spacing between the two plots. One
c               direction sign places the two plots directly side by
c               side, additional direction signs add intermediate
c               space.
c
c               Written in standard FORTRAN-77.
c               Peter G"untert, 1-3-1989
c
c               DI-3000 routines replaced by Postscript
c                 Tai-he Xia 27-12-93
c               Revised version, Peter G"untert, 18-08-1994
c
c     ------------------------------------------ example input data file
c     # example file
c
c     # This example file contains the data for two plots with four and
c     # two curves, respectively.
c
c     #VIEWPORT -0.8 0.8 0.0 0.8
c     #MINMAX    0.0 8.0 -10.0 100.0
c     #REPRESENTATIONS .--
c     #CONNECTIONS -SL
c     3
c     1.00    3.456   5.786   3.900
c     1.40    4.394   4.786  53.000
c     3.20    5.176   7.786  73.700
c     5.95    6.238  -9.786  -3.400
c     6.00    7.490   0.786  40.500
c     #NEXTCURVE
c     #REPRESENTATIONS -
c     #CONNECTIONS S
c     1
c     1.50     12.89
c     1.79      9.70
c     3.50      4.45
c     #NEXTPLOT
c     #VIEWPORT -0.8 0.8 -0.9 -0.1
c     #REPRESENTATIONS -
c     #CONNECTIONS S
c     2
c     2.3     12.89   66.1
c     3.8      9.70   12.8
c     6.0      4.45   17.3
c     9.1     17.47   19.5
c     ------------------------------------------------------------------
c
      program graf
c
      parameter (maxn=43000,maxcrv=20,maxplt=20,maxlin=132,maxla=500,
     *           pi=3.141592653589793,twopi=2.0*pi)
      dimension x(maxn),y(maxn,maxcrv),v(4),ytmp(maxn),indx(maxn),
     *          xsplin(maxn),ysplin(maxn),xla(maxla),yla(maxla),
     *          vxmini(maxplt),vxmaxi(maxplt),
     *          vymini(maxplt),vymaxi(maxplt)
      logical yes,cross,warned,nocnct,firfil,firplt,fircur,firdat,
     *        previn,inside,xticso,yticso,rotate,illexp,readit,
     *        userep,usecon
      character*(maxcrv) repres,conect,linstl
      character*20 str
      character*100 xlatxt(maxla),ylatxt(maxla),text(maxla),deffnt
c      character*(maxlin) line,filnam
      character*(maxlin) line
c
      parameter (jrec=1,jvie=2,jmin=3,jcro=4,jlet=14,jlab=5,jdvi=6,
     *           jxtx=7,jytx=8,jrep=9,jcon=10,jlin=11,jnxc=12,jnxp=13,
     *           jxla=15,jyla=16,jrtx=17,jltx=18,jctx=19,
     *           jxto=20,jyto=21,jrst=22,jlst=23,jcst=24,
     *           jpla=25,jsrt=26,jfnt=27,jrot=28,nkey=28)
      dimension ldef(nkey)
      logical havdef,cnfdef,nxtcur,nxtplt,pagain,sort
      character key(nkey)*20,def(nkey)*80
      common /defdat/cnfdef
      common/transf/xmin,xmax,ymin,ymax,pxmin,pxmax,pymin,pymax
      data key(jrec),def(jrec)/'RECTANGLE'        ,'n'/
      data key(jvie),def(jvie)/'VIEWPORT'         ,'*'/
      data key(jmin),def(jmin)/'MINMAX'           ,'*'/
      data key(jcro),def(jcro)/'CROSS'            ,'n'/
      data key(jlet),def(jlet)/'LETTERSIZE'       ,'0.025'/
      data key(jlab),def(jlab)/'LABELS'           ,'y'/
      data key(jdvi),def(jdvi)/'DVIFONTS'         ,'n'/
      data key(jxtx),def(jxtx)/'XTEXT'            ,' '/
      data key(jytx),def(jytx)/'YTEXT'            ,' '/
      data key(jrep),def(jrep)/'REPRESENTATIONS'  ,'.'/
      data key(jcon),def(jcon)/'CONNECTIONS'      ,'-'/
      data key(jlin),def(jlin)/'LINESTYLES'       ,'0'/
      data key(jxla)          /'XLABEL'/
      data key(jyla)          /'YLABEL'/
      data key(jnxc)          /'NEXTCURVE'/
      data key(jnxp)          /'NEXTPLOT'/
      data key(jrtx)          /'RIGHTTEXT'/
      data key(jltx)          /'LEFTTEXT'/
      data key(jctx)          /'CENTERTEXT'/
      data key(jrst)          /'RIGHTSTRING'/
      data key(jlst)          /'LEFTSTRING'/
      data key(jcst)          /'CENTERSTRING'/
      data key(jxto)          /'XTICSONLY'/
      data key(jyto)          /'YTICSONLY'/
      data key(jrot)          /'ROTATE'/
      data key(jpla)          /'PLOTAGAIN'/
      data key(jsrt)          /'SORT'/
      data key(jfnt)          /'FONT'/
c
c      print '(/A/)',
c     *      'GRAF: make graphs and histograms (postscript version)'
      nplot=0
      havdef=.false.
      cnfdef=.true.
      fircur=.true.
      firplt=.true.
      firfil=.true.
      readit=.true.
      do j=1,nkey
        ldef(j)=max(1,lenstr(def(j)))
      end do
      deffnt='Times-Roman'
c     --------------------------------------------------- open data file
  150 continue
c  150 print '(a)','input data file (default extension: grf):'
c      read (*,'(A)') line
c  150 call readfn ('input data file','.grf',' ',line)
c      if (line.eq.' ') go to 930
c      call appext (line,'grf')
c      open (unit=1,file=line,status='OLD',err=180)
c      go to 190
c  180 print '(a)','error: input file cannot be opened'
c      go to 150
c     --------------------------------------------- open postscript file
  190 if (firfil) then
c        l=lenstr(line)
c        if (line(max(1,l-3):l).eq.'.grf') then
c          line(l-3:l)='.ps'
c        else
c          line(l+1:)='.ps'
c        end if
c        call readfn ('output postscript file','.ps',
c     *              line(1:lenstr(line)),filnam)
c        open (11,file=filnam,status='unknown')
        open (2,status='SCRATCH')
        call jbegin
      else
        yes=cnfdef
        cnfdef=.true.
c        call yesno ('new plot','n',fircur)
        fircur=.false.
        cnfdef=yes
        if (fircur) firplt=.false.
      end if
c     --------------------------------------------------- read data file
  220 rewind (2)
      nxtcur=.false.
      nxtplt=.false.
      pagain=.false.
      sort=.false.
      xticso=.false.
      yticso=.false.
      rotate=.false.
      ntxt=0
      if (fircur) then
        nxla=0
        nyla=0
      end if
      firdat=.true.
  250 continue
        read (*,'(A)',end=300) line
        do i=1,maxlin
          if (line(i:i).gt.' ') go to 280
        end do
  280   if (i.gt.maxlin) go to 250
        if (line(i:i).eq.'#') then
          do j=i+1,maxlin
            if (line(j:j).ge.'a' .and. line(j:j).le.'z')
     *        line(j:j)=char(ichar(line(j:j))-ichar('a')+ichar('A'))
            if (line(j:j).le.' ') go to 283
          end do
  283     if (j.gt.i+1) then
            call fulnam (line(i+1:j),key,nkey,inam)
c            print *,line(1:70),line(i+1:j),inam
            if (inam.eq.jnxc) then
              nxtcur=.true.
            else if (inam.eq.jnxp) then
              nxtplt=.true.
            else if (inam.eq.jpla) then
              nxtplt=.true.
              pagain=.true.
            else if (inam.eq.jsrt) then
              sort=.true.
            else if (inam.eq.jfnt) then
              deffnt=line(j+1:maxlin)
              call lefstr (deffnt)
            else if (inam.eq.jxla) then
              nxla=nxla+1
              xlatxt(nxla)=line(j+1:maxlin)
            else if (inam.eq.jyla) then
              nyla=nyla+1
              ylatxt(nyla)=line(j+1:maxlin)
            else if (inam.eq.jrtx) then
              ntxt=ntxt+1
              text(ntxt)='3 '//line(j+1:maxlin)
            else if (inam.eq.jltx) then
              ntxt=ntxt+1
              text(ntxt)='1 '//line(j+1:maxlin)
            else if (inam.eq.jctx) then
              ntxt=ntxt+1
              text(ntxt)='2 '//line(j+1:maxlin)
            else if (inam.eq.jrst) then
              ntxt=ntxt+1
              text(ntxt)='-3 '//line(j+1:maxlin)
            else if (inam.eq.jlst) then
              ntxt=ntxt+1
              text(ntxt)='-1 '//line(j+1:maxlin)
            else if (inam.eq.jcst) then
              ntxt=ntxt+1
              text(ntxt)='-2 '//line(j+1:maxlin)
            else if (inam.eq.jxto) then
              xticso=.true.
            else if (inam.eq.jyto) then
              yticso=.true.
            else if (inam.eq.jrot) then
              rotate=.true.
            else if (inam.gt.0) then
              havdef=.true.
              def(inam)=' '
              do k=j+1,maxlin
                if (line(k:k).gt.' ') go to 285
              end do
  285         if (k.le.maxlin) def(inam)=line(k:maxlin)
              ldef(inam)=lenstr(def(inam))
            else if (inam.eq.-2) then
              call errmsg ('GRAF: Ambiguous keyword '//line(i+1:j))
            end if
          end if
          if (nxtcur .or. nxtplt) go to 300
        else if (readit) then
          l=index(line,'#')
          if (l.gt.0) line(l:)=' '
          if (firdat) then
            firdat=.false.
            call lefstr (line)
            l=lenstr(line)
            do k=1,l
              if (line(k:k).lt.'0' .or. line(k:k).gt.'9') go to 291
            end do
            go to 295
  291       ncrv=-1
            do k=1,l
              if (line(k:k).gt.' ' .and. line(k+1:k+1).le.' ')
     *          ncrv=ncrv+1
            end do
            write (2,'(I4)') ncrv
          end if
  295     write (2,'(A)') line
        end if
      go to 250
c     ------------------------------------------------- read data points
  300 if (readit) then
        rewind (2)
        n=0
        read (2,*,end=308) ncrv
        if (ncrv.gt.maxcrv) call errmsg ('GRAF: Too many curves')
  306   continue
          n=n+1
          read (2,*,end=308) x(n),(y(n,i),i=1,ncrv)
        go to 306
c  308   if (.not.(nxtcur .or. nxtplt)) close (unit=1)
  308   n=n-1
c        if (n.le.0) go to 150
        if (n.le.0) go to 930
        if (n.gt.maxn) call errmsg ('GRAF: Too many data points')
c       ------------------------------------------- use always defaults?
        if (firplt .and. fircur .and. havdef) then
c          call yesno ('use always defaults','n',yes)
          yes=.true.
          cnfdef=.not.yes
        end if
c       ----------------------------------------------- sort data points
        if (sort) then
          call rsort (x,indx,n)
          do j=1,ncrv
            do i=1,n
              ytmp(i)=y(i,j)
            end do
            do i=1,n
              y(i,j)=ytmp(indx(i))
            end do
          end do
        end if
      end if
c     --------------------------------------- display data point summary
c      print '(A,I12)','data file read: number of curves     :',ncrv
c      print '(A,I12)','                number of data points:',n
c     --------------------------------------- set extreme display values
      if (fircur) then
        nplot=nplot+1
        if (nplot.gt.maxplt) call errmsg ('GRAF: Too many plots')
        xmin=x(1)
        xmax=x(n)
        ymin= 1.0e+30
        ymax=-1.0e+30
        call lefstr (def(jrep))
        call lefstr (def(jcon))
        do j=1,ncrv
          if (def(jrep)(j:j).eq.'-' .or.
     *        def(jrep)(j:j).eq.'y' .or. def(jrep)(j:j).eq.'Y') then
            userep=.false.
          else if (def(jrep)(j:j).gt.' ') then
            userep=.true.
          end if
          if (def(jcon)(j:j).eq.'-') then
            usecon=.false.
          else if (def(jcon)(j:j).gt.' ') then
            usecon=.true.
          end if
          if (userep .or. usecon) then
            do i=1,n
              ymin=min(ymin,y(i,j))
              ymax=max(ymax,y(i,j))
            end do
          end if
        end do
c        print '(A,G12.4),'                minimal x value      :',xmin
c        print '(A,G12.4),'                maximal x value      :',xmax
c        print '(A,G12.4),'                minimal y value      :',ymin
c        print '(A,G12.4),'                maximal y value      :',ymax
        if (xmax.eq.xmin) then
          if (xmin.eq.0.0) then
            xmin=-1.0
            xmax=+1.0
          else
            xmin=xmin*0.99
            xmax=xmax*1.01
          end if
        else if (ymax.eq.ymin) then
          if (ymin.eq.0.0) then
            ymin=-1.0
            ymax=+1.0
          else
            ymin=ymin*0.99
            ymax=ymax*1.01
          end if
        end if
        if (def(jmin).eq.'*') then
c          print '(2a)','minimal and maximal values for x and y ',
c     *                 '(default: as in data file):'
          call getlin (line)
          if (line.eq.' ') write (line,'(4F13.4)') xmin,xmax,ymin,ymax
        else
c          print '(4a)','minimal and maximal values for x and y ',
c     *                 '(default: ',def(jmin)(1:ldef(jmin)),'):'
          call getlin (line)
          if (line.eq.' ') line=def(jmin)
        end if
        rewind (2)
        write (2,'(A)') line
        rewind (2)
        read (2,*) xmin,xmax,ymin,ymax
        dx=xmax-xmin
        dy=ymax-ymin
c       ---------------------------------- use rectangular display area?
        if (firplt) then
          vxmin=-1.0
          vxmax= 1.0
          vymin=-1.0
          vymax= 1.0
          call yesno ('use device dependent rectangular display area',
     *                def(jrec)(1:ldef(jrec)),yes)
          if (yes) then
            aspek=1.414
            if (rotate) aspek=1.0/aspek
            call jaspek (1,aspek)
            if (aspek.lt.1.0) then
              vymin=-aspek
              vymax= aspek
            else
              vxmin=-1.0/aspek
              vxmax= 1.0/aspek
            end if
          end if
          call jvspac (vxmin,vxmax,vymin,vymax)
        end if
c       ------------------------------------------------ select viewport
c        print '(A,4F7.3)','maximal viewport size (full display area):',
c     *        vxmin,vxmax,vymin,vymax
        if (def(jvie).eq.'*') then
c          print '(a)','viewport size (default: full display area):'
          call getlin (line)
          if (line.eq.' ')
     *      write (line,'(4F10.5)') vxmin,vxmax,vymin,vymax
        else
c          print '(3a)','viewport size (default: ',
c     *                 def(jvie)(1:ldef(jvie)),'):'
          call getlin (line)
          if (line.eq.' ') line=def(jvie)
        end if
        if (index(line,'>')+index(line,'<')+
     *      index(line,'v')+index(line,'^').gt.0) then
          illexp=.false.
          iplot=0
          l=0
          k=0
          do i=1,lenstr(line)
            if (line(i:i).eq.'>') then
              illexp=illexp .or. k.ne.0 .and. k.ne.1
              l=l+1
              k=1
            else if (line(i:i).eq.'<') then
              illexp=illexp .or. k.ne.0 .and. k.ne.2
              l=l+1
              k=2
            else if (line(i:i).eq.'v') then
              illexp=illexp .or. k.ne.0 .and. k.ne.3
              l=l+1
              k=3
            else if (line(i:i).eq.'^') then
              illexp=illexp .or. k.ne.0 .and. k.ne.4
              l=l+1
              k=4
            else if (line(i:i).ge.'A' .and. line(i:i).le.'Z') then
              illexp=illexp .or. iplot.gt.0
              iplot=ichar(line(i:i))-ichar('A')+1
              illexp=illexp .or. iplot.ge.nplot
            else
              illexp=line(i:i).gt.' '
            end if
          end do
          if (illexp)
     *      call errmsg ('GRAF: Illegal viewport specification')
          if (iplot.eq.0) iplot=nplot-1
          if (k.eq.1) then
            s=(vxmax-vxmin)*0.05*(l-1)
            v(1)=vxmaxi(iplot)+s
            v(2)=vxmaxi(iplot)+s+(vxmaxi(iplot)-vxmini(iplot))
            v(3)=vymini(iplot)
            v(4)=vymaxi(iplot)
          else if (k.eq.2) then
            s=(vxmax-vxmin)*0.05*(l-1)
            v(1)=vxmini(iplot)-s-(vxmaxi(iplot)-vxmini(iplot))
            v(2)=vxmini(iplot)-s
            v(3)=vymini(iplot)
            v(4)=vymaxi(iplot)
          else if (k.eq.3) then
            s=(vymax-vymin)*0.05*(l-1)
            v(1)=vxmini(iplot)
            v(2)=vxmaxi(iplot)
            v(3)=vymini(iplot)-s-(vymaxi(iplot)-vymini(iplot))
            v(4)=vymini(iplot)-s
          else if (k.eq.4) then
            s=(vymax-vymin)*0.05*(l-1)
            v(1)=vxmini(iplot)
            v(2)=vxmaxi(iplot)
            v(3)=vymaxi(iplot)+s
            v(4)=vymaxi(iplot)+s+(vymaxi(iplot)-vymini(iplot))
          end if
        else
          rewind (2)
          write (2,'(A)') line
          rewind (2)
          read (2,*) (v(k),k=1,4)
        end if
c        print *,v
        vxmini(nplot)=max(vxmin,v(1))
        vxmaxi(nplot)=min(vxmax,v(2))
        vymini(nplot)=max(vymin,v(3))
        vymaxi(nplot)=min(vymax,v(4))
        if (vxmini(nplot).ge.vxmaxi(nplot) .or.
     *      vymini(nplot).ge.vymaxi(nplot))
     *    call errmsg ('GRAF: Illegal viewport specification')
        call jvport (vxmini(nplot),vxmaxi(nplot),
     *               vymini(nplot),vymaxi(nplot))
        aspek=(vymaxi(nplot)-vymini(nplot))/
     *        (vxmaxi(nplot)-vxmini(nplot))
c       ------------------------------------------------- cross of axes?
        cross=.false.
        if (xmin*xmax.lt.0.0 .and. ymin*ymax.lt.0.0)
     *    call yesno ('cross of axes instead of frame',
     *                def(jcro)(1:ldef(jcro)),cross)
c       --------------------------------------------- draw axes or frame
        call jps ('%')
        call jps ('% ------ Axes or frame ------')
        call jps ('%')
        defwid=0.65
        call jlwide (defwid)
        xmaxi=xmax+dx/50.0
        ymaxi=ymax+dy/50.0
        if (cross) then
          xmini=min(xmin,-dx/10.0)
          ymini=min(ymin,-dy/15.0)
          call jmove (xmin,0.0)
          call jdraw (xmax,0.0)
          call jmove (0.0,ymin)
          call jdraw (0.0,ymax)
          yxaxis=0.0
          xyaxis=0.0
        else
          xmini=xmin-dx/10.0
          ymini=ymin-dy/15.0
          call jpidex (0,0)
          call jrect (xmin,ymin,xmax,ymax)
          yxaxis=ymin
          xyaxis=xmin
        end if
        deltax=xmaxi-xmini
        deltay=ymaxi-ymini
        dxmark=deltax/200.0
        dymark=deltay/200.0
c       ------------------------------------------------ set letter size
c        print '(4a)','letter size relative to full x range ',
c     *               '(default: ',def(jlet)(1:ldef(jlet)),'):'
        call getlin (line)
        if (line.eq.' ') line=def(jlet)
        rewind (2)
        write (2,'(A)') line
        rewind (2)
        read (2,*) z
        dxwrit=deltax*z
c        dywrit=deltay*z*1.4
        dywrit=deltay*z
        if (aspek.gt.1.0) then
          dxmark=dxmark*aspek
          dxwrit=dxwrit*aspek
        else if (aspek.lt.1.0) then
          dymark=dymark/aspek
          dywrit=dywrit/aspek
        end if
        call jsize (dxwrit,dywrit)
        call jfont (deffnt)
c        call jfont
c       ---------------------------------------- set up default x labels
        if (nxla.eq.0) then
          z=log10(dx/10.0)
          if (z.lt.0) then
            u=10.0**(aint(z)-1)
          else
            u=10.0**aint(z)
          end if
          z=dx/10.0/u
          if (z.lt.1.5) then
            z=1.0
          else if (z.lt.3.5) then
            z=2.0
          else if (z.lt.7.5) then
            z=5.0
          else
            z=10.0
          end if
          z=z*u
          i1=int(xmin/z)
          if (xmin.gt.0.0 .and. xmin.gt.z*i1*1.0001) i1=i1+1
          i2=int(xmax/z)
          if (xmax.lt.0.0 .and. xmax.lt.z*i2*1.0001) i2=i2-1
          kk=0
          do i=i1,i2
            nxla=nxla+1
            xla(nxla)=z*i
            if (mod(i,2).eq.0 .and. .not.(cross .and. i.eq.0)) then
              write (str,'(F20.5)') xla(nxla)
              do k=20,1,-1
                if (str(k:k).ne.'0') go to 418
              end do
  418         if (str(k:k).eq.'.') k=k-1
              kk=max(kk,k)
              xlatxt(nxla)=str
            else
              xlatxt(nxla)=' '
            end if
          end do
          do i=1,nxla
            do j=1,20
              if (xlatxt(i)(j:j).gt.' ') go to 430
            end do
            j=1
  430       if (xlatxt(i).ne.' ') xlatxt(i)=xlatxt(i)(j:kk)
          end do
c       ------------------------------------ set up non-default x labels
        else
          rewind (2)
          write (2,'(A)') (xlatxt(i),i=1,nxla)
          rewind (2)
          do i=1,nxla
            read (2,*) xla(i)
            do j=1,80
              if (xlatxt(i)(j:j).gt.' ') go to 522
            end do
  522       do k=j+1,80
              if (xlatxt(i)(k:k).le.' ') go to 524
            end do
  524       do j=k+1,80
              if (xlatxt(i)(j:j).gt.' ') go to 526
            end do
  526       if (j.le.80) then
              xlatxt(i)=xlatxt(i)(j:80)
            else
              xlatxt(i)=' '
            end if
          end do
        end if
c       ------------------------------------------------- write x labels
        call jps ('%')
        call jps ('% ------ labels ------')
        call jps ('%')
        call jjust (2,3)
        do i=1,nxla
          if (xla(i).ge.xmin .and. xla(i).le.xmax) then
            if (xlatxt(i).ne.' ') then
              call jmove (xla(i),yxaxis-1.5*dymark)
              call jdraw (xla(i),yxaxis+1.5*dymark)
              call jmove (xla(i),yxaxis-2.5*dymark)
              if (.not.xticso)
     *          call j1strg (xlatxt(i)(1:lenstr(xlatxt(i))))
            else
              call jmove (xla(i),yxaxis-dymark)
              call jdraw (xla(i),yxaxis+dymark)
            end if
          end if
        end do
c       ---------------------------------------- set up default y labels
        if (nyla.eq.0) then
          z=log10(dy/10.0)
          if (z.lt.0) then
            u=10.0**(aint(z)-1)
          else
            u=10.0**aint(z)
          end if
          z=dy/10.0/u
          if (z.lt.1.5) then
            z=1.0
          else if (z.lt.3.5) then
            z=2.0
          else if (z.lt.7.5) then
            z=5.0
          else
            z=10.0
          end if
          z=z*u
          i1=int(ymin/z)
          if (ymin.gt.0.0 .and. ymin.gt.z*i1*1.0001) i1=i1+1
          i2=int(ymax/z)
          if (ymax.lt.0.0 .and. ymax.lt.z*i2*1.0001) i2=i2-1
          kk=0
          do i=i1,i2
            nyla=nyla+1
            yla(nyla)=z*i
            if (mod(i,2).eq.0 .and. .not.(cross .and. i.eq.0)) then
              write (str,'(F20.5)') yla(nyla)
              do k=20,1,-1
                if (str(k:k).ne.'0') go to 537
              end do
  537         if (str(k:k).eq.'.') k=k-1
              kk=max(kk,k)
              ylatxt(nyla)=str
            else
              ylatxt(nyla)=' '
            end if
          end do
          do i=1,nyla
            do j=1,20
              if (ylatxt(i)(j:j).gt.' ') go to 540
            end do
            j=1
  540       if (ylatxt(i).ne.' ') ylatxt(i)=ylatxt(i)(j:kk)
          end do
c       ------------------------------------ set up non-default y labels
        else
          rewind (2)
          write (2,'(A)') (ylatxt(i),i=1,nyla)
          rewind (2)
          do i=1,nyla
            read (2,*) yla(i)
            do j=1,80
              if (ylatxt(i)(j:j).gt.' ') go to 542
            end do
  542       do k=j+1,80
              if (ylatxt(i)(k:k).le.' ') go to 544
            end do
  544       do j=k+1,80
              if (ylatxt(i)(j:j).gt.' ') go to 546
            end do
  546       if (j.le.80) then
              ylatxt(i)=ylatxt(i)(j:80)
            else
              ylatxt(i)=' '
            end if
          end do
        end if
c       ------------------------------------------------- write y labels
        call jjust (3,2)
        do i=1,nyla
          if (yla(i).ge.ymin .and. yla(i).le.ymax) then
            if (ylatxt(i).ne.' ') then
              call jmove (xyaxis-1.5*dxmark,yla(i))
              call jdraw (xyaxis+1.5*dxmark,yla(i))
              call jmove (xyaxis-2.5*dxmark,yla(i))
              if (.not.yticso)
     *          call j1strg (ylatxt(i)(1:lenstr(ylatxt(i))))
            else
              call jmove (xyaxis-dxmark,yla(i))
              call jdraw (xyaxis+dxmark,yla(i))
            end if
          end if
        end do
c       -------------------------------------------------- text for axes
c        print '(3a)','text for x-axis (default: "',
c     *               def(jxtx)(1:ldef(jxtx)),'"):'
        call getlin (line)
        if (line.eq.' ') line=def(jxtx)
        if (line.ne.' ') then
          call jjust (3,3)
          call jmove (xmax-dx*0.05,yxaxis-5.0*dymark-dywrit)
          do i=80,1,-1
            if (line(i:i).gt.' ') go to 630
          end do
  630     continue
          call j1strg (line(1:i))
        end if
c        print '(3a)','text for y-axis (default: "',
c     *               def(jytx)(1:ldef(jytx)),'"):'
        call getlin (line)
        if (line.eq.' ') line=def(jytx)
        if (line.ne.' ') then
          call jjust (3,2)
          yt=ymax-dy*0.05
          if (abs(yt-yy).lt.1.2*dywrit) yt=yy-1.6*dywrit
          call jmove (xyaxis-5.0*dxmark-dxwrit,yt)
          do i=80,1,-1
            if (line(i:i).gt.' ') go to 660
          end do
  660     continue
          call j1strg (line(1:i))
        end if
      end if
c     ----------------------------------------------------- special text
      if (ntxt.gt.0) then
        call jps ('%')
        call jps ('% ------ text ------')
        call jps ('%')
        rewind (2)
        write (2,'(A)') (text(i),i=1,ntxt)
        rewind (2)
        do i=1,ntxt
          read (2,*) just,xx,yy,size
          call jsize (dxwrit*size,dywrit*size)
          call jjust (abs(just),2)
          if (just.lt.0) then
            xx=xmin+xx*(xmax-xmin)
            yy=ymin+yy*(ymax-ymin)
          end if
          call jmove (xx,yy)
          k=0
          do j=2,80
            if (text(i)(j-1:j-1).le.' ' .and. text(i)(j:j).gt.' ')
     *      then
              k=k+1
              if (k.eq.4) go to 690
            end if
          end do
          go to 700
  690     continue
          call j1strg (text(i)(j:lenstr(text(i))))
  700     continue
        end do
        ntxt=0
      end if
c     ---------------------------------------------- set representations
c      print '(a)','possible data point representations:',
c     *  '  - ... nothing',
c     *  '  . ... small point        p ... point',
c     *  '  + ... +-shaped cross     x ... x-shaped cross',
c     *  '  o ... small circle       d ... filled circle (dot)',
c     *  '  s ... small square       f ... filled square',
c     *  '  ^ ... rotated square',
c     *  '  e ... error bar (2 col.) y ... symmetric error bar (2 col.)',
c     *  '  r ... range (2 col.)     w ... filled range (2 columns)',
c     *  '  c ... connected squares (2 columns)',
c     *  '  h ... histogram, not hatched',
c     *  '  v ... histogram, vertically hatched',
c     *  '  = ... histogram, horizontally hatched',
c     *  '  / ... histogram, //// hatched',
c     *  '  b ... histogram, filled'
c      print '(3a)','data point representations (default: ',
c     *             def(jrep)(1:ldef(jrep)),'):'
      call getlin (line)
      if (line.eq.' ') line=def(jrep)
      do i=1,80
        if (line(i:i).gt.' ') go to 740
      end do
  740 repres=line(i:80)
      do i=1,ncrv
        if (repres(i:i).eq.' ') repres(i:i)=repres(i-1:i-1)
        if (repres(i:i).ge.'a' .and. repres(i:i).le.'z')
     *    repres(i:i)=char(ichar(repres(i:i))+ichar('A')-ichar('a'))
        if (index('-&.+XOSFCHV=/BRPD^EWY\\',repres(i:i)).eq.0)
     *    call errmsg ('GRAF: Illegal representation type')
      end do
c     -------------------------------------------------- set connections
c      print '(a)','possible data point connections:',
c     *      '  - ... nothing',
c     *      '  l ... lines   h ... heavy lines   t ... thin lines',
c     *      '  s ... spline  b ... bold spline   f ... fine spline'
c      print '(3a)','data point connections (default: ',
c     *             def(jcon)(1:ldef(jcon)),'):'
      call getlin (line)
      if (line.eq.' ') line=def(jcon)
      do i=1,80
        if (line(i:i).gt.' ') go to 800
      end do
  800 conect=line(i:80)
      nocnct=.true.
      do i=1,ncrv
        if (conect(i:i).eq.' ') conect(i:i)=conect(i-1:i-1)
        if (conect(i:i).ge.'a' .and. conect(i:i).le.'z')
     *    conect(i:i)=char(ichar(conect(i:i))+ichar('A')-ichar('a'))
        if (index('-LHTSBF',conect(i:i)).eq.0)
     *    call errmsg ('GRAF: Illegal connection type')
        nocnct=(nocnct .and. conect(i:i).eq.'-')
      end do
c     -------------------------------------------------- set line styles
      if (.not.nocnct) then
c        print '(4a)','line styles for connections (0[solid]..9, ',
c     *               'default: ',def(jlin)(1:ldef(jlin)),'):'
        call getlin (line)
        if (line.eq.' ') line=def(jlin)
        do i=1,80
          if (line(i:i).gt.' ') go to 823
        end do
  823   linstl=line(i:80)
        do i=1,ncrv
          if (linstl(i:i).eq.' ') linstl(i:i)=linstl(i-1:i-1)
          if (linstl(i:i).lt.'0' .or. linstl(i:i).gt.'9')
     *      call errmsg ('GRAF: Illegal line style')
        end do
      end if
c     ------------------------------------------------- draw data points
      do i=1,ncrv
        if (repres(i:i).ne.'-' .and. repres(i:i).ne.'&') then
          write (line,'(I4)') i
          call lefstr (line)
          call jps ('%')
          call jps ('% ------ curve '//line(1:lenstr(line))//
     *              ' (representatation '//repres(i:i)// ') ------')
          call jps ('%')
          if (repres(i:i).eq.'P') call jlwide (2.5)
c          if (repres(i:i).eq.'D') call jlwide (5.0)
          do j=1,n
            if (x(j).ge.xmin .and. x(j).le.xmax .and.
     *          ((repres(i:i).ne.'Y' .and.
     *            y(j,i).ge.ymin .and. y(j,i).le.ymax) .or.
     *           (repres(i:i).eq.'Y' .and.
     *            y(j,i-1).ge.ymin .and. y(j,i-1).le.ymax))) then
              if (index('.P',repres(i:i)).gt.0) then
                call jmark (x(j),y(j,i))
              else if (repres(i:i).eq.'+') then
                call jmove (x(j)-dxmark,y(j,i))
                call jdraw (x(j)+dxmark,y(j,i))
                call jmove (x(j),y(j,i)-dymark)
                call jdraw (x(j),y(j,i)+dymark)
              else if (repres(i:i).eq.'X') then
                call jmove (x(j)-dxmark,y(j,i)-dymark)
                call jdraw (x(j)+dxmark,y(j,i)+dymark)
                call jmove (x(j)-dxmark,y(j,i)+dymark)
                call jdraw (x(j)+dxmark,y(j,i)-dymark)
              else if (repres(i:i).eq.'^') then
                call jmove (x(j)-dxmark*1.414,y(j,i))
                call jdraw (x(j),y(j,i)-dymark*1.414)
                call jdraw (x(j)+dxmark*1.414,y(j,i))
                call jdraw (x(j),y(j,i)+dymark*1.414)
                call jdraw (x(j)-dxmark*1.414,y(j,i))
              else if (repres(i:i).eq.'O') then
                call jcirc (x(j),y(j,i),dxmark)
              else if (repres(i:i).eq.'D') then
                call jpidex (0,47)
                call jcirc (x(j),y(j,i),dxmark)
              else if (repres(i:i).eq.'S') then
                call jrect (x(j)-dxmark,y(j,i)-dymark,
     *                      x(j)+dxmark,y(j,i)+dymark)
              else if (repres(i:i).eq.'F') then
                call jpidex (0,47)
                call jrect (x(j)-dxmark,y(j,i)-dymark,
     *                      x(j)+dxmark,y(j,i)+dymark)
              else if (repres(i:i).eq.'C') then
                if (y(j,i).gt.y(j,i-1)) then
                  call jmove (x(j),y(j,i)-dymark)
                  call jdraw (x(j),y(j,i-1)+dymark)
                else
                  call jmove (x(j),y(j,i)+dymark)
                  call jdraw (x(j),y(j,i-1)-dymark)
                end if
                call jpidex (0,0)
                call jrect (x(j)-dxmark,y(j,i)-dymark,
     *                      x(j)+dxmark,y(j,i)+dymark)
                call jpidex (0,47)
                call jrect (x(j)-dxmark,y(j,i-1)-dymark,
     *                      x(j)+dxmark,y(j,i-1)+dymark)
              else if (repres(i:i).eq.'E') then
                call jmove (x(j),y(j,i))
                call jdraw (x(j),y(j,i-1))
                call jmove (x(j)-dxmark,y(j,i))
                call jdraw (x(j)+dxmark,y(j,i))
                call jmove (x(j)-dxmark,y(j,i-1))
                call jdraw (x(j)+dxmark,y(j,i-1))
              else if (repres(i:i).eq.'Y') then
                call jmove (x(j),y(j,i-1)+y(j,i))
                call jdraw (x(j),y(j,i-1)-y(j,i))
                call jmove (x(j)-dxmark,y(j,i-1)+y(j,i))
                call jdraw (x(j)+dxmark,y(j,i-1)+y(j,i))
                call jmove (x(j)-dxmark,y(j,i-1)-y(j,i))
                call jdraw (x(j)+dxmark,y(j,i-1)-y(j,i))
              else if (index('HV=/BRW\\',repres(i:i)).gt.0) then
                if (j.eq.1) then
                  x1=xmin
                else
                  x1=0.5*(x(j-1)+x(j))
                end if
                if (j.eq.n) then
                  x2=xmax
                else
                  x2=0.5*(x(j)+x(j+1))
                end if
                if (repres(max(1,i-1):max(1,i-1)).eq.'&') then
                  y1=min(y(j,i-1),y(j,i))
                  y2=max(y(j,i-1),y(j,i))
                else
                  y1=min(0.0,y(j,i))
                  y2=max(0.0,y(j,i))
                end if
                if (repres(i:i).eq.'R') then
                  call jrect (x1,y(j,i-1),x2,y(j,i))
                else if (repres(i:i).eq.'W') then
                  call jpidex (0,47)
                  call jrect (x1,y(j,i-1),x2,y(j,i))
                else if (repres(i:i).eq.'B') then
                  call jpidex (0,47)
                  call jrect (x1,y1,x2,y2)
                else
                  call jrect (x1,y1,x2,y2)
                end if
                if (repres(i:i).eq.'V') then
                  z=x2-x1
                  k=max(2,nint(z/dxmark))
                  do l=1,k-1
                    xx=x1+l*z/k
                    call jmove (xx,y1)
                    call jdraw (xx,y2)
                  end do
                else if (repres(i:i).eq.'=') then
                  k=max(2,nint(abs(y2-y1)/dymark))
                  do l=1,k-1
                    yy=y1+l*(y2-y1)/k
                    call jmove (x1,yy)
                    call jdraw (x2,yy)
                  end do
                else if (repres(i:i).eq.'/') then
                  call hatch (x1,y1,x2,y2,
     *                        deltay/deltax/aspek,2.0*dymark)
                else if (repres(i:i).eq.'\\') then
                  call hatch (x1,y1,x2,y2,
     *                        -deltay/deltax/aspek,2.0*dymark)
                end if
              end if
            end if
          end do
          call jlwide (defwid)
          call jpidex (0,0)
        end if
      end do
c     ------------------------------------------------- draw connections
      warned=.false.
      do i=1,ncrv
        if (conect(i:i).ne.'-') then
          write (line,'(I4)') i
          call lefstr (line)
          call jps ('%')
          call jps ('% ------ curve '//line(1:lenstr(line))//
     *              ' (connection '//conect(i:i)// ') ------')
          call jps ('%')
        end if
        if (index('SBF',conect(i:i)).gt.0) then
          if (conect(i:i).eq.'B') call jlwide (2.0*defwid)
          if (conect(i:i).eq.'F') call jlwide (0.0)
          call jlstyl (ichar(linstl(i:i))-ichar('0'))
          nsplin=0
          i1=n+1
          i2=0
          do j=1,n
            if (x(j).ge.xmin .and. x(j).le.xmax) then
              i1=min(i1,j)
              i2=max(i2,j)
            end if
            nsplin=nsplin+1
            xsplin(nsplin)=x(j)
            ysplin(nsplin)=y(j,i)
            if (j.gt.1) then
              if (x(j).le.x(j-1)) nsplin=nsplin-1
            end if
          end do
          if (nsplin.lt.n .and. .not. warned) then
            write (*,'(A,I4,2A)') '% Warning:',n-nsplin,' points with ',
     *            'non-unique x-values skipped for spline interpolation'
            warned=.true.
          end if
          if (i1.lt.i2) then
            call spline (xsplin,ysplin,nsplin,1.0e30,1.0e30,ytmp)
            call splint (xsplin,ysplin,ytmp,nsplin,x(i1),yy)
            call jmove (x(i1),yy)
            xxprev=x(i1)
            yyprev=yy
            dxx=(x(i2)-x(i1))/200.1
            do xx=x(i1)+dxx,x(i2),dxx
              call splint (xsplin,ysplin,ytmp,nsplin,xx,yy)
              if (yy.gt.ymax .and. yyprev.lt.ymax) then
                call jdraw (xxprev+dxx*(ymax-yyprev)/(yy-yyprev),ymax)
              else if (yy.lt.ymax .and. yyprev.gt.ymax) then
                call jmove (xxprev+dxx*(ymax-yyprev)/(yy-yyprev),ymax)
                call jdraw (xx,yy)
              else if (yy.lt.ymin .and. yyprev.gt.ymin) then
                call jdraw (xxprev+dxx*(ymin-yyprev)/(yy-yyprev),ymin)
              else if (yy.gt.ymin .and. yyprev.lt.ymin) then
                call jmove (xxprev+dxx*(ymin-yyprev)/(yy-yyprev),ymin)
               call jdraw (xx,yy)
              else if (yy.ge.ymin .and. yy.le.ymax) then
                call jdraw (xx,yy)
              end if
              xxprev=xx
              yyprev=yy
            end do
            call splint (xsplin,ysplin,ytmp,nsplin,x(i2),yy)
            if (yy.ge.ymin .and. yy.le.ymax) call jdraw (x(i2),yy)
          end if
        else if (index('LHT',conect(i:i)).gt.0) then
          if (conect(i:i).eq.'H') call jlwide (2.0*defwid)
          if (conect(i:i).eq.'T') call jlwide (0.0)
          call jlstyl (ichar(linstl(i:i))-ichar('0'))
          previn=.false.
          do j=1,n
            inside=x(j).ge.xmin .and. x(j).le.xmax .and.
     *             y(j,i).ge.ymin .and. y(j,i).le.ymax
            if (previn) then
              if (inside) then
                call jdraw (x(j),y(j,i))
                if (mod(j,500).eq.0) call jmove (x(j),y(j,i))
              end if
            else
              if (inside) call jmove (x(j),y(j,i))
            end if
            previn=inside
          end do
        end if
        call jlwide (defwid)
      end do
      call jlstyl (0)
c     ------------------------------------------- another curve or plot?
      if (nxtcur) then
        fircur=.false.
        call jps ('%')
        call jps ('% ------ next set of curves ------')
        call jps ('%')
        go to 220
      else if (nxtplt) then
        firplt=.false.
        fircur=.true.
        readit=.not.pagain
        call jps ('%')
        call jps ('% ------ next plot ------')
        call jps ('%')
        go to 220
      else
        firfil=.false.
c        go to 150
      end if
c     ------------------------------------------------------------------
  930 continue
      call jend
c      close (11)
      close (2)
      end
c     ==================================================================
c               Read line from standard input, if necessary.
c
      subroutine getlin (line)
c
      character*(*) line
      logical cnfdef
      common /defdat/cnfdef
c
      line=' '
      if (cnfdef) read (*,'(A)') line
      end
c     ==================================================================
      subroutine readfn (text,defext,defnam,filnam)
c
      character*(*) text,filnam,defnam,defext
c
      if (defnam.eq.' ') then
c        print '(4A)',text,' (',defext(2:),'):'
      else
c        print '(6A)',text,' (',defext(2:),', default: ',defnam,'):'
      end if
      read (*,'(A)') filnam
      i=index(filnam,'#')
      if (i.gt.0) filnam(i:)=' '
      if (filnam.eq.' ') then
        if (defnam.eq.' ' .or. defnam.eq.'none') return
        filnam=defnam
      end if
      l=lenstr(filnam)
      do i=1,l
        if (filnam(i:i).gt.' ') go to 20
      end do
   20 do j=i,l
        filnam(j-i+1:j-i+1)=filnam(j:j)
      end do
      l=l-i+1
      if (filnam(l:l).eq.'|') then
        filnam(l:l)=' '
        return
      end if
      if (filnam.eq.'*') return
      do i=l,1,-1
        if (filnam(i:i).eq.defext(1:1)) go to 40
      end do
   40 do j=l,1,-1
        if (index(':]/',filnam(j:j)).gt.0) go to 60
      end do
   60 if (i.le.j) filnam(l+1:)=defext
      end
c     ==================================================================
c     YESNO:    Get yes/no answer from standard input.
c
c               On input, TEXT is a character string that is written
c               to standard output, followed by an indication of the
c               default answer. This default has to be specified in
c               the character string DEFLT: if the first character of
c               this string equals 'y' or 'Y' the default value for
c               the logical variable ANSWER will be .TRUE. otherwise
c               it will be .FALSE.. On output, ANSWER will have this
c               default value if a blank line is read from standard
c               input. If a non-blank line is read from standard in-
c               put the value of ANSWER on output is determined by
c               the first non-blank character on this line: 'y' or
c               'Y' yield .TRUE., and 'n' or 'N' yield .FALSE., re-
c               spectively; otherwise the default (according to DEFLT)
c               is used.
c
c               Written in standard FORTRAN-77.
c               Peter G"untert, 25-4-1989
c     ------------------------------------------------------------------
      subroutine yesno (text,deflt,answer)
c
      character*(*) text,deflt
      logical answer
      character*80 line
c
      answer=deflt(1:1).eq.'Y' .or. deflt(1:1).eq.'y'
c      print '(4a)',text,' (default: ',deflt,')?'
      call getlin (line)
      do i=1,80
        if (line(i:i).gt.' ') then
          if (answer) then
            answer=.not.(line(i:i).eq.'N' .or. line(i:i).eq.'n')
          else
            answer=line(i:i).eq.'Y' .or. line(i:i).eq.'y'
          end if
          return
        end if
      end do
      end
c     ==================================================================
c     FULNAM:   Determine full name by comparing a (possibly)
c               truncated name with a list of all possible
c               names.
c
c               NAM is input as the truncated name, the array
c               NAMLST as the list of all possible names, N
c               as the number of entries in this list. On out-
c               put, the follwing cases may occur:
c
c                                   INAM  NAM
c                 unambiguous name:   >0  replaced by full name
c                 name not found  :   -1  unchanged
c                 ambiguous name  :   -2  unchanged
c
c               INAM is set to the array index of NAMLST if an
c               unambiguous name is found.
c
c               Written in standard FORTRAN-77.
c               Peter G"untert, 28-11-1988
c     ------------------------------------------------------------------
      subroutine fulnam (nam,namlst,n,inam)
c
      character*(*) nam,namlst(n)
      logical match
c
      l=len(nam)
      lc=lenstr(nam)
      if (lc.eq.0) then
        if (n.eq.1) then
          inam=1
          nam=namlst(1)
          return
        else
          inam=-2
          return
        end if
      end if
      if (nam(l:l).le.' ') nam=nam(1:lc)//'*'
      inam=0
      do i=1,n
        if (match(namlst(i),nam)) then
          if (inam.eq.0) then
            inam=i
          else
            inam=-2
            nam=nam(1:lc)
            return
          end if
        end if
      end do
      if (inam.eq.0) then
        inam=-1
        nam=nam(1:lc)
      else
        nam=namlst(inam)
      end if
      end
c     ==================================================================
c     ERRMSG:   Display error messages and stop.
c
c               Written in standard FORTRAN-77.
c               Peter G"untert, 20-1-1988
c     ------------------------------------------------------------------
      subroutine errmsg (text)
c
      character*(*) text
c
      print '(3a)','ERROR in ',text(1:lenstr(text)),'.'
      stop
      end
c     ==================================================================
c     LENSTR:   Return index of last non-blank character of the
c               string S, or zero if S is a blank string.
c
c               Written in standard FORTRAN-77.
c               Peter G"untert, 28-11-1988
c     ------------------------------------------------------------------
      function lenstr(s)
c
      character*(*) s
c
      do l=len(s),1,-1
        if (s(l:l).gt.' ') go to 20
      end do
   20 lenstr=l
      end
c     ==================================================================
      subroutine lefstr (s)
      character*(*) s
      logical flag
      if (s(1:1).gt.' ' .or. s.eq.' ') return
      i=0
      flag=.false.
      do l=1,len(s)
        flag=flag .or. s(l:l).gt.' '
        if (flag) then
          i=i+1
          s(i:i)=s(l:l)
        end if
      end do
      s(i+1:)=' '
      end
c     ==================================================================
c     MATCH:    Does the string S match the string O which may
c               contain wildcards?
c
c               The two possible wildcards are:
c                 "%" ... matches exactly one character
c                 "*" ... matches any number (including zero)
c                         of characters
c               The two strings S and O will not be changed.
c
c               Written in standard FORTRAN-77.
c               Peter Guentert, 28-11-1988
c     ------------------------------------------------------------------
      logical function match(s,o)
c
      character*(*) s,o
      logical       fix
c
      fix=.true.
      io=lenstr(o)
      is=lenstr(s)
      if (io.eq.0) then
        match=.true.
        return
      else if (is.eq.0) then
        match=.true.
        return
      end if
      jo=1
      js=1
   10 if (jo.le.io) then
        if (o(jo:jo).eq.'*') then
          fix=.false.
          jo=jo+1
        else
          j=index(o(jo:io),'*')-1
          if (j.lt.0) j=io-jo+1
          if (fix) then
            imax=js
          else
            imax=is-j+1
          end if
          do i=js,imax
            do l=jo,jo+j-1
              if (o(l:l).ne.'%') then
                k=i+l-jo
                if (o(l:l).ne.s(k:k)) go to 30
              end if
            end do
            jo=jo+j
            js=i+j
            fix=.true.
            go to 10
   30       continue
          end do
          match=.false.
          return
        end if
      go to 10
      end if
      match=(js.gt.is .or. (.not. fix))
      end
c     ==================================================================
c               subroutines for spline interpolation
c               (from "Numerical Recipes")
c
      subroutine spline (x,y,n,yp1,ypn,y2)
c
      parameter (nmax=8200)
      dimension x(n),y(n),y2(n),u(nmax)
c
      if (yp1.gt.0.99e30) then
        y2(1)=0.0
        u(1)=0.0
      else
        y2(1)=-0.5
        u(1)=(3.0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      end if
      do i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.0
        y2(i)=(sig-1.0)/p
        u(i)=(6.0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
     *       /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
      end do
      if (ypn.gt.0.99e30) then
        qn=0.0
        un=0.0
      else
        qn=0.5
        un=(3.0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      end if
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.0)
      do k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
      end do
      end
c     ==================================================================
      subroutine splint (xa,ya,y2a,n,x,y)
c
      dimension xa(n),ya(n),y2a(n)
c
      klo=1
      khi=n
    1 if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if (xa(k).gt.x) then
          khi=k
        else
          klo=k
        end if
        go to 1
      end if
      h=xa(khi)-xa(klo)
      if (h.eq.0.0) call errmsg ('SPLINT: Bad XA input')
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+
     *  ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6
      end
c     ==================================================================
c               quicksort algorithm
c
      subroutine rsort (x,indx,n)
c
      parameter (nstack=30)
      dimension x(n),indx(n),isl(nstack),isr(nstack)
c
      do i=1,n
        indx(i)=i
      end do
      js=1
      isl(1)=1
      isr(1)=n
   50 continue
        jl=isl(js)
        jr=isr(js)
        js=js-1
   52   continue
          i=jl
          j=jr
          xx=x((jl+jr)/2)
   53     continue
   54       if (x(i).lt.xx) then
              i=i+1
            go to 54
            end if
   55       if (x(j).gt.xx) then
              j=j-1
            go to 55
            end if
            if (i.gt.j) go to 58
            t=x(i)
            x(i)=x(j)
            x(j)=t
            k=indx(i)
            indx(i)=indx(j)
            indx(j)=k
            i=i+1
            j=j-1
          if (i.le.j) go to 53
   58     if (j-jl.lt.jr-i) then
            if (i.lt.jr) then
              js=js+1
              isl(js)=i
              isr(js)=jr
            end if
            jr=j
          else
            if (jl.lt.j) then
              js=js+1
              isl(js)=jl
              isr(js)=j
            end if
            jl=i
          end if
        if (jl.lt.jr) go to 52
      if (js.gt.0) go to 50
      end
c     ==================================================================
c               hatch rectangular area
c
c               (x1,y1)   lower left edge of rectangle
c               (x2,y2)   upper right edge of rectangle (x2>x1, y2>y1)
c               slope     slope of hatching lines
c               deltab    vertical distance between hatching lines,
c                         deltab>0
c
      subroutine hatch (x1,y1,x2,y2,slope,deltab)
c
      if (slope.gt.0.0) then
        bll=y1-slope*x1
        bur=y2-slope*x2
        do b=y1-slope*x2,y2-slope*x1,deltab
          if (b.le.bll) then
            call jmove ((y1-b)/slope,y1)
         else
            call jmove (x1,slope*x1+b)
          end if
          if (b.le.bur) then
            call jdraw (x2,slope*x2+b)
          else
            call jdraw ((y2-b)/slope,y2)
          end if
        end do
      else
        bll=y1-slope*x2
        bur=y2-slope*x1
        do b=y1-slope*x1,y2-slope*x2,deltab
          if (b.le.bll) then
            call jmove ((y1-b)/slope,y1)
         else
            call jmove (x2,slope*x2+b)
          end if
          if (b.le.bur) then
            call jdraw (x1,slope*x1+b)
          else
            call jdraw ((y2-b)/slope,y2)
          end if
        end do
      end if
      end
cbegin Tai-he Xia
c     ==================================================================
      subroutine jbegin
c
      parameter (maxfnt=10)
      character*80 font,fonts(maxfnt),deffnt
      common/text/font,fonts,deffnt
      common/view/dixmin,dixmax,diymin,diymax
c
      dixmin=-1.0
      dixmax= 1.0
      diymin=-1.0
      diymax= 1.0
      deffnt='Times-Roman'
      font=' '
      do i=1,maxfnt
        fonts(i)=' '
      end do
      write (*,'(a)') '%!PS-Adobe-3.0 EPSF-3.0',
     *      '%%Creator: Graf','%%Title: Plot','%%Pages: 1',
     *      '%%DocumentFonts: (atend)',
     *      '%%BoundingBox:   40   40  540  790','%%EndComments',
     *      '%%BeginProlog','%','% ------ Definitions ------','%',
     *      '/grafdict 20 dict def','grafdict begin',
     *      '/m {stroke moveto} def',
     *      '/l {lineto} def',
     *      '/s {stroke} def',
     *      '/ff {findfont} def',
     *      '/mf {makefont setfont} def',
     *      '/lw {stroke setlinewidth} def'
      write (*,'(a)')
     *      '/t {currentpoint stroke moveto',
     *      '    dup stringwidth pop fy 0.7 mul justv mul exch',
     *      '    justh mul exch rmoveto show} def',
     *      '/p {stroke 2 copy moveto currentlinewidth 0 le',
     *      '    {0.4} {currentlinewidth 0.8 mul} ifelse',
     *      '    0 360 arc fill} def',
     *      '/c {stroke 3 copy pop moveto dup 0 rmoveto 0 360 arc',
     *      '    drawmode} def',
     *      '/r {stroke 4 copy moveto pop exch lineto 3 copy pop',
     *      '    lineto exch lineto pop closepath drawmode} def',
     *      '/drawmode {stroke} def',
     *      'end','%%EndProlog'
      write (*,'(a)') '%%BeginSetup',
     *      '%','% ------ Initializations ------','%',
     *      '0 setlinewidth 0 setlinecap 2 setlinejoin',
     *      '[1 0] 0 setdash 0.0 setgray',
     *      '%%EndSetup','%%Page: 1 1',
     *      '/sobj save def','grafdict begin'
      end
c     ==================================================================
      subroutine jend
c
      parameter (maxfnt=10)
      character*80 font,fonts(maxfnt),deffnt
      common/text/font,fonts,deffnt
c
      write (*,'(a)') '%','% ------ Trailer ------','%',
     *      's','end','sobj restore','showpage','%%Trailer'
      do i=1,maxfnt
        if (fonts(i).eq.' ') go to 20
      end do
   20 write (*,'(30a)') '%%DocumentFonts:',
     *      (' ',fonts(k)(1:lenstr(fonts(k))),k=1,i-1)
      end
c     ==================================================================
      subroutine jps (line)
      character*(*) line
      write (*,'(a)') line
      end
c     ==================================================================
      subroutine jrect (x1,y1,x2,y2)
c
      cx1=calx(x1)
      cy1=caly(y1)
      cx2=calx(x2)
      cy2=caly(y2)
      if (min(cx1,cy1,cx2,cy2).gt.-99.99) then
        write (*,'(f5.1,3f6.1,'' r'')') cx1,cy1,cx2,cy2
      else
        write (*,'(f6.1,3f7.1,'' r'')') cx1,cy1,cx2,cy2
      end if
      end
c     ==================================================================
      subroutine jcirc (x,y,r)
c
      cx=calx(x)
      cy=caly(y)
      cr=calx(x+r)-cx
      if (min(cx,cy,cr).gt.-99.99) then
        write (*,'(f5.1,f6.1,f5.1,'' c'')') cx,cy,cr
      else
        write (*,'(f6.1,f7.1,f5.1,'' c'')') cx,cy,cr
      end if
      end
c     ==================================================================
      subroutine jmove (x,y)
      cx=calx(x)
      cy=caly(y)
      if (cx.gt.-99.99 .and. cy.gt.-99.99) then
        write (*,'(f5.1,f6.1,'' m'')') cx,cy
      else
        write (*,'(f6.1,f7.1,'' m'')') cx,cy
      end if
      end
c     ==================================================================
      subroutine jdraw (x,y)
      cx=calx(x)
      cy=caly(y)
      if (cx.gt.-99.99 .and. cy.gt.-99.99) then
        write (*,'(f5.1,f6.1,'' l'')') cx,cy
      else
        write (*,'(f6.1,f7.1,'' l'')') cx,cy
      end if
      end
c     ==================================================================
      subroutine jsize (x,y)
c
      common/fxy/fx,fy
c
      fx=calx(2*x)-calx(x)
      fy=caly(2*y)-caly(y)
      end
c     ==================================================================
      subroutine jfont (dfont)
c
      character*(*) dfont
c
      parameter (maxfnt=10)
      character*80 font,fonts(maxfnt),deffnt
      common/text/font,fonts,deffnt
c
      deffnt=dfont
      end
c     ==================================================================
      subroutine j1strg (str)
c
      character*(*) str,oldfnt*50
c
      parameter (maxfnt=10)
      character*80 font,fonts(maxfnt),deffnt
      common/text/font,fonts,deffnt
      common/fxy/fx,fy
      save fxp,fyp
      data fxp,fyp/-1.0,-1.0/
c
      k=0
      m=len(str)
      oldfnt=font
      if (str(1:1) .eq. '\\' .and. str(2:2).ne.'\\') then
        do k=2,m
          if (str(k:k).le.' ') go to 20
        end do
   20   font=str(2:k-1)
      else
        font=deffnt
      end if
c      print *,'deffnt=',deffnt
c      print *,'font=',font
      if (font.ne.oldfnt .or. fxp.ne.fx .or. fyp.ne.fy) then
        if (fx.ne.fxp) write (*,'(a,f6.2,a)') '/fx {',fx,'} def'
        if (fy.ne.fyp) write (*,'(a,f6.2,a)') '/fy {',fy,'} def'
        write (*,'(3a)')
     *        '/',font(1:lenstr(font)),' ff [fx 0 0 fy 0 0] mf'
        do i=1,maxfnt
          if (fonts(i).eq.' ') fonts(i)=font
          if (fonts(i).eq.font) go to 40
        end do
   40   fxp=fx
        fyp=fy
      end if
      if (k+1.le.m) write (*,'(3a)') '(',str(k+1:m),') t'
      end
c     ==================================================================
      subroutine jmark (x,y)
      cx=calx(x)
      cy=caly(y)
      if (cx.gt.-99.99 .and. cy.gt.-99.99) then
        write (*,'(f5.1,f6.1,'' p'')') cx,cy
      else
        write (*,'(f6.1,f7.1,'' p'')') cx,cy
      end if
      end
c     ==================================================================
      real function calx(x)
      common/transf/xmin,xmax,ymin,ymax,pxmin,pxmax,pymin,pymax
      calx=(x-xmin)/(xmax-xmin)*(pxmax-pxmin)+pxmin
      end
c     ==================================================================
      real function caly(y)
      common/transf/xmin,xmax,ymin,ymax,pxmin,pxmax,pymin,pymax
      caly=(y-ymin)/(ymax-ymin)*(pymax-pymin)+pymin
      end
c     ==================================================================
      subroutine jvport (vxmini,vxmaxi,vymini,vymaxi)
c
      parameter (inch=72,cm=1.0/2.54*inch)
      common/transf/xmin,xmax,ymin,ymax,pxmin,pxmax,pymin,pymax
      common/view/dixmin,dixmax,diymin,diymax
c
      if (dixmin.eq.-1.0 .and. dixmax.eq.1.0 .and.
     *    diymin.eq.-1.0 .and. diymax.eq.1.0) then
        psxmin=2.0*cm
        psxmax=19.0*cm
        psymin=(15.0+diymin*8.5)*cm
        psymax=(15.0+diymax*8.5)*cm
      else if (dixmin.eq.-1.0 .and. dixmax.eq.1.0) then
        psxmin=3.0*cm
        psxmax=27.0*cm
        psymin=(10.5+diymin*12.0)*cm
        psymax=(10.5+diymax*12.0)*cm
      else
        psxmin=(10.5+dixmin*12.0)*cm
        psxmax=(10.5+dixmax*12.0)*cm
        psymin=3.0*cm
        psymax=27.0*cm
      end if
      pxmin=(vxmini-dixmin)/(dixmax-dixmin)*(psxmax-psxmin)+psxmin
      pxmax=(vxmaxi-dixmin)/(dixmax-dixmin)*(psxmax-psxmin)+psxmin
      pymin=(vymini-diymin)/(diymax-diymin)*(psymax-psymin)+psymin
      pymax=(vymaxi-diymin)/(diymax-diymin)*(psymax-psymin)+psymin
      end
c     ==================================================================
      subroutine jlstyl (k)
c
      character*20 pat(0:9)
c
c     pattern can be changed by defining alternatively on off ps units
c     e.g., '5 2' means 5 units on 2 units off, ' ' means sodid line.
c
      data (pat(i),i=0,9)/' ','5 4','1','5 2 1 2','3 3','8 4',
     *                    '8 2 1 2 1 2 1 2','8 3 3 3','2 2',
     *                    '8 2 1 2 1 2'/
c
      write (*,'(3a)') 's [',pat(k)(1:lenstr(pat(k))),
     *                  '] 0 setdash'
      end
c     ==================================================================
      subroutine jlwide (width)
      save wid
      data wid/-1.0/
      if (width.ne.wid) then
        write (*,'(F5.2,'' lw'')') width
        wid=width
      end if
      end
c     ==================================================================
      subroutine jaspek (idev,aspek)
c     A4-specific!
      if (aspek.lt.1.0) write (*,'(''595 0 translate 90 rotate'')')
      end
c     ==================================================================
      subroutine jvspac (vxmin,vxmax,vymin,vymax)
c
      common/view/dixmin,dixmax,diymin,diymax
c
      dixmin=vxmin
      dixmax=vxmax
      diymin=vymin
      diymax=vymax
      call jvport (vxmin,vxmax,vymin,vymax)
      end
c     ==================================================================
      subroutine jjust (justh,justv)
      save justx,justy
      data justx,justy/-1,-1/
c
      if (justh.ne.justx) then
        justx=justh
        write (*,'(''/justh {'',f4.1,''} def'')') -0.5*(justh-1)
      end if
      if (justv.ne.justy) then
        justy=justv
        write (*,'(''/justv {'',f4.1,''} def'')') -0.5*(justv-1)
      end if
      end
c     ==================================================================
      subroutine jpidex (icol,intens)
      logical black
      save black
      data black/.false./
c
      if (intens.ne.0 .and. .not.black) then
        black=.true.
        write (*,'(''/drawmode {fill} def'')')
      else if (intens.eq.0 .and. black) then
        black=.false.
        write (*,'(''/drawmode {stroke} def'')')
      end if
      end