c***********************************************************************
c                                                                      *
c    Program    S Y S T E R P L O T    version  0.9      Aug 24, 2000  *
c                                                                      *
c***********************************************************************
c                                                                      *
c    J.M.M. Smits  and  R. de Gelder                                   *
c    Department of Inorganic Chemistry                                 *
c    University of Nijmegen                                            *
c    Toernooiveld 1                                                    *
c    6525 ED  Nijmegen                                                 *
c    The Netherlands                                                   *
c    smits@sci.kun.nl   and   rdg@sci.kun.nl                           *
c                                                                      *
c    This program is the property of the abovementioned authors.       *
c    However, herewith permission is granted to users to adapt the     *
c    program to their local platform as far as necessary.              *
c    Proper credits are due, though.                                   *
c                                                                      *
c    The program uses OpenGL and GLUT with its FORTRAN bindings for    *
c    which the authors cannot claim ownership in any way.              *
c    The same holds for the way the program is structured to handle    *
c    window events, which can be found in the example program sources  *
c    which are part of the GLUT distribution.                          *
c                                                                      *
c***********************************************************************
c
c
      block data
      common /strdat/ ncurves, npoint, curve(15,20000)
      common /layout/ nscrnx, nscrny, nxpix, nypix, numfree,
     *                maxpar, maxpix, margin, nxorg, nyorg,
     *                icsm, icob, icross, ibackgr
      common /button/ ibutx(38), ibuty(38), ibutz(38), ibx, iby
      common /curves/ active, smoothed, reverse, sqrooted, overall,
     *                icolor, maxstep, nstep, step, nummax, nf, nl,
     *                valmin, valmax, scale, tickval, clipp, clipm,
     *                siglim, npcur, current(20000), cursm(20000),
     *                ncurrent(20000), curtmp(20000), itmp(20000)
      logical*4       active, smoothed, reverse, sqrooted, overall
      data ibutx /    0,   80,  160,  240,  320,  400,  480,  560,
     *              640,  720,  800,  880,    0,   80,  160,  240,
     *              320,  400,  480,  560,  640,  720,  800,  880,
     *              960,  960,  960,  960,  960,  960,  960,  960,
     *              960,  960,  960,  960,  960,  960 /
      data ibuty /  100,  100,  100,  100,  100,  100,  100,  100,
     *              100,  100,  100,  100,   60,   60,   60,   60,
     *               60,   60,   60,   60,   60,   60,   60,   60,
     *              580,  540,  500,  460,  420,  380,  340,  300,
     *              260,  220,  180,  140,  100,   60 /
      data ibutz           / 38 * 1 /
      data ibx, iby        / 0, 0 /
      data nscrnx, nscrny  / 1280, 1024 /
      data nxpix, nypix    / 1040,  620 /
c     data margin, numfree /   44,  140 /
      data margin, numfree /   64,  140 /
c     data nxorg, nyorg    /   46,  140 /
      data nxorg, nyorg    /   66,  140 /
      data maxpix, maxpar  /  480,  479 /
      data maxstep, nstep  /   20,   99 /
c     data tickval         /  2.0 /
      data tickval         /  0.02 /
      data icsm, icob, icross, ibackgr / 1, 2, 3, 4 /
      data active          / .false. /
      data smoothed        / .false. /
      data reverse         / .false. /
      data sqrooted        / .false. /
      data overall         / .true.  /
      data siglim          / 3.0 /
      data clipp, clipm    / 100, -100 /
      end
c
c
      program systerplot
c
      include         "systerplot.inc"
      character*64    argument, ccode, filename
      integer*4       argm_length, ccode_length
      character*96    title
      external        display
      external        reshape
      external        mouse
c
      call callcode ( argument, argm_length, ccode, ccode_length )
      filename = argument(1:argm_length) // '.syster'
      call oldfile ( 1, filename )
c
      read (1,10,end=60) ncurves
   10 format (i6)
      if ( ncurves .gt. 15 ) then
         write (6,20) ncurves
   20    format ('warning: number of curves (ncurves) .gt. 15: ',i5,/,
     *           'the program is dimensioned for 15 curves only',/,
     *           'only the first 15 curves will be used...')
         ncurves = 15
      end if
      npoint = 0
   30 npoint = npoint + 1
      read (1,40,end=50) ( curve(j,npoint), j = 1, ncurves )
   40 format (15f10.0)
      goto 30
c
   50 npoint = npoint - 1
      iwx = ( nscrnx - nxpix ) / 2
      iwy = ( nscrny - nypix ) / 2
c     call glutInitDisplayMode( GLUT_DOUBLE + GLUT_RGB )
      call glutInitDisplayMode( GLUT_RGB )
      call glutInitWindowSize( nxpix-1, nypix-1 )
      call glutInitWindowPosition( iwx, iwy )
      call glutInit
c now that the overall size has been set, reset nxpix, nypix to the
c values for the 'drawing board'
      nxpix = 960
      nypix = 620
      title = 'Systematic error curve for'
      title = title(1:leng( title )) // ' ' // ccode(1:ccode_length)
      iw = glutCreateWindow( title )
      call glutSetIconTitle( 'SystErrPlot' )
      call fglClearColor( 0.0, 0.0, 0.25, 1.0 )
      call glutReshapeFunc( reshape )
      call glutDisplayFunc( display )
      call glutMouseFunc( mouse )
      call glutMainLoop
c
   60 stop
      end
c
c
      subroutine display
      include "GL/fgl.h"
      include "GL/fglut.h"
c
      call fglClear( GL_COLOR_BUFFER_BIT )
      call gfxdrawwindow
      return
      end
c
c
      subroutine reshape ( iwidth, iheight )
      include "GL/fgl.h"
      include "GL/fglut.h"
c
      width = float( iwidth )
      height = float( iheight )
      call fglViewport ( 0, 0, iwidth, iheight )
      call fglMatrixMode( GL_PROJECTION )
      call fglLoadIdentity ()
      call fgluOrtho2D( dble(0.0), dble(width),
     *                  dble(0.0), dble(height) )
      return
      end
c
c
      subroutine gfxdrawwindow
      include         "systerplot.inc"
