C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/ticks.f90
C
      SUBROUTINE TICKS2(XMIN,XMAX,NTICKS,XMINP,XMAXP,DIST,LABS,IDIML,
     *                 MBEG,MPOINT,MEND,ENOTA,ME,ISCALE,NSIG,NDEC,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
C==========GIVEN RAW DATA MINIMUM AND MAXIMUM VALUES AND THE NUMBER OF
C          DESIRED GRID INTERVALS, TICKS2 RETURNS "NICE" CHARACTER
C          TICK MARK LABELS AND ASSOCIATED SCALING INFORMATION FOR
C          LABELING LINEAR GRAPH AXES.  TICK MARK LABELS  WILL  BE
C          ACCURATE  UP  TO  MAXSIG=12  SIGNIFICANT DIGITS.  GIVEN
C          NSIG, THE NUMBER OF SIGNIFICANT DIGITS ASSOCIATED  WITH
C          THE  RAW DATA, AND NDEC, THE NUMBER OF PLACES AFTER THE
C          DECIMAL POINT, CHARACTER LABELS WILL BE FORMATTED  WITH
C          THE  FORTRAN FORMAT F(NSIG+2).(NDEC) IF THEY FIT SUCH A
C          FORMAT  WITHOUT  LOOSING  EITHER  HIGH   OR   LOW-ORDER
C          SIGNIFICANCE;  OTHERWISE  THEY  WILL  BE FORMATTED WITH
C          1PE(NSIG+6).(NSIG-1).  AN  OPTION  (SEE  ENOTA,  ISCALE
C          BELOW)  IS PROVIDED TO ALLOW POSSIBLE EXPONENT NOTATION
C          TO BE ELIMINATED FROM INDIVIDUAL LABELS; THAT IS,  EACH
C          LABEL  WILL BE OF CONSISTENT LENGTH AND SCALE SUCH THAT
C          ONE SCALE FACTOR (10**ISCALE) HOLDS FOR ALL.  AN OPTION
C          IS ALSO PROVIDED (SEE MBEG,MEND) TO MINIMIZE THE LENGTH
C          OF THE CHARACTER LABELS BY ELIMINATING  LEADING  BLANKS
C          AND  UNNECESSARY  TRAILING ZEROS AND DECIMAL POINT.  IF
C          NSIG WOULD CAUSE SIGNIFICANT DIGITS OF ANY LABEL TO  BE
C          DROPPED  BY  THE 1PE FORMAT ABOVE, THEN SUCH LABELS ARE
C          INACCURATE; INACCURATE LABELS ARE RETURNED AS A  STRING
C          OF  ASTERISKS(*'S)  AND AN ERROR INDICATOR (SEE IER) IS
C          SET.  LIKEWISE AN ERROR INDICATOR IS SET IF XMIN,  XMAX
C          ARE  TOO  CLOSE  TO OVERFLOW OR UNDERFLOW LIMITS; OR IF
C          THE EXPRESSION XMAX-XMIN, NEEDED FOR SCALING  THE  DATA
C          FOR   PLOTTING,   IS  SO  LARGE  THAT  LESS  THAN  NSIG
C          SIGNIFICANT DIGITS OF EITHER OPERAND PARTICIPATE IN THE
C          SUBTRACTION; OR IF XMAX-XMIN IS SO SMALL THAT LESS THAN
C          NSIG DIGITS REMAIN AFTER CANCELLING.
C
C-----PARAMETERS:
C     XMIN,XMAX  ARE GIVEN AS MINIMUM AND MAXIMUM DATA VALUES.
C                XMIN<XMAX.
C
C     NTICKS IS GIVEN AS THE NUMBER OF TICK MARK LABELS REQUESTED.
C            NTICKS-1 IS THE NUMBER OF GRID(AXIS) INTERVALS.
C
C     XMINP,XMAXP  ARE RETURNED AS "NICE" MINIMUM AND MAXIMUM ENDPOINTS
C                  SUCH THAT BOTH ARE ESSENTIALLY INTEGER MULTIPLES OF
C                  OF DIST; THAT IS, NTICKS-1=FLOOR((XMAXP-XMINP)/DIST)
C                  =CEIL((XMAXP-XMINP)/DIST).  ALSO XMINP.LE.XMIN AND
C                  XMAXP.GE.XMAX.  ASSUMING A CONSISTENT UNIT OF
C                  MEASUREMENT, RAW DATA, X, MAY BE SCALED TO X' FOR
C                  PLOTTING BY APPLYING THE FORMULA:
C                  X'=(X-XMINP)*(AXIS LENGTH)/(XMAXP-XMINP) OR
C                  X'=(X-XMINP)*(AXIS LENGTH)/((NTICKS-1)*DIST) .
C
C     DIST  IS RETURNED AS A "NICE" GRID INTERVAL;  THAT IS DIST WILL
C           BE WITHIN FUZZ OF 1, 2, OR 5 TIMES A POWER OF 10 A LA ACM
C           ALGORITHM 463, SCALE2.
C
C     LABS  IS RETURNED AS THE RIGHT-JUSTIFIED CHARACTER REPRESENTATION
C           OF THE REAL NUMBERS XMINP+(I-1)*DIST FOR I=1,2,...,NTICKS.
C           EACH LABEL LABS(J,I),J=1,2,...,IDIML, WILL BE ROUNDED TO
C           MAXSIG SIGNIFICANT DIGITS AND WILL NOT CONTAIN EXPLICIT
C           EXPONENT NOTATION (SEE ISCALE AND ENOTA BELOW).  WHEN IER=0
C           EACH NUMERIC LABEL WILL BE A CHARACTER STRING AS FOLLOWS:
C           FOR I=1,2,...,NTICKS,
C           LABS(1,2,...,IDIML, I)='BB...SDD...DPDD...D00...0'    WHERE
C           BB... ARE MBEG-1 POSSIBLY UNNECESSARY LEADING BLANKS.
C           S     IS THE SIGN, BLANK(' ') OR MINUS('-').
C           DD... ARE NSIG DIGITS.
C           P     IS A DECIMAL POINT AT LABS(MPOINT,I).
C           00... ARE IDIML-MEND POSSIBLY UNNECESSARY TRAILING ZEROS.
C
C     IDIML  IS GIVEN AS THE FIRST DIMENSION OF LABS. IDIML MUST BE AT
C            LEAST AS LARGE AS NSIG+2 (TO ALLOW FOR A DECIMAL POINT AND
C            SIGN).   IDIML MUST NEVER BE GREATER THAN 32.
C
C     MBEG, MPOINT, MEND  ARE RETURNED AS POINTERS SUCH THAT FOR ALL I:
C            LABS(MBEG,I) IS THE HIGHEST-ORDER TRUE DIGIT OR SIGN.
C            LABS(MPOINT,I) IS A DECIMAL POINT('.').
C            LABS(MEND,I) IS THE LOWEST-ORDER NON-ZERO DIGIT.
C            THE NUMBER OF CHARACTERS IN THE TICK  MARK  LABELS  CAN  BE
C            MINIMIZED  WHILE  MAINTAINING  CONSISTENT  LENGTH LABELS BY
C            LABELING WITH ((LABS(J,I),J=MBEG,MEND),I=1,NTICKS), OR
C            OUTPUT EXACTLY AS FORMATED BY LABELING WITH
C            ((LABS(J,I),J=MBEG,IDIML),I=1,NTICKS).
C
C     ENOTA  IS RETURNED AS A 4 CHARACTER STRING REPRESENTING 10**ISCALE
C            IN  SCIENTIFIC  E-NOTATION.   THAT  IS,   THE   LAYOUT   IS
C            ENOTA='ESDD'  IF ME=4, OR ENOTA='ESDB' IF ME=3.  ENOTA WILL
C            CONTAIN  ME  SIGNIFICANT  CHARACTERS  AND  4-ME  TRAILING
C            BLANKS. WHEN ISCALE=0 ENOTA='E+0 '.
C
C     ME  IS RETURNED AS THE NUMBER OF SIGNIFICANT CHARACTERS IN ENOTA.
C         ME=3 OR ME=4.
C
C     ISCALE  IS RETURNED AS THE INTEGER VALUED EXPONENT OF 10; EACH
C           EACH LABEL IS ADJUSTED SO THAT THEY ALL HAVE THE COMMON
C           SCALE FACTOR, ISCALE.
C
C     NSIG  IS GIVEN AS THE PRECISION OF THE RAW DATA, THAT IS THE
C           NUMBER OF SIGNIFICANT DIGITS REPRESENTING THE RAW DATA. NSIG
C           MUST BE GREATER THAN 1 AND LESS THAN 13.
C
C     NDEC  IS GIVEN AS THE NUMBER OF PLACES AFTER THE DECIMAL POINT.
C           NDEC MUST BE .GE.0 AND .LE. NSIG.
C
C     IER  IS RETURNED AS AN ERROR INDICATOR:
C          IER=0  FOR NORMAL RETURN.
C          IER=1  FOR XMIN.GE.XMAX OR NTICKS.LT.2 .
C          IER=2  FOR NSIG.LT.1 OR NSIG.GT.12, OR NDEC.LT.0 OR
C                 NDEC.GT.NSIG .
C          IER=3  FOR IDIML .LT. NSIG+2, OR IDIML .GT. 32.
C          IER=4  FOR XMIN, XMAX TOO CLOSE TO OVERFLOW LIMITS.
C          IER=5  FOR XMIN, XMAX TOO CLOSE TO UNDERFLOW LIMITS.
C          IER=6  IF THE DATA IS SPREAD (XMIN-XMAX) TOO FAR OR TOO CLOSE
C                 SUCH THAT SUBSEQUENT SCALING FOR PLOTTING WOULD CAUSE
C                 AN NSIG (OR MORE) DECIMAL DIGIT PRECISION LOSS.
C                 (E.G. FOR MACH10=16: LET XMIN=1, XMAX=1.E17,NSIG=3 OR
C                 LET XMIN=1.0000, XMAX=1.0001,NTICKS=10,NSIG=12)
C          IER=7  IF NSIG CAUSES SIGNIFICANT DIGITS TO BE TRUNCATED FOR
C                 FOR ANY LABEL.  SUCH LABELS WILL BE RETURNED AS IDIML
C                 ASTERISKS(E.G. XMINP=1.0, DIST=0.1, NSIG=1,
C                 IDIML=10 WILL YIELD **********).
C          IF IER IS RETURNED 1,2 3,4,5, OR 6 THEN OTHER RETURNED
C          ARGUMENT VALUES WILL BE UNDEFINED.
C
C     METHOD: A DOUBLE PRECISION VERSION OF ACM ALGORITHM 463 BY C. R.
C             LEWART IS USED TO COMPUTE NICE NUMERIC ENDPOINTS AND
C             GRID INTERVAL.  IF NOT CLEARLY IMPOSSIBLE, THE LABELS ARE
C             THEN FORMATTED F(NSIG+2).(NDEC).  IF ANY LABEL LOOSES
C             SIGNIFICANCE THEN ALL LABELS ARE REFORMATTED TO
C             1PE(NSIG+6).(NSIG-1) WITH FOUR OF THE NSIG+6 CHARACTERS
C             REPRESENTING THE EXPONENT STORED SEPARATELY IN ENOTA.
C             LABELS CAN THEN BE OUTPUT FOR AN AXIS WITH OR WITHOUT
C             INDIVIDUAL EXPONENT NOTATION; IN THE LATER CASE, WHEN
C             LABEL SPACE IS AT A PREMIUM, THE AXIS TITLE NEEDS TO BE
C             CONCATENATED WITH A CHARACTER STRING EQUIVALENT TO THE
C             PHRASE "(TIMES 10**ISCALE)".
C
C-----MACHINE DEPENDENT ROUTINES: CEQ/CNE, CHMOVE, DPTOCH, INTOCH.
C     MACHINE PORTABLE ROUTINES: DCHECK, SCALE2, TFLOOR, TCEIL, D1MACH.
C
C  H. D. KNOBLE  - JANUARY 1979  - PENN STATE COMPUTATION CENTER.
C
      LOGICAL*1 LABS(IDIML,NTICKS),ENOTA(4),WORK(26)
C                                                |_2*MAXSIG+2
      LOGICAL CEQ,CNE, P
C---------MAXEXP IS THE MAXIMUM ABSOLUTE BASE 10 EXPONENT ALLOWED;
C         MACH10 IS THE MAXIMUM NUMBER OF BASE 10
C                FLOATING-POINT SIGNIFICANT DIGITS. THAT IS
C                MACH10=FLOOR(LOG10(2**(M-2)+.5)) WHERE M IS THE
C                NUMBER OF BITS IN THE FLOATING-POINT MANTISSA.
C         CT IS THE COMPARISON TOLERANCE (SEE TFLOOR).
      DATA MAXEXP/75/, MACH10/16/
      MAXSIG=MACH10-4
      CT=3.D0*D1MACH(54321)
      IER=0
C--------CHECK INPUT PARAMETERS.
      IF(XMIN.GE.XMAX.OR.NTICKS.LT.2) IER=1
      IF(NSIG.LT.1.OR.NSIG.GT.MAXSIG.OR.NDEC.LT.0.OR.NDEC.GT.NSIG)IER=2
      IF(IDIML.LT.NSIG+2 .OR. IDIML.GT.32)IER=3
      IF(IER.NE.0) RETURN
C-----------CHECK FOR POSSIBLE PRECISION LOSS DURING SCALE2.
      CALL DCHECK(XMAX,XMIN,MACH10,NSIG,IER)
      IER=IER*6
      IF(IER.NE.0) RETURN
C--------GET "NICE" MIN, MAX AND GRID INTERVAL A LA ACM
C        ALGORITHM 463.
      CALL SCALE2(XMIN,XMAX,NTICKS-1,XMINP,XMAXP,DIST,CT)
C----------CHECK FOR PRECISION LOSS DURING SUBSEQUENT PLOT SCALING.
      CALL DCHECK(XMAXP,XMINP,MACH10,NSIG,IER)
      IER=IER*6
      IF(IER.NE.0) RETURN
C-------SCALE "NICE" AXIS PARAMETERS.
      POW1=0.D0
      POW2=0.D0
      IF(XMINP.NE.0.D0) POW1=DLOG10(DABS(XMINP))
      IF(XMAXP.NE.0.D0) POW2=DLOG10(DABS(XMAXP))
      POW10=DMAX1(POW1,POW2)
      MAGP=TFLOOR(POW10,CT)
C------COMPUTE SCALE TO FIT AXIS PARAMETERS IN F(NSIG+2).(NDEC).
      P=.FALSE.
      KDEC=NDEC
      ISCALE=0
      IF((MAGP.GE.0.AND.MAGP.LT.NSIG-NDEC).OR.
     *   (MAGP.LT.0.AND.MAGP.GE.-NDEC))       GO TO 45
C---------COMPUTE SCALE TO FIT AXIS PARAMETERS IN
C         1PE(NSIG+6).(NSIG-1).
42    ISCALE=MAGP
      IER=0
      P=.TRUE.
      KDEC=NSIG-1
      IF(IABS(ISCALE).LE.MAXEXP) GO TO 45
C--------CHECK FOR POSSIBLE UNDERFLOW/OVERFLOW.
      IF(ISCALE.GE.MAXEXP) IER=4
      IF(ISCALE.LE.-MAXEXP) IER=5
      RETURN
C-----------APPLY FINAL SCALE TO AXIS PRARMETERS.
45    SCLFAC=10.D0**(-ISCALE)
      SDIST=DIST*SCLFAC
      SMINP=XMINP*SCLFAC
      SMAXP=XMAXP*SCLFAC
      KWEND=2*MAXSIG+2
      NS2=NSIG+2
      IRIGHT=IDIML-NS2
C-----------BEGIN BUILDING CHARACTER LABELS.
      DO 50 I=1,NTICKS
      XLAB=(I-1)*SDIST+SMINP
C----------XLAB CLOSE TO ZERO RELATIVE TO 10**ISCALE MAY CAUSE
C          FORMAT PROCESSING TO GET *'S.
      IF(TFLOOR(XLAB,CT).EQ.TCEIL(XLAB,CT)) XLAB=TFLOOR(XLAB,CT)
C-------CALCULATE KW, KD AND USE F(KW.KD) FORMAT TO ELIMINATE
C       FUZZ ON END OF NUMBERS.
      MAG=0
      IF(XLAB.NE.0.D0) MAG=TFLOOR(DLOG10(DABS(XLAB)),CT)
      IF(MAG.GE.0) MAG=MAG+1
      KD=MAXSIG-MAG
      KW=KWEND
C---------CONVERT TICK VALUES TO CHARACTER LABELS.
55    CALL DPTOCH(XLAB,WORK,1,KW,KD,IERR)
      IF(IERR.NE.0) GO TO 90
C-------- MAKE SURE PRECISION WAS NOT LOST BY NSIG  CHOICE.
      NB=KW-KD+KDEC+1
      NE=KW
      IF(NB.GT.NE) GO TO 100
      DO 60 K=NB,NE
      IF(CNE(WORK(K),'0',1)) GO TO 90
60    CONTINUE
      GO TO 100
C----------MORE THAN NSIG SIGNIFICANT DIGITS - SET LABELS TO
C          ASTERISKS.
90    IER=7
      IF(P)
     *CALL CHMOVE('*********************************',LABS(1,I),IDIML)
      GO TO 50
100   CALL CHMOVE('                                 ',LABS(1,I),IDIML)
C--------RIGHT-JUSTIFY XLAB IN LABS(1,2,...,IDIML, I).
      CALL CHMOVE(WORK(NB-NS2),LABS(1+IRIGHT,I),NS2)
50    CONTINUE
      IF(IER.EQ.7.AND. .NOT. P) GO TO 42
C---------COMPUTE MPOINT, POINTER TO DECIMAL POINT.
      MPOINT=NS2-KDEC+IRIGHT
C---------COMPUTE MEND, POINTER TO LAST TRAILING DIGIT AFTER THE
C         DECIMAL POINT WHICH IS NOT AN UNNECESSARY ZERO.
      MEND=MPOINT+KDEC
      NB=MPOINT+1
      NE=MEND
      IF(NB.GT.NE) GO TO 115
      DO 110 J=NB,NE
      DO 112 I=1,NTICKS
      IF(CNE(LABS(MEND,I),'0',1)) GO TO 120
112   CONTINUE
      MEND=MEND-1
110   CONTINUE
C-------FIX MEND TO ENABLE DECIMAL POINT TO BE IGNORED FOR WHOLE
C       NUMBERS.
115   IF(MEND.EQ.MPOINT) MEND=MEND-1
C--------COMPUTE MBEG, POINTER TO FIRST LEADING SIGNIFICANT DIGIT OR
C        SIGN.
120   MPM1=MPOINT-1
      DO 130 MBEG=1,MPM1
      DO 132 I=1,NTICKS
      IF(CNE(LABS(MBEG,I),' ',1)) GO TO 140
132   CONTINUE
130   CONTINUE
C-----------FORMAT ENOTA, THE EXPONENT PORTION.
140   CALL CHMOVE('E+0 ',ENOTA,4)
      ME=3
      IF(ISCALE.EQ.0) RETURN
      IF(ISCALE.LT.0)CALL CHMOVE('-',ENOTA(2),1)
C-------CONVERT INTEGER SCALE TO CHARACTER.
      CALL INTOCH(IABS(ISCALE),ENOTA(3),1,IERR)
      IF(IERR.EQ.0) RETURN
      ME=ME+1
      CALL INTOCH(IABS(ISCALE),ENOTA(3),2,IERR)
      RETURN
      END
      SUBROUTINE DCHECK(XMAX,XMIN,MACH10,NSIG,IER)
C=========DCHECK CHECKS THE DIFFERENCE XMAX-XMIN FOR PRECISION LOSS
C         AND RETURNS IER=1 IF:
C  (A) LESS THAN NSIG SIGNIFICANT DECIMAL DIGITS OF EITHER OPERAND
C      PARTICIPATE IN THE OPERATION; THAT IS XMAX IS VERY LARGE RELATIVE
C      TO XMIN.
C  (B) MACH10-NSIG OR MORE DECIMAL DIGITS CANCEL DURING THE SUBTRACTION;
C      THAT IS XMAX IS ABOUT EQUAL TO XMIN RELATIVE TO NSIG.
C  IER=0 OTHERWISE.
C
C  MACH10 IS GIVEN THE NUMBER OF DECIMAL DIGITS CARRIED IN DOUBLE
C         PRECISION FOR THE HOST COMPUTER. MACH10 MAY BE COMPUTED AS:
C         FLOOR(DLOG10(2.D0**(M-2)+.5D0)) WHERE M IS THE NUMBER OF BITS
C         IN THE FLOATING-POINT MANTISSA (SEE MATULA, D. W., "IN-AND-OUT
C         CONVERSIONS, CACM 11(1):47-50, JANUARY 1968.).
C  NSIG IS GIVEN AS THE DECIMAL PRECISION OF XMAX,XMIN.
C
C    H. D. KNOBLE  -  DECEMBER 1972  - PENN STATE COMPUTATION CENTER.
C
      DOUBLE PRECISION XMAX,XMIN
      IER=0
      IF(XMAX.EQ.XMIN.OR.XMAX.EQ.0.OR.XMIN.EQ.0) RETURN
C------(A)
      IF(DABS(1.D0+((XMAX-XMIN)-XMAX)/XMIN).GT.10.D0**(1-NSIG))IER=1
C------(B)
      IF(DABS((XMAX-XMIN)/XMAX).LT.10.D0**(-(MACH10-(NSIG+1))))IER=1
      RETURN
      END
      SUBROUTINE SCALE2 (XMIN,XMAX,N,XMINP,XMAXP,DIST,CT)
C==========REFERENCE: ACM ALGOIRTHM 463, ALGORITHMS SCALE1
C                     SCALE2, AND SCALE3 FOR DETERMINATION OF SCALES
C                     ON COMPUTER GENERATED PLOTS, BY C. R. LEWART,
C                     CACM 16(10):639-640, OCTOBER 1973.
C
C=========GIVEN XMIN,XMAX AND N SCALE2 FINDS A NEW RANGE XMINP
C         AND XMAXP DIVISIBLE INTO EXACTLY N LINEAR INTERVALS OF SIZE
C         DIST, WHERE N IS GREATER THAN 1.
C         DIST WILL BE A "NICE" NUMBER, THAT IS 1, 2, OR 5 TIMES
C         A POWER OF 10. THE CLOSED INTERVAL [XMINP,XMAXP] WILL
C         CONTAIN THE ENDPOINTS XMIN, XMAX.
C
C-------THIS IS A DOUBLE PRECISION VERSION AND WILL SUPPORT UP
C       TO 12 SIGNIFICANT DIGITS ON THE IBM SERIES.
C
C   MODIFIED TO USE FUZZY FLOOR.
C   H. D. KNOBLE  - JANUARY 1979  - PENN STATE COMPUTATION CENTER.
C
      IMPLICIT REAL*8 (A-H, O-Z,$)
      DIMENSION VINT(5)
      DOUBLE PRECISION CT,M1,M2,NP,NX
      DATA VINT(1),VINT(2),VINT(3),VINT(4),VINT(5)/
     *      1.D0,   2.D0,   5.D0,   10.D0,  20.D0 /
   10 DEL=0.000000000000002D0
      FN=N
C--------FIND APPROXIMATE INTERVAL SIZE, A.
      A=(XMAX-XMIN)/FN
      AL=DLOG10(A)
      NAL=AL
      IF (A.LT.1.0D0) NAL=NAL-1
C----------A IS SCALED INTO A VARIABLE NAMED B BEWTEEN 1 AND 10.
      B=A/10.0D0**NAL
C-----------COMPUTE THE CLOSEST PERMISSIBLE VALUE FOR B.
      DO 20  I=1,3
      IF (B.LT.(VINT(I)+DEL)) GO TO 30
   20 CONTINUE
      I=4
C-------COMPUTE THE INTERVAL SIZE.
   30 DIST=VINT(I)*10.0D0**NAL
      FM1=XMIN/DIST
      M1=TFLOOR(FM1,CT)
      IF(DABS(M1+1.-FM1).LT.DEL)M1=M1+1.
C---------COMPUTE THE NEW MINIMUM AND MAXIMUM LIMITS.
      XMINP=DIST*M1
      FM2=XMAX/DIST
      M2=TFLOOR(FM2+1.0D0,CT)
      IF (DABS(FM2+1.0D0-M2).LT.DEL) M2=M2-1.0D0
      XMAXP=DIST*M2
C----------CHECK WHETHER A SECOND PASS IS NEEDED.
      NP=M2-M1
      IF (NP.LE.DFLOAT(N)) GO TO 40
      I=I+1
      GO TO 30
   40 NX=TFLOOR((DFLOAT(N)-NP)/2.0D0,CT)
      XMINP=XMINP-NX*DIST
      XMAXP=XMINP+DFLOAT(N)*DIST
C---------ADJUST LIMITS TO ACCOUNT FOR ROUND-
      IF (XMINP.GT.XMIN) XMINP=XMIN
      IF (XMAXP.LT.XMAX) XMAXP=XMAX
      RETURN
      END
      DOUBLE PRECISION FUNCTION TFLOOR(X,CT)
C===========TOLERANT FLOOR FUNCTION.
C
C    X  -  IS GIVEN AS A DOUBLE PRECISION ARGUMENT TO BE OPERATED ON.
C          IT IS ASSUMED THAT X IS REPRESENTED WITH M MANTISSA BITS.
C    CT -  IS   GIVEN   AS   A   COMPARISON   TOLERANCE   SUCH   THAT
C          0.LT.CT.LE.3-SQRT(5)/2. IF THE RELATIVE DIFFERENCE BEWTEEN
C          X AND A WHOLE NUMBER IS  LESS  THAN  CT,  THEN  TFLOOR  IS
C          RETURNED   AS   THIS   WHOLE   NUMBER.   BY  TREATING  THE
C          FLOATING-POINT NUMBERS AS A FINITE ORDERED SET  NOTE  THAT
C          THE  HEURISTIC  EPS=2.**(-(M-1))   AND   CT=3*EPS   CAUSES
C          ARGUMENTS  OF  TFLOOR/TCEIL TO BE TREATED AS WHOLE NUMBERS
C          IF THEY ARE  EXACTLY  WHOLE  NUMBERS  OR  ARE  IMMEDIATELY
C          ADJACENT TO WHOLE NUMBER REPRESENTATIONS.  SINCE EPS,  THE
C          "DISTANCE"  BETWEEN  FLOATING-POINT  NUMBERS  ON  THE UNIT
C          INTERVAL, AND M, THE NUMBER OF BITS IN X'S MANTISSA, EXIST
C          ON  EVERY  FLOATING-POINT   COMPUTER,   TFLOOR/TCEIL   ARE
C          CONSISTENTLY DEFINABLE ON EVERY FLOATING-POINT COMPUTER.
C
C          FOR MORE INFORMATION SEE THE FOLLOWING REFERENCES:
C    [1] P. E. HAGERTY, "MORE ON FUZZY FLOOR AND CEILING," APL  QUOTE
C        QUAD 8(4):20-24, JUNE 1978. NOTE THAT TFLOOR=FL5.
C    [2] L. M. BREED, "DEFINITIONS FOR FUZZY FLOOR AND CEILING",  APL
C        QUOTE QUAD 8(3):16-23, MARCH 1978.
C
C   H. D. KNOBLE, PENN STATE UNIVERSITY.
C=====================================================================
      DOUBLE PRECISION X,Q,RMAX,EPS5,CT,FLOOR,DINT
C---------FLOOR(X) IS THE LARGEST INTEGER ALGEBRAICALLY LESS THAN
C         OR EQUAL TO X; THAT IS, THE UNFUZZY FLOOR FUNCTION.
      DINT(X)=X-DMOD(X,1.D0)
      FLOOR(X)=DINT(X)-DMOD(2.D0+DSIGN(1.D0,X),3.D0)
C---------HAGERTY'S FL5 FUNCTION FOLLOWS...
      Q=1.D0
      IF(X.LT.0)Q=1.D0-CT
      RMAX=Q/(2.D0-CT)
      EPS5=CT/Q
      TFLOOR=FLOOR(X+DMAX1(CT,DMIN1(RMAX,EPS5*DABS(1.D0+FLOOR(X)))))
      IF(X.LE.0 .OR. (TFLOOR-X).LT.RMAX)RETURN
      TFLOOR=TFLOOR-1.D0
      RETURN
      END
      DOUBLE PRECISION FUNCTION TCEIL(X,CT)
C==========TOLERANT CEILING FUNCTION.
C    SEE TFLOOR.
      DOUBLE PRECISION X,CT,TFLOOR
      TCEIL= -TFLOOR(-X,CT)
      RETURN
      END
      DOUBLE PRECISION FUNCTION D1MACH (IDUM)
      INTEGER IDUM
C=======================================================================
C THIS ROUTINE COMPUTES THE UNIT ROUNDOFF OF THE MACHINE IN DOUBLE
C PRECISION.  THIS IS DEFINED AS THE SMALLEST POSITIVE MACHINE NUMBER
C EPS SUCH THAT (1.0D0+EPS > 1.0D0) & (1.D0-EPS < 1.D0).
C
C THIS COMPUTATION OF EPS IS THE WORK OF ALAN C. HINDMARSH,
C LAWRENCE LIVERMORE NATIONAL LAB.
C-----------------------------------------------------------------------
      DOUBLE PRECISION EPS, COMP
      EPS = 1.0D0
 10   EPS = EPS*0.5D0
      COMP = 1.0D0 + EPS
      IF (COMP .NE. 1.0D0) GO TO 10
      D1MACH = EPS*2.0D0
      RETURN
      END
      IMPLICIT REAL*8 (A-H,O-Z)
C
C===========SHORT TEST DRIVER FOR SUBROUTINE TICKS2.
C
      LOGICAL*1 LABS(14,20),ENOTA(4),TITLE(80)
      IDIML=14
C--------INPUT A TEST CASE.
55    READ(5,10,END=99) TITLE
10    FORMAT(80A1)
      WRITE(6,11) TITLE
11    FORMAT(////80A1)
      READ(5,*) XMIN,XMAX,N,NSIG,NDEC
      WRITE(6,1) XMIN,XMAX,N,NSIG,NDEC
1     FORMAT('0XMIN =',G20.12,'  XMAX =',G20.12,' N=',I3,' (NSIG,NDEC)='
     *       ,'(',I2,',',I2,')')
C--------CHECK OUT TICKS2.
      CALL TICKS2(XMIN,XMAX,N,XMINP,XMAXP,DIST,LABS,14,MBEG,MPOINT,
     *           MEND,ENOTA,ME,ISCALE,NSIG,NDEC,IER)
      IF(IER.NE.0 .AND. IER.NE.7) GO TO 7
C--------OUTPUT RESULTS.
      WRITE(6,2) XMINP,XMAXP,DIST
2     FORMAT('0XMINP=',G20.12,'  XMAXP=',G20.12,'  GRID INTERVAL=',
     *       G20.12)
      WRITE(6,3) MBEG,MEND,MPOINT,ISCALE,(ENOTA(J),J=1,ME)
3     FORMAT('0MBEG,MEND,MPOINT=',3I4,'  TIMES 10**',I3,
     *       '    CHARACTER EXPONENT=',4A1)
C--------OUTPUT FORMATTED LABELS FOLLOWED BY MINIMALLY FORMATED LABELS.
      WRITE(6,9)
9     FORMAT('0',10X,'CHARACTER LABELS'/'   AS REQUESTED',9X,
     *       'OF MINIMAL LENGTH'/'+',40('_'))
      DO 5 I=1,N
5     WRITE(6,4) (LABS(J,I),J=1,IDIML),(LABS(J,I),J=MBEG,MEND)
4     FORMAT(' ->',14A1,5X,'->',14A1)
C--------ERROR RETURN FROM TICKS2.
      IF(IER.EQ.0) GO TO 55
7     WRITE(6,8) IER
8     FORMAT('0ERROR RETURN CODE=',I2)
C--------GO TRY ANOTHER TEST CASE.
      GO TO 55
C--------WRAP UP.
99    STOP
      END
//DATA.INPUT DD *
1EXAMPLE 1, TYPICAL AXES PARAMETERS, NEGATIVE TO POSITIVE, SEE CACM.
-3.1 11.1 6 3 1
 EXAMPLE 1, VERY SMALL EXPONENTS.
-3.1D-65 11.1D-65 6 3 1
1EXAMPLE 2, TYPICAL AXES PARAMETERS, ALL POSITIVE, EXAMPLE FROM CACM.
5.2 10.1 6 3 1
 EXAMPLE 2, VERY LARGE EXPONENTS.
5.2D+65 10.1D+65 6  3 1
1EXAMPLE 3, LARGE WHOLE NUMBERS WITH DECREASING NSIG.
-120000000000 -1000000000 10 12 0
 EXAMPLE 3, LARGE WHOLE NUMBERS WITH DECREASING NSIG.
-120000000000 -1000000000 10 11 1
 EXAMPLE 3, LARGE WHOLE NUMBERS WITH DECREASING NSIG.
-120000000000 -1000000000 10  8 1
 EXAMPLE 3, LARGE WHOLE NUMBERS WITH DECREASING NSIG.
-120000000000 -1000000000 10  2 0
1EXAMPLE 4, LOTS OF DECIMAL PLACES WITH DECREASING NSIG.
.000000000012 .0000000001 10 12 12
 EXAMPLE 4, LOTS OF DECIMAL PLACES WITH DECREASING NSIG.
.000000000012 .0000000001 10  8 0
 EXAMPLE 4, LOTS OF DECIMAL PLACES WITH DECREASING NSIG.
.000000000012 .0000000001 10  2 0
1EXAMPLE 5, TYPICAL AXES PARAMETERS, F12.3 FORMAT.
.37 1.92 5  12 3
1EXAMPLE 6, TYPICAL AXES PARAMETERS, F12.2, ELIMINATE DECIMAL POINT.
-25.6 -5.1 7   12 2
1EXAMPLE 7, TEST ERROR CODE 7, WON'T FIT WITHOUT LOOSING SIGNIFICANCE.
1.0 2.0 11 1 1
 EXAMPLE 7, JUST FITS IN 1PE FORMAT.
1.0 2.0 11 2 0
1EXAMPLE 8, TEST ERROR RETURN 6, XMAX-XMIN TOO LARGE.
1.000 1.000E+17 10 4 3
 EXAMPLE 8, TEST ERROR RETURN 6, XMAX-XMIN TOO SMALL.
1.000 1.0001 10 12 3