C
C    Routine to draw a period axis, based on an adapted version of PGBOX
C

      SUBROUTINE PAXIS (FLO, FHI)                        
C
C-----------------------------------------------------------------------
      INCLUDE  'pgplot.inc'
      CHARACTER*20  CLBL
      LOGICAL  IRANGE, MAJOR, PGNOTO, DEBUG
      INTEGER  I, I1, I2, J, NC, NP, NV, KI
      INTEGER  NSUBX, JMAX, XNFORM
      REAL     TIKL, TIKL1, TIKL2, XC
      REAL     XINT, XINT2, XVAL
      REAL     PGRND
      REAL     A, B, C
      REAL     XNDSP, XMDSP 
      INTRINSIC ABS, INDEX, INT, LOG10, MAX, MIN, MOD, NINT, SIGN, REAL
      REAL     XBLC, XTRC, YBLC, YTRC
C
      IRANGE(A,B,C) = (A.LE.B.AND.B.LE.C) .OR. (C.LE.B.AND.B.LE.A)
C
      DEBUG = .FALSE.
      NXSUB = 1
      IF (PGNOTO('PGBOX')) RETURN
      CALL PGBBUF
      CALL PGQWIN(XBLC, XTRC, YBLC, YTRC)
C
      XNFORM = 1
C
C Displacement of labels from edge of box
C (for X bottom/top, Y left/right, and Y left/right with V option).
C
      XNDSP = 1.2
      XMDSP = 0.7
      YNDSP = 0.7
      YMDSP = 1.2
      YNVDSP = 0.7
      YMVDSP = 0.7
C
C Remove window.
C
      CALL GRAREA(PGID,0.,0.,-1.,-1.)
C
C Draw line for period axis.
C
      CALL GRMOVA(XTRC, YTRC)
      CALL GRLINA(XBLC, YTRC)
C
C Length of X tick marks.
C
      TIKL1 = PGXSP(PGID)*0.6*(YTRC-YBLC)/PGYLEN(PGID)
      TIKL2 = TIKL1*0.5
C
C Start at high frequency end and loop until the entire  
C period axis is annotated.    
C
      FREQ = FHI
      IF(FREQ .LE. FLO .OR. FREQ .LE. 0.0)GO TO 99
      PER = 1.0/(FREQ/1.0e6)           ! Frequencies are in microHz 
C
C Choose the frequency interval, FINT, and corresponding period 
C interval, PINT.
C
      FINT = MAX(0.05, 5.0*PGXSP(PGID)/PGXLEN(PGID))
     1           *(XTRC-XBLC)

      IF(DEBUG)PRINT 9,FLO,FHI,FREQ,PER,FINT,PGXSP(PGID),PGXLEN(PGID),
     *                 XTRC,XBLC
    9 FORMAT(' FLO= ',F10.1,'  FHI=',F10.1,'  FREQ=',F10.1,'   PER=',
     *       F10.2/
     *       ' FINT=',F10.1,'  XSP=',F10.3,'  XLEN=',F10.3,'  XTRC=',
     *       F10.2,'  XBLC=',F10.1//)

      PEROLD = 0.0
   10 PINT = FINT/1.0e6 * PER * PER
      PINT = RND(PINT)

      IF(DEBUG)PRINT 11,PINT,PER,PINT*NINT(PER/PINT),
     *                  1.E6/(PINT*NINT(PER/PINT))
   11 FORMAT(' PINT=',F7.1,'  PER=',F8.2,'  PER''',F8.2,'  FREQ=',F10.1)

C
C Now make PER an integral multiple of PINT
C
      PER = PINT * NINT(PER/PINT)
      IF(PER .EQ. 0.0)GO TO 99
      IF(PER.LT.PEROLD+PINT)PER = PEROLD+PINT
      FREQ = 1.0e6 * 1.0/PER
      IF(FREQ .LE. FLO .OR. FREQ .LE. 0.0)GO TO 99
C
C Check FREQ is >= FLO and <= FHI
C
   20 IF(FREQ .GT. FHI)THEN
        PER = PER + PINT
        FREQ = 1.0e6 * 1.0/PER
        IF(FREQ .LE. FLO .OR. FREQ .LE. 0.0)GO TO 99
        GO TO 20
      ELSE IF(FREQ .LT. FLO)THEN
        GO TO 99
      END IF
C
C Split PER into mantissa and exponent 
C
      NP = INT(LOG10(ABS(PER)))-4
      NV = NINT(PER/10.**NP)
C
C Draw X ticks and write the X label
C
      MAJOR = .TRUE.                        
      TIKL = TIKL2
      IF (MAJOR) TIKL = TIKL1
      XVAL = FREQ   
      IF (IRANGE(XBLC,XVAL,XTRC)) THEN
          CALL GRMOVA(XVAL, YTRC-TIKL)
          CALL GRLINA(XVAL, YTRC)
      END IF
      XC = (XVAL-XBLC)/(XTRC-XBLC)
      CALL PGNUMB(NV,NP,XNFORM,CLBL,NC)
      CALL PGMTXT('T', XMDSP, XC, 0.5, CLBL(1:NC))

C
C     That's the first tick mark and label done. Now decrement
C     FREQ, recalculate the intervals and repeat
C
      PEROLD = PER
      FREQ = FREQ - FINT
      IF(FREQ .LE. FLO .OR. FREQ .LE. 0.0)GO TO 99
      PER = 1./(FREQ/1.0e6)
      GO TO 10

   99 CONTINUE
C
C     Restore window: interior of box and return
C
      CALL GRAREA(PGID,PGXOFF(PGID),PGYOFF(PGID),
     1            PGXLEN(PGID),PGYLEN(PGID))
      CALL PGEBUF
      RETURN
      END


C
C Routine to find the nearest "round" number to x, a
C "round" number being 1, 2 or 5 times a power of 10. If X is negative,
C RND(X) = -RND(ABS(X)). eg RND(8.7) = 10.0, RND(-0.4) = -0.5.  
C If X is zero, the value returned is zero.
C
      REAL FUNCTION RND (X)
      REAL X
C
C-----------------------------------------------------------------------
      INTEGER  I,ILOG
      REAL     FRAC,NICE(4),PWR,XLOG,XX
      INTRINSIC ABS, LOG10, SIGN
      DATA     NICE/1.0,2.0,5.0,10.0/
C
      IF (X.EQ.0.0) THEN
          RND = 0.0
          RETURN
      END IF
      XX   = ABS(X)
      XLOG = LOG10(XX)
      ILOG = XLOG
      IF (XLOG.LT.0) ILOG=ILOG-1
      PWR  = 10.0**ILOG
      FRAC = XX/PWR
      I = 4
      IF (FRAC.LE.7.0) I = 3
      IF (FRAC.LE.3.0) I = 2
      IF (FRAC.LE.1.5) I = 1
      RND = SIGN(PWR*NICE(I),X)
      END