c
      call draw_coloured_rectangle ( 0,   0, 1039,  59, 0 )
      call draw_coloured_rectangle ( 0, 140,  959, 619, ibackgr )
      call makebutton ( 0 )
      do 10 i = 2, 38
         call fglRasterPos2i( ibutx(i), ibuty(i) )
         call fglCopyPixels( 0, 100, 80, 40, GL_COLOR )
   10 continue
      call chardraw ( 0, 0, 'Fo', 1 )
      call chardraw ( 0, 0, 'Fc', 2 )
      call chardraw ( 0, 0, 'dF / Fo', 3 )
      call chardraw ( 0, 0, 'dF / s(Fo)', 4 )
      call chardraw ( 0, 0, 'Fo / Fc', 5 )
      call chardraw ( 0, 0, 'Fo / s(Fo)', 6 )
      call chardraw ( 0, 0, 'abscor', 7 )
      call chardraw ( 0, 0, '    drift', 8 )
      call chardraw ( 0, 0, 'weight', 9 )
      call chardraw ( 0, 0, 'BckgrRat', 10 )
      call chardraw ( 0, 0, 'dF / <Fo>', 11 )
      call chardraw ( 0, 0, 'dF/<s(Fo)>', 12 )
      call chardraw ( 0, 0, 'vs. Fo', 13 )
      call chardraw ( 0, 0, 'vs. Fc', 14 )
      call chardraw ( 0, 0, 'vs. theta', 15 )
      call chardraw ( 0, 0, 'vs. phi', 16 )
      call chardraw ( 0, 0, 'vs. omega', 17 )
      call chardraw ( 0, 0, 'vs. kappa', 18 )
      call chardraw ( 0, 0, 'vs. abscor', 19 )
      call chardraw ( 0, 0, '     vs. drift', 20 )
      call chardraw ( 0, 0, 'vs. weight', 21 )
      call chardraw ( 0, 0, 'vs. XRtime', 22 )
      call chardraw ( 0, 0, '<dF> / Fo', 23 )
      call chardraw ( 0, 0, '<dF>/s(Fo)', 24 )
      call chardraw ( 0, 0, '  >   >>  >>>', 25 )
      call chardraw ( 0, 0, '    ||<--  -->||', 26 )
      call chardraw ( 0, 0, '  <   <<  <<<', 27 )
      call chardraw ( 0, 0, ' << step >>', 28 )
      call chardraw ( 0, 0, 'sc frame', 29 )
      call chardraw ( 0, 0, '  << +clip >>', 30 )
      call chardraw ( 0, 0, '  << -clip >>', 31 )
      call chardraw ( 0, 0, '  reverse x', 32 )
      call chardraw ( 0, 0, '   sqrt(y)', 33 )
      call chardraw ( 0, 0, 'smooth', 34 )
      call chardraw ( 0, 0, '    << slim >>', 35 )
      call chardraw ( 0, 0, '<< <Fo> >>', 36 )
      call chardraw ( 0, 0, ' reset view', 37 )
      call chardraw ( 0, 0, 'EXIT', 38 )
c restore if window is redrawn
      do 20 i = 1, 38
         if ( ibutz( i ) .eq. -1 ) call makebutton( -i )
   20 continue
      if ( active ) then
         call draw_axes
         call gfxdrawboard
      end if
      return
      end
c
c
      subroutine mouse ( ibut, istate, ix, iy )
      include         "systerplot.inc"
      external timer
c
      if ( istate .ne. GLUT_DOWN ) goto 900
      if ( ibut .ne. GLUT_LEFT_BUTTON ) goto 900
      iy = 618 - iy
      if ( iy .lt. 60 ) goto 900
      if ( iy .ge. 140 .and. ix .le. 959 ) goto 900
      ibut = 0
      if ( ix .le. 959 ) then
         ibut = ( ix / 80 ) + 1
         if ( iy .le. 99 ) ibut = ibut + 12
      else
         ibut = 38 - ( ( iy - 60 ) / 40 )
      end if
      if ( ibut .ne. 29 .and. ibut .ne. 32 .and. ibut .ne. 33 .and.
     *     ibut .ne. 34 ) then
         if ( ibutz( ibut ) .eq. -1 ) return
      end if
      if ( ( ibut .ge.  1 .and. ibut .le. 12 ) .or.
     *     ( ibut .ge. 23 .and. ibut .le. 24 ) ) then
         do 1 i = 1, 12
            if ( ibutz( i ) .eq. -1 ) call makebutton( i )
    1    continue
         if ( ibutz( 23 ) .eq. -1 ) call makebutton( 23 )
         if ( ibutz( 24 ) .eq. -1 ) call makebutton( 24 )
         call draw_coloured_rectangle ( 0, 140,  959, 619, ibackgr )
         if ( ibutz( 29 ) .eq. -1 ) call makebutton ( 29 )
         active = .false.
      end if
      if ( ibut .ge. 13 .and. ibut .le. 22 ) then
         do 2 i = 13, 22
            if ( ibutz( i ) .eq. -1 ) call makebutton( i )
    2    continue
         call draw_coloured_rectangle ( 0, 140,  959, 619, ibackgr )
         if ( ibutz( 29 ) .eq. -1 ) call makebutton ( 29 )
         active = .false.
      end if
      if ( ibut .ge.  1 .and. ibut .le. 24 ) then
         call makebutton ( -ibut )
         if ( ( ibut .eq.  1 .and. ibutz( 13 ) .eq. -1 ) .or.
     *        ( ibut .eq.  2 .and. ibutz( 14 ) .eq. -1 ) .or.
     *        ( ibut .eq.  7 .and. ibutz( 19 ) .eq. -1 ) .or.
     *        ( ibut .eq.  8 .and. ibutz( 20 ) .eq. -1 ) .or.
     *        ( ibut .eq.  9 .and. ibutz( 21 ) .eq. -1 ) ) then
            call makebutton ( ibut + 12 )
            ibx = 0
         end if
         if ( ( ibut .eq. 13 .and. ibutz(  1 ) .eq. -1 ) .or.
     *        ( ibut .eq. 14 .and. ibutz(  2 ) .eq. -1 ) .or.
     *        ( ibut .eq. 19 .and. ibutz(  7 ) .eq. -1 ) .or.
     *        ( ibut .eq. 20 .and. ibutz(  8 ) .eq. -1 ) .or.
     *        ( ibut .eq. 21 .and. ibutz(  9 ) .eq. -1 ) ) then
            call makebutton( ibut - 12 )
            iby = 0
         end if
         if ( ibut .le. 12 .or. ibut .eq. 23 .or. ibut .eq. 24 ) then
            iby = ibut
          else
            ibx = ibut
          end if
          if ( ibx .ne. 0 .and. iby .ne. 0 ) then
             call fill_current
          end if
      else if ( ibut .eq.  25 ) then
c shift viewframe gradually to the right, i.e. to higher x-values
c over one third (ix .le. 985), one whole (ix .ge. 1013) or one half
c (otherwise) of the viewframe
         if ( active ) then
            if ( nl .lt. npcur ) then
               if ( ix .le. 985 ) then
                  nl = nl - 1 + nummax / 3
               else if ( ix .ge. 1013 ) then
                  nl = nl - 1 + nummax
               else
                  nl = nl - 1 + nummax / 2
               end if
               if ( nl .gt. npcur ) nl = npcur
               nf = nl - nummax + 1
               if ( nf .lt. 1 ) then
                  nf = 1
                  nl = nummax
                  if ( nl .gt. npcur ) nl = npcur
               end if
               call draw_axes
               call gfxdrawboard
            end if
         end if
      else if ( ibut .eq. 26 ) then
c shift viewframe to the extreme left, i.e. to the lowest x-values
c (ix .le. 998) or to the extreme right, i.e. to the highest x-values
c (otherwise)
         if ( active ) then
            if ( ix .le. 998 ) then
               if ( nf .gt. 1 ) then
                  nf = 1
                  nl = nummax
                  if ( nl .gt. npcur ) nl = npcur
               end if
            else
               if ( nl .lt. npcur ) then
                  nl = npcur
                  nf = nl - nummax + 1
                  if ( nf .lt. 1 ) nf = 1
               end if
            end if
            call draw_axes
            call gfxdrawboard
         end if
      else if ( ibut .eq. 27 ) then
c shift viewframe gradually to the left, i.e. to lower x-values
c over one third (ix .le. 985), one whole (ix .ge. 1013) or one half
c (otherwise) of the viewframe
         if ( active ) then
            if ( nf .gt. 1 ) then
               if ( ix .le. 985 ) then
                  nf = nf + 1 - nummax / 3
               else if ( ix .ge. 1013 ) then
                  nf = nf + 1 - nummax
               else
                  nf = nf + 1 - nummax / 2
               end if
               if ( nf .lt. 1 ) nf = 1
               nl = nf - 1 + nummax
               if ( nl .gt. npcur ) then
                  nl = npcur
                  nf = nl + 1 - nummax
                  if ( nf .lt. 1 ) nf = 1
               end if
               call draw_axes
               call gfxdrawboard
            end if
         end if
      else if ( ibut .eq. 28 ) then
c decrease step size per data point (ix .le. 998) or increase step size
c per data point (otherwise)
         if ( active ) then
            nstadd = int ( 0.4 * float( nstep ) ) + 1
            if ( ix .le. 998 ) then
               if ( nstep .ge. 1 ) then
                  nstep = nstep - nstadd
c                 if ( nstep .lt. 1 ) nstep = 1
                  if ( nstep .lt. 0 ) nstep = 0
               end if
            else
               if ( nstep .lt. maxstep ) then
                  nstep = nstep + nstadd
                  if ( nstep .gt. maxstep ) nstep = maxstep
               end if
            end if
            if ( nstep .ge. 1 ) then
                nummax = ( nxpix - margin - 2 ) / nstep
            else
                nummax = npcur
                step = float( nxpix - margin - 2 ) / float( npcur )
            end if
            nf = 1
            nl = nummax
            call draw_axes
            call gfxdrawboard
         end if
      else if ( ibut .eq. 29 ) then
c on/off switch: scale plot to overall extreme values (off, default) or
c scale plot to extreme values within current viewframe (on)
         if ( active ) then
            if ( ibutz( 29 ) .eq. -1 ) then
               call makebutton( 29 )
               overall = .true.
            else
               call makebutton( -29 )
               overall = .false.
            end if
            call draw_axes
            call gfxdrawboard
         end if
      else if ( ibut .eq. 30 ) then
c decrease clipp (ix .le. 998) or increase clipp (otherwise)
c datapoints with values larger than clipp are skipped
         if ( active ) then
            clpad = ( 0.4 * clipp ) + 1.0
            if ( ix .le. 998 ) then
               clipp = clipp - clpad
            else
               clipp = clipp + clpad
            end if
            clipp = float( nint( clipp ) )
            if ( clipp .lt. 1.0 ) clipp = 1.0
            call fill_current
         end if
      else if ( ibut .eq. 31 ) then
c decrease clipm (ix .le. 998) or increase clipm (otherwise)
c datapoints with values smaller than clipm are skipped
         if ( active ) then
            clpad = ( -0.4 * clipm ) + 1.0
            if ( ix .le. 998 ) then
               clipm = clipm - clpad
            else
               clipm = clipm + clpad
            end if
            clipm = float( nint( clipm ) )
            if ( clipm .gt. -1.0 ) clipm = -1.0
            call fill_current
         end if
      else if ( ibut .eq. 32 ) then
c on/off switch: plot from lowest to highest value (off, default), or
c plot from highest to lowest value (on)
         if ( active ) then
            if ( ibutz( 32 ) .eq. -1 ) then
               call makebutton( 32 )
               reverse = .false.
            else
               call makebutton( -32 )
               reverse = .true.
            end if
            call fill_current
         end if
      else if ( ibut .eq. 33 ) then
c on/off switch: plot the normal y-values (off, default), or plot
c the square root of the y-values (on)
         if ( active ) then
            if ( ibutz( 33 ) .eq. -1 ) then
               call makebutton( 33 )
               sqrooted = .false.
            else
               call makebutton( -33 )
               sqrooted = .true.
            end if
            call fill_current
         end if
      else if ( ibut .eq. 34 ) then
c on/off switch: do not add smoothed curve (off, default), or add
c 5-point smoothed curve (on) to plot
         if ( active ) then
            if ( ibutz( 34 ) .eq. -1 ) then
               call makebutton( 34 )
               smoothed = .false.
            else
               call makebutton( -34 )
               smoothed = .true.
            end if
            call fill_current
         end if
      else if ( ibut .eq. 35 ) then
c decrease siglim (ix .le. 998) or increase siglim (otherwise)
c refls with abs(Fo/sig(Fo)) below siglim are skipped
         if ( active ) then
            if ( ix .le. 998 ) then
               siglim = siglim - 1.0
            else
               siglim = siglim + 1.0
            end if
            if ( siglim .lt. 0.0 ) siglim = 0.0
            call fill_current
         end if
      else if ( ibut .eq. 36 ) then
c decrease the nuber of <Fobs> intervals (ix .le. 998) or increase
c the number of <Fobs> intervals (otherwise)
c this number is used to calculate mean values for Fobs and/or
c sigma(Fobs) over evenly divided intervals of Fobs
         if ( active ) then
            if ( ix .le. 998 ) then
            else
            end if
         end if
      else if ( ibut .eq. 37 ) then
c reset view to default values
         if ( active ) then
            reverse = .false.
            sqrooted = .false.
            smoothed = .false.
            overall  = .true.
            call makebutton( 29 )
            call makebutton( 32 )
            call makebutton( 33 )
            siglim = 3.0
            nstep = 99
            clipp =  100.0
            clipm = -100.0
            call fill_current
         end if
      else if ( ibut .eq. 38 ) then
c stop the program
         stop
      else
         write (6,400) ix, iy, ibut
  400    format ('illegal button detected: ',i4,':',i4,', button: ',i2)
      end if
  900 call gluttimerfunc( 1000, timer, 0 )
      return
      end
c
c
      subroutine fill_current
      include         "systerplot.inc"
      dimension ixpos( 10 ), iypos( 12 )
      data ixpos /  2,  3,  6,  7,  8,  9, 11, 12, 13,  5 /
      data iypos /  2,  3,  0,  0,  0,  0, 11, 12, 13, 10,  0,  0 /
c
c     curve( 1,..) = reflection number
c     curve( 2,..) = observed structure factor Fobs
c     curve( 3,..) = calculated structure factor, scaled to Fobs
c     curve( 4,..) = sigma( Fobs )
c     curve( 5,..) = xraytime
c     curve( 6,..) = theta
c     curve( 7,..) = phik
c     curve( 8,..) = omk
c     curve( 9,..) = kappa
c     curve(10,..) = background ratio
c     curve(11,..) = abscor
c     curve(12,..) = drift
c     curve(13,..) = weight
c
      ixval = ixpos( ibx - 12 )
      call sortf ( ixval )
      iyval = 0
      if ( iby .le. 12 ) iyval = iypos ( iby )
      if ( iyval .ne. 0 ) then
         npcur = 0
         do 10 i = 1, npoint
            if ( abs(curve(2,itmp(i))/curve(4,itmp(i))).ge. siglim) then
               npcur = npcur + 1
               curtmp(npcur) = curve(iyval,itmp(i))
            end if
   10    continue
      end if
      if ( iyval .eq. 0 ) then
         if      ( iby .eq. 3 ) then
c           ( Fobs - Fcalc ) / Fobs : (2-3)/2
            npcur = 0
            do 20 i = 1, npoint
               if ( abs(curve(2,itmp(i))/curve(4,itmp(i))).ge. siglim)
     *                                                             then
                  npcur = npcur + 1
                  curtmp(npcur) =
     *            (curve(2,itmp(i))-curve(3,itmp(i))) / curve(2,itmp(i))
               end if
   20       continue
         else if ( iby .eq. 4 ) then
c           ( Fobs - Fcalc ) / sigma( Fobs ) : (2-3)/4
            npcur = 0
            do 30 i = 1, npoint
               if ( abs(curve(2,itmp(i))/curve(4,itmp(i))) .ge. siglim)
     *                                                             then
                  npcur = npcur + 1
                  curtmp(npcur) =
     *            (curve(2,itmp(i))-curve(3,itmp(i))) / curve(4,itmp(i))
               end if
   30       continue
         else if ( iby .eq. 5 ) then
c           Fobs / Fcalc : 2/3
            npcur = 0
            do 40 i = 1, npoint
               if ( abs(curve(2,itmp(i))/curve(4,itmp(i))) .ge. siglim)
     *                                                             then
                  npcur = npcur + 1
                  curtmp(npcur) = curve(2,itmp(i)) / curve(3,itmp(i))
               end if
   40       continue
         else if ( iby .eq. 6 ) then
c           Fobs / sigma( Fobs ) : 2/4
            npcur = 0
            do 50 i = 1, npoint
               if ( abs(curve(2,itmp(i))/curve(4,itmp(i))) .ge. siglim)
     *                                                             then
                  npcur = npcur + 1
                  curtmp(npcur) = curve(2,itmp(i)) / curve(4,itmp(i))
               end if
   50       continue
         else if ( iby .eq. 11 ) then
            npcur = 0
         else if ( iby .eq. 12 ) then
            npcur = 0
         else if ( iby .eq. 23 ) then
            npcur = 0
         else if ( iby .eq. 24 ) then
            npcur = 0
         end if
      end if
      if ( npcur .eq. 0 ) return
      do 60 i = 1, npcur
         if ( curtmp(i) .gt. clipp ) curtmp(i) = clipp
         if ( curtmp(i) .lt. clipm ) curtmp(i) = clipm
         if ( reverse ) then
             current(i) = curtmp(i)
         else
             current(npcur+1-i) = curtmp(i)
         end if
   60 continue
      if ( sqrooted ) then
         do 70 i = 1, npcur
            if ( current(i) .ge. 0.0 ) then
               current(i) = sqrt( current(i) )
            else
               current(i) = - sqrt( - current(i) )
            end if
   70    continue
      end if
      if ( smoothed ) then
         do 80 i = 3, npcur - 2
            cursm(i) = ( current(i-2)       +
     *                   current(i-1) * 2.0 +
     *                   current(i  ) * 4.0 +
     *                   current(i+1) * 2.0 +
     *                   current(i+2)       ) / 10.0
   80    continue
         cursm(1) = ( current(1) * 2.0 +
     *                current(2)       ) / 3.0
         cursm(2) = ( current(1)       +
     *                current(2) * 2.0 +
     *                current(3)       ) / 4.0
         cursm(npcur-1) = ( current(npcur-2)       +
     *                      current(npcur-1) * 2.0 +
     *                      current(npcur  )       ) / 4.0
         cursm(npcur)   = ( current(npcur-1) +
     *                      current(npcur  ) * 2.0 ) / 3.0
      end if
      if ( nstep .eq. 99 ) then
         nstep = float( nxpix - margin - 2 ) / float( npcur )
         if ( nstep .gt. maxstep ) nstep = maxstep
         if ( nstep .lt.       2 ) nstep = 2
         nummax = ( nxpix - margin - 2 ) / nstep
         nf = 1
         nl = nummax
      end if
      call draw_axes
      call gfxdrawboard
      return
      end
c
c
      subroutine timer( ival )
      return
      end
c
c
      subroutine draw_axes
      include         "systerplot.inc"
c     character*4 tick
      character*6 tick
c
      valmax = -10000000.0
      valmin =  10000000.0
      if ( overall ) then
         ib = 1
         ie = npcur
      else
         ib = nf
         ie = nl
      end if
      do 5 i = ib, ie
         if ( current(i) .gt. valmax ) valmax = current(i)
         if ( current(i) .lt. valmin ) valmin = current(i)
    5 continue
      order = alog10( valmax - valmin )
c     if ( order .ge. 0.0 ) then
         range = 10.0 ** float( int( order ) )
         imult = int( ( valmax - valmin ) / range ) + 1
         im = -1
         if ( imult .gt. 2 ) then
            im = imult / 2
            if ( im * 2 .ne. imult ) imult = imult + 1
         end if
         tickval = ( float( imult ) * range ) / 10.0
         if      ( ( valmax - valmin ) / tickval .le. 1.0 ) then
            tickval = tickval / 10.0
         else if ( ( valmax - valmin ) / tickval .lt. 2.0 ) then
            tickval = tickval / 5.0
         else if ( ( valmax - valmin ) / tickval .lt. 2.5 ) then
            tickval = tickval / 4.0
         else if ( ( valmax - valmin ) / tickval .lt. 5.0 ) then
            tickval = tickval / 2.0
         end if
c     end if
      p = float( int( abs( valmax ) / tickval ) )
      if ( valmax .gt. 0.1 ) then
         valmax = ( p + 1.5 ) * tickval
      else if ( valmax .lt. -0.1 ) then
         valmax = -1.0 * ( p - 1.5 ) * tickval
         if ( valmax .lt. 0.0 ) valmax = 0.0
      else
         valmax = tickval * 0.5
      end if
      p = float( int( abs( valmin ) / tickval ) )
      if ( valmin .gt. 0.1 ) then
         valmin = ( p - 0.5 ) * tickval
         if ( valmin .lt. 0.0 ) valmin = tickval * -0.5
      else if ( valmin .lt. -0.1 ) then
         valmin = -1.0 * ( p + 1.5 ) * tickval
      else
         valmin = tickval * -0.5
      end if
      call writeval ( tick, valmax, -1 )
      call writeval ( tick, valmin, -2 )
      scale  = float( maxpix ) / ( valmax - valmin )
      call draw_coloured_rectangle ( 0, 140, 959, 619, ibackgr )
      ix1 = margin
      iy1 = numfree
      ix2 = ix1 + 1
      iy2 = nypix - 1
      call draw_coloured_rectangle ( ix1, iy1, ix2, iy2, icross )
      nwan = nint( (0.0 - valmin) * scale )
      ix1 = margin + 2
      iy1 = numfree + nwan
      ix2 = nxpix - 1
      iy2 = iy1 + 1
      if ( iy1 .ge. numfree .and. iy2 .lt. nypix ) then
         call draw_coloured_rectangle ( ix1, iy1, ix2, iy2, icross )
         val = 0.0
         call writeval ( tick, val, 0 )
         ix1 = 5
         iy1 = iy1 - 3
         call chardraw ( ix1, iy1, tick, 0 )
      end if
      val = 0.0
   20 val = val + tickval
      if ( val .ge. valmax ) goto 30
      ix1 = margin + 2
      iy1 = numfree + nint( (val - valmin) * scale )
      ix2 = nxpix - 1
      iy2 = iy1
      if ( iy1 .ge. numfree .and. iy2 .lt. nypix ) then
         call draw_coloured_rectangle ( ix1, iy1, ix2, iy2, icross )
         call writeval ( tick, val, 0 )
         ix1 = 5
         iy1 = iy1 - 3
         call chardraw ( ix1, iy1, tick, 0 )
      end if
      goto 20
   30 val = 0.0
   40 val = val - tickval
      if ( val .le. valmin ) goto 60
      ix1 = margin + 2
      iy1 = numfree + nint( (val - valmin) * scale )
      ix2 = nxpix - 1
      iy2 = iy1
      if ( iy1 .ge. numfree .and. iy2 .lt. nypix ) then
c        if ( val .gt. 0.99 .and. val .lt. 1.01 ) goto 50
         call draw_coloured_rectangle ( ix1, iy1, ix2, iy2, icross )
   50    call writeval ( tick, val, 0 )
         ix1 = 5
         iy1 = iy1 - 3
         call chardraw ( ix1, iy1, tick, 0 )
      end if
      goto 40
   60 return
      end
c
c
      subroutine gfxdrawboard
      include         "systerplot.inc"
c
      do 10 i = nf, nl
         ncurrent(i) = nint( ( current(i) - valmin ) * scale )
         if ( ncurrent(i) .lt.      0 ) ncurrent(i) = 0
         if ( ncurrent(i) .gt. maxpar ) ncurrent(i) = maxpar
   10 continue
      call fglLineWidth( 2.0 )
      icolor = icob
      call draw_curve
c
      if ( smoothed ) then
         do 20 i = nf, nl
            ncurrent(i) = nint( ( cursm(i) - valmin ) * scale )
            if ( ncurrent(i) .lt.      0 ) ncurrent(i) = 0
            if ( ncurrent(i) .gt. maxpar ) ncurrent(i) = maxpar
   20    continue
c        call fglLineWidth( 3.0 )
         icolor = icsm
         call draw_curve
      end if
c
      active = .true.
      return
      end
c
c
      subroutine draw_coloured_rectangle ( ixb, iyb, ixe, iye, ifill )
      include         "systerplot.inc"
c
      call kleur ( ifill )
      call fglRecti ( ixb, iyb, ixe+1, iye+1 )
c
      return
      end
C
      subroutine draw_curve
      include         "systerplot.inc"
c
      call kleur ( icolor )
      call fglBegin( GL_LINE_STRIP )
         ix = nxorg
         iy = nyorg + ncurrent( nf )
         call fglVertex2i( ix, iy )
         ip = margin + 2
         fip = float( ip )
         do 10 i = nf+1, nl
            if ( nstep .ge. 1 ) then
               ip = ip + nstep
            else
               fip = fip + step
               ip = int( fip )
            end if
            ix = ip
            iy = nyorg + ncurrent(i)
            call fglVertex2i( ix, iy )
   10    continue
      call fglEnd
      return
      end
c
c
      subroutine kleur ( ikleur )
      include "GL/fgl.h"
      include "GL/fglut.h"
      dimension   red(0:7), green(0:7), blue(0:7)
      data red   / 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 /
      data green / 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0 /
      data blue  / 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /
      data falpha / 1.0 /
c
      if ( ikleur .ge. 0 .and. ikleur .le. 7 ) then
         fred   = red( ikleur )
         fgreen = green( ikleur )
         fblue  = blue( ikleur )
      else
         fred   = 1.0
         fgreen = 0.0
         fblue  = 0.0
         write (6,10) ikleur
   10    format ('undefined color, ikleur = ',i5)
      end if
      call fglColor3f( fred, fgreen, fblue )
      return
      end
c
c
      subroutine chardraw ( ix, iy, string, ibut )
      include         "systerplot.inc"
      character*(*)  string
      real*4         curcol(4)
      character*1    ch
c  Stroke and font names must be explicitly declared as
c  external instead of relying on "GL/fglut.h" because
c  the IRIX Fortran compiler does not know how to only
c  link in used external data symbols
c     integer*4 GLUT_STROKE_ROMAN
c     integer*4 GLUT_STROKE_MONO_ROMAN
c     integer*4 GLUT_BITMAP_9_BY_15
c     integer*4 GLUT_BITMAP_8_BY_13
c     integer*4 GLUT_BITMAP_TIMES_ROMAN_10
c     integer*4 GLUT_BITMAP_TIMES_ROMAN_24
c     integer*4 GLUT_BITMAP_HELVETICA_10
c     integer*4 GLUT_BITMAP_HELVETICA_12
c     integer*4 GLUT_BITMAP_HELVETICA_18
c
c     external GLUT_STROKE_ROMAN
c     external GLUT_STROKE_MONO_ROMAN
      external GLUT_BITMAP_9_BY_15
c     external GLUT_BITMAP_8_BY_13
c     external GLUT_BITMAP_TIMES_ROMAN_10
c     external GLUT_BITMAP_TIMES_ROMAN_24
c     external GLUT_BITMAP_HELVETICA_10
      external GLUT_BITMAP_HELVETICA_12
c     external GLUT_BITMAP_HELVETICA_18
c
      nc = leng( string )
      if ( nc .eq. 0 ) return
      if ( ix .eq. 0 .and. iy .eq. 0 .and. ibut .ne. 0 ) then
         call fglGetFloatv( GL_CURRENT_COLOR, curcol )
         call kleur( 0 )
c
c     function glutBitmapWidth doesn't work properly, it seems.
c     use rather unsatisfactory constant character width instead.
c        iw = 0
c        do 10 i = 1, nc
c           ch = string(i:i)
c           ich = ichar( ch )
c           iw = iw + glutBitmapWidth( GLUT_BITMAP_HELVETICA_12,
c    *                                 ichar(ch) )
c  10    continue
         iw = nc * 6
         call fglRasterPos2i( ibutx( ibut ) + 40 - ( iw / 2 ),
     *                        ibuty( ibut ) + 17 )
         do 20 i = 1, nc
            ch = string(i:i)
            call glutBitmapCharacter( GLUT_BITMAP_HELVETICA_12,
     *                                ichar(ch) )
   20    continue
         call fglColor3f( curcol(1), curcol(2), curcol(3) )
      else
         if ( iy .gt. nypix   - 11 .or.
     *        iy .lt. numfree +  2 ) return
         call fglRasterPos2i( ix, iy )
         do 30 i = 1, nc
            ch = string(i:i)
            call glutBitmapCharacter( GLUT_BITMAP_9_BY_15, ichar(ch) )
   30    continue
      end if
      return
      end
c
c
      subroutine makebutton ( kind )
      include         "systerplot.inc"
      byte      ival1a(240), ival1b(240), ival1c(240), ival1d(240)
      byte      ival1e(240), ival1f(240), ival1g(240), ival1h(240)
      integer*4 ival4a(60),  ival4b(60),  ival4c(60),  ival4d(60)
      integer*4 ival4e(60),  ival4f(60),  ival4g(60),  ival4h(60)
      equivalence ( ival1a(1), ival4a(1) )
      equivalence ( ival1b(1), ival4b(1) )
      equivalence ( ival1c(1), ival4c(1) )
      equivalence ( ival1d(1), ival4d(1) )
      equivalence ( ival1e(1), ival4e(1) )
      equivalence ( ival1f(1), ival4f(1) )
      equivalence ( ival1g(1), ival4g(1) )
      equivalence ( ival1h(1), ival4h(1) )
      byte      ival1i(240), ival1j(240), ival1k(240), ival1l(240)
      byte      ival1m(240), ival1n(240), ival1o(240)
      integer*4 ival4i(60),  ival4j(60),  ival4k(60),  ival4l(60)
      integer*4 ival4m(60),  ival4n(60),  ival4o(60)
      equivalence ( ival1i(1), ival4i(1) )
      equivalence ( ival1j(1), ival4j(1) )
      equivalence ( ival1k(1), ival4k(1) )
      equivalence ( ival1l(1), ival4l(1) )
      equivalence ( ival1m(1), ival4m(1) )
      equivalence ( ival1n(1), ival4n(1) )
      equivalence ( ival1o(1), ival4o(1) )
      equivalence ( ival1h(1), ival4h(1) )
      data ival4a / 60 * 690563369 /
      data ival4b / 690563550, -555853479, 56 * 1499027801,
     *              1499015465, 690563369 /
      data ival4c / 690563550, -555819298, -564569767,
     *              54 * 1499027801, 1499027753, 2 * 690563369 /
      data ival4d / 690563550, 2 * -555819298, 54 * 1499027801,
     *              3 * 690563369 /
      data ival4e / 690563550, 2 * -555819298, 54 * -1583242847,
     *              3 * 690563369 /
      data ival4f / 690563550, -555819298, -557003572,
     *              54 * -858993460, -858993623, 2 * 690563369 /
      data ival4g / 690563550, -555823924, 56 * -858993460,
     *              -859035351, 690563369 /
      data ival4h / 690563532, 58 * -858993460, -869717719 /
      data ival4i / 60 * -555819298 /
      data ival4j / 690563369, 690609886, 58 * -555819298 /
      data ival4k / 2 * 690563369, 702471902, 57 * -555819298 /
      data ival4l / 2 * 690563369, 698458291, 54 * -1583242847,
     *              -1583242786, 2 * -555819298 /
      data ival4m / 690563369, 690575705, 56 * 1499027801,
     *              1499061982, -555819298 /
      data ival4n / 690563417, 58 * 1499027801, 1507778270 /
      data ival4o / 60 * 1499027801 /
c
      if ( kind .eq. 0 ) then
         call fglRasterPos2i( 0, 100 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1a )
         call fglRasterPos2i( 0, 101 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1b )
         call fglRasterPos2i( 0, 102 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1c )
         call fglRasterPos2i( 0, 103 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1d )
         do 10 i = 1, 32
            call fglRasterPos2i( 0, 103 + i )
            call fglDrawPixels( 80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1e)
   10    continue
         call fglRasterPos2i( 0, 136 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1f )
         call fglRasterPos2i( 0, 137 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1g )
         call fglRasterPos2i( 0, 138 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1h )
         call fglRasterPos2i( 0, 139 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1a )
      else if ( kind .gt. 0 .and. kind .le.  38 ) then
         ix = ibutx( kind )
         iy = ibuty( kind )
         call fglRasterPos2i( ix, iy )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1a )
         call fglRasterPos2i( ix, iy +  1 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1b )
         call fglRasterPos2i( ix, iy +  2 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1c )
         call fglRasterPos2i( ix, iy +  3 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1d )
         do 20 i = 4, 35
            call fglRasterPos2i( ix, iy + i )
            call fglDrawPixels( 3, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1e )
            call fglRasterPos2i( ix + 78, iy + i )
            call fglDrawPixels( 3, 1, GL_RGB, GL_UNSIGNED_BYTE,
     *                          ival1e(232) )
   20    continue
         call fglRasterPos2i( ix, iy + 36 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1f )
         call fglRasterPos2i( ix, iy + 37 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1g )
         call fglRasterPos2i( ix, iy + 38 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1h )
         call fglRasterPos2i( ix, iy + 39 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1a )
         ibutz( kind ) = 1
      else if ( kind .lt. 0 .and. kind .ge. -38 ) then
         ix = ibutx( -kind )
         iy = ibuty( -kind )
         call fglRasterPos2i( ix, iy )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1i )
         call fglRasterPos2i( ix, iy +  1 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1j )
         call fglRasterPos2i( ix, iy +  2 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1k )
         call fglRasterPos2i( ix, iy +  3 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1l )
         do 30 i = 4, 35
            call fglRasterPos2i( ix, iy + i )
            call fglDrawPixels( 3, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1l )
            call fglRasterPos2i( ix + 78, iy + i )
            call fglDrawPixels( 3, 1, GL_RGB, GL_UNSIGNED_BYTE,
     *                          ival1l(232) )
   30    continue
         call fglRasterPos2i( ix, iy + 36 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1l )
         call fglRasterPos2i( ix, iy + 37 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1m )
         call fglRasterPos2i( ix, iy + 38 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1n )
         call fglRasterPos2i( ix, iy + 39 )
         call fglDrawPixels(  80, 1, GL_RGB, GL_UNSIGNED_BYTE, ival1o )
         ibutz( -kind ) = -1
      end if
      return
      end
c
c
      subroutine callcode ( argument, argm_length, ccode, ccode_length )
      character*(*)   argument, ccode
      integer*4       argm_length, ccode_length
      character*26    lowercase, uppercase
      data lowercase / 'abcdefghijklmnopqrstuvwxyz' /
      data uppercase / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /
c
c     subroutine callcode gets the first argument given in the program
c     calling sequence; if no argument is given it is asked for.
c     it is assumed that this argument is the compound code, which is
c     subsequently used to construct several file names.
c     the compound code itself is extracted and converted to uppercase
c     to be used for program use.
c
c     this subroutine uses two special statements which are / may be
c     platform dependent:
c     iargc()      integer function, returns the number of arguments
c                  given in the calling sequence, excluding the program
c                  name itself; returns 0 if no arguments are given.
c     getarg(i,s)  this subroutine gets the i-th argument from the
c                  calling sequence, excluding the program name, and
c                  stores the argument value in character string s.
c
      numarg = iargc()
      if ( numarg .eq. 0 ) then
         write (6,10)
   10    format ('Enter the compound code: ',$)
         read (5,20) argument
   20    format (a)
      else
         call getarg (1, argument)
      end if
      argm_length = leng( argument )
      ccode_length = 0
      do 30 i = 1, argm_length
         if ( argument(i:i) .ge. 'a' .and. argument(i:i) .le. 'z' .or.
     *        argument(i:i) .ge. 'A' .and. argument(i:i) .le. 'Z' .or.
     *        argument(i:i) .ge. '0' .and. argument(i:i) .le. '9' ) then
            ccode_length = ccode_length + 1
            goto 30
         else
            goto 40
         end if
   30 continue
   40 if ( ccode_length .lt. 1 ) then
         write (6,50) argument(1:argm_length)
   50    format ('Illegal compound code in calling argument ',a)
         stop
      else
         ccode = ' '
         ccode(1:ccode_length) = argument(1:ccode_length)
      end if
      do 70 i = 1, ccode_length
         do 60 j = 1, 26
            if ( ccode(i:i) .eq. lowercase(j:j) ) then
               ccode(i:i) = uppercase(j:j)
               goto 70
            end if
   60    continue
   70 continue
      return
      end
c
c
      integer function leng ( string )
      character*(*) string
c
c     integer function leng returns the length in bytes of the
c     meaningful part of a string, i.e. excluding trailing blanks;
c     leng = 0 if an empty string is supplied.
c     to be used if no such intrinsic is supplied by the compiler.
c     n.b. the standard intrinsic integer function len returns the
c     defined length of a string in bytes.
c
      ls = len( string )
      do 10 i = ls, 1, -1
         if ( string(i:i) .ne. ' ' ) goto 20
   10 continue
      leng = 0
      return
c
   20 leng = i
      return
      end
c
c
      subroutine oldfile ( iunit, filename )
      character*(*) filename
      logical*1 exists
c
c     subroutine oldfile tests the presence of file 'filename', assumed
c     to be present, and subsequently opens it on unit iunit.
c
      lenf = leng( filename )
      inquire ( file = filename(1:lenf), exist = exists,
     *          err = 20, iostat = ios )
      if ( exists ) then
         open ( file = filename(1:lenf), unit = iunit, status = 'old',
     *          err = 40, iostat = ios )
      return
      else
         write (6,10) filename(1:lenf)
   10    format ('Fatal error, file ',a,' does not exist')
      stop
      end if
c
   20 write (6,30) filename(1:lenf), ios
   30 format ('Fatal error, processor error during inquire on file ',
     * a,'; ios = ',i3)
      stop
c
   40 write (6,50) filename(1:lenf), ios
   50 format ('Fatal error, processor error during opening of file ',
     * a,'; ios = ',i3)
      stop
      end
c
c
c
c
      subroutine sortf ( n )
      include         "systerplot.inc"
c
c  sort curve(15,20000) on curve(n,..) creating pointers array itmp
c  on completion curve(n,itmp(1)) will be the largest value of curve(n,..)
c  only the first npoint elements will be sorted
c
      do 10 i = 1, npoint
   10 itmp(i) = i
      int = 2
   20 int = int + int
      if ( int .lt. npoint ) goto 20
      int = min0( npoint, (3*int)/4-1 )
   30 int = int/2
      ifin = npoint - int
      do 70 ii = 1, ifin
      i = ii
      j = i + int
      if ( curve(n,itmp(i)) .ge. curve(n,itmp(j)) ) goto 70
      k = itmp(j)
   40 itmp(j) = itmp(i)
      j = i
      i = i - int
      if ( i ) 60, 60, 50
   50 if ( curve(n,itmp(i) ) .lt. curve(n,k) ) goto 40
   60 itmp(j) = k
   70 continue
      if ( int .gt. 1 ) goto 30
      return
      end
c
c
      subroutine  writeval ( tick, v, key )
      character*(*) tick
      data km, kp / 0, 0 /
      if ( v .lt. -9999.0 ) v = -9999.0
      if ( v .gt. 99999.0 ) v = 99999.0
      if ( key .eq. -1 ) then
         if ( v .gt. -10000.0 .and. v .le.  -1000.0 ) kp = 0
         if ( v .gt.  -1000.0 .and. v .le.   -100.0 ) kp = 1
         if ( v .gt.   -100.0 .and. v .le.    -10.0 ) kp = 2
         if ( v .gt.    -10.0 .and. v .lt.      0.0 ) kp = 3
         if ( v .ge.      0.0 .and. v .lt.     10.0 ) kp = 4
         if ( v .ge.     10.0 .and. v .lt.    100.0 ) kp = 3
         if ( v .ge.    100.0 .and. v .lt.   1000.0 ) kp = 2
         if ( v .ge.   1000.0 .and. v .lt.  10000.0 ) kp = 1
         if ( v .ge.  10000.0 .and. v .lt. 100000.0 ) kp = 0
         return
      else if ( key .eq. -2 ) then
         if ( v .gt. -10000.0 .and. v .le.  -1000.0 ) km = 0
         if ( v .gt.  -1000.0 .and. v .le.   -100.0 ) km = 1
         if ( v .gt.   -100.0 .and. v .le.    -10.0 ) km = 2
         if ( v .gt.    -10.0 .and. v .lt.      0.0 ) km = 3
         if ( v .ge.      0.0 .and. v .lt.     10.0 ) km = 4
         if ( v .ge.     10.0 .and. v .lt.    100.0 ) km = 3
         if ( v .ge.    100.0 .and. v .lt.   1000.0 ) km = 2
         if ( v .ge.   1000.0 .and. v .lt.  10000.0 ) km = 1
         if ( v .ge.  10000.0 .and. v .lt. 100000.0 ) km = 0
         return
      end if
      if ( key .eq. 0 ) then
         k = min( kp, km )
         k = min( k, 2 )
         if ( k .eq. 4 ) write (tick,10) v
         if ( k .eq. 3 ) write (tick,20) v
         if ( k .eq. 2 ) write (tick,30) v
         if ( k .eq. 1 ) write (tick,40) v
         if ( k .eq. 0 ) write (tick,50) v
   10    format (f6.4)
   20    format (f6.3)
   30    format (f6.2)
   40    format (f6.1)
   50    format (f6.0)
      end if
      return
      end
