```c   This program performs calculations
c   for a proposed ASTM duration of load standard.
c
c   Steve Verrill, 5/6/97
c   Steve Verrill, 10/22/97 --- modified to reflect
c   my new understanding of the test protocol
c
c
c

real red(5),cv(5),probf(5,5)
real lnp6
real sigma(5)

real samp1(200),samp2(200)
real sort1(200),sort2(200)

real strrat(200)

real av(5),bv(5)

integer iorigord(200)

integer numb(5,5),numf(5,5)

character*80 results

character*80 argv,intwrite

call getarg(11,results)

open(unit=7,file=results,access="append",status="unknown")

call getarg(12,results)

write(6,9900) results
9900 format("The results from this run have been appended to",
x       " the file `",a80,"`To see how to download",
x       " this file,",
x       ' read this description of anonymous ftp',
x       '.','Alternatively, you should be able to do a',
x       ' "Save As" or "Print" in your browser.')

write(7,9700)
9700 format(////,1x,"**************************************",//,1x,
x               "Output from astm1097.f",//,1x,
x               "**************************************",/)

c
c   reductions in properties
c

call getarg(1,argv)
write(intwrite,7) argv
7    format(a80)

call getarg(2,argv)
write(intwrite,7) argv

c  av(j) and bv(j) must be calculated from the fact that the lines go
c  through point (ln(2160 = 24x90 hrs),1.0 - red(j)) and point
c  (ln(.1 hr), 1.0)

testlen = log(2160.0)
tdiff = log(2160.0) - log(.1)
xbegin = log(.1)

do 8 i = 1,numred

red(i) = red(i)/100.0

bv(i) = -red(i)/tdiff
av(i) = 1.0 - bv(i)*xbegin

8    continue

c
c   cvs
c

call getarg(3,argv)
write(intwrite,7) argv

call getarg(4,argv)
write(intwrite,7) argv

do 55 i = 1,numcv

cv(i) = abs(cv(i)/100.0)

if ((cv(i) .eq. 0.0) .or. (cv(i) .ge. 1.0)) then

write(6,54)
54         format("A coefficient of variation value",
x      " did not lie between 0 and 100.  Please correct",
x      " the form.")

stop

endif

55   continue

c
c   sample sizes
c

call getarg(5,argv)
write(intwrite,7) argv

write(6,61) num1
write(7,62) num1
61   format('There were ',i3,' unaged observations',
x       '.')
62   format(/,1x,'There were',i3,' unaged observations',
x       '.',/)

call getarg(6,argv)
write(intwrite,7) argv

write(6,63) num2
write(7,64) num2
63   format('There were ',i3,' aged observations',
x       '.')
64   format(/,1x,'There were',i3,' aged observations',
x       '.',/)

call getarg(7,argv)
write(intwrite,7) argv

write(6,65) idos
write(7,66) idos
65   format('The ID of the order statistic',
x' from the unaged sample that was used to',
x' determine the load was ',i3,'.')
66   format(/,1x,'The ID of the order statistic',
x' from the unaged sample that was used',/,1x,
x'to determine the load was ',i3,'.',/)

call getarg(8,argv)
write(intwrite,7) argv

write(6,67) numbok
write(7,68) numbok
67   format('The number of failures',
x' permitted',
x' in the second sample was ',i3,'.')
68   format(/,1x,'The number of failures',
x' permitted in the second',
x' sample was ',i3,'.',/)

call getarg(9,argv)
write(intwrite,7) argv

write(6,69) numtrial
write(7,70) numtrial
69   format('The number of simulation',
x' trials requested was ',i6,'.')
70   format(/,1x,'The number of simulation',
x' trials requested was ',i6,'.',/)

call getarg(10,argv)
write(intwrite,7) argv

write(6,71) istart
write(7,72) istart
71   format('The start value for the random',
x' number generator was ',i8,'.')
72   format(/,1x,'The start value for the random',
x' number generator was ',i8,'.',/)

c
c   perform the simulation
c

do 110 i = 1,5

do 100 j = 1,5

numf(i,j) = 0
probf(i,j) = 0.0

100     continue

110  continue

dum = rnor(istart)

do 1000 iii = 1,numtrial

c
c   generate the first sample
c

do 120 j = 1,num1

samp1(j) = rnor(0)

120     continue

call rsort(samp1,sort1,iorigord,num1)

c
c   generate the second sample
c

do 130 j = 1,num2

samp2(j) = rnor(0)

130     continue

call rsort(samp2,sort2,iorigord,num2)

do 137 i = 1,numcv

do 135 j = 1,numred

numb(i,j) = 0

135        continue

137     continue

do 900 i = 1,numcv

do 140 j = 1,num2

140        continue

do 800 j = 1,numred

a = av(j)
b = bv(j)

do 700 k = 1,num2

xmu = (strrat(k) - av(j))/bv(j)
c
c  The standard deviation of the ln(time to failure) distribution
c  must be based on data --- variability of lifetimes
c  given a fixed stress ratio.  From the plot of the Wood data on
c  page 3 of Chuck Gerhards' "Effect of Duration and Rate of Loading
c  on Strength of Wood and Wood-based Materials", I obtained an estimate
c  of sigma (from ranges) as 10.4 mm where a log base 10 (hours) unit is 15.2 mm.
c  Thus the sd in log base 10 (hours) units was 10.4/15.2 = .684 which leads
c  to a sd of 1.57 in log base e (hours) units.
c

xlnt = xmu + 1.57*rnor(0)

if (xlnt .lt. testlen) then

numb(i,j) = numb(i,j) + 1

endif

700           continue

if (numb(i,j) .gt. numbok) numf(i,j) =
x                                    numf(i,j) + 1

800        continue

900     continue

1000 continue

xnumt = numtrial

do 1007 i = 1,numcv

do 1005 j = 1,numred

probf(i,j) = numf(i,j)/xnumt

1005    continue

1007 continue

write(6,1010)
write(7,1012) (red(i), i = 1,numred)
1010 format('Normal',
x'Nonparametric Approach',
x'Reductions are on top.  Coefficients',
x' of variation are along the left side.  Probabilities',
x' of failing the test are in the table.')
1012 format(///,1x,'Normal',//,1x,'Nonparametric Approach',
x//,1x,'Reductions are on top.  Coefficients of',
x' variation',/,1x,'are along the left side.  ',
x'Probabilities of failing the test',/,1x,'are in the table.',
x//,
x5x,5(3x,f4.2),/)

do 1016 i = 1,numcv
write(7,1014) cv(i),(probf(i,j),j = 1,numred)
1014    format(1x,f4.2,5(2x,f5.3))
1016 continue

write(6,1020)
1020 format("")

write(6,1030) numred
1030 format("",
x"")
do 1034 i = 1,numred
write(6,1032) red(i)
1032    format("")
1034 continue
write(6,1036)
1036 format("")

do 1050 i = 1,numcv

write(6,1042) cv(i),(probf(i,j), j = 1,numred)
1042    format("",(""))

write(6,1044)
1044    format("")

1050 continue

write(6,1052)
1052 format("Coefficientofvariation Reduction  ",f4.2,"",f5.3,"")

c
c   The lognormal case
c

lnp6 = log(.6)

do 3020 i = 1,numcv

sigma(i) = sqrt(log(cv(i)*cv(i) + 1.0))

3020 continue

c
c   perform the simulation
c

do 3110 i = 1,5

do 3100 j = 1,5

numf(i,j) = 0
probf(i,j) = 0.0

3100    continue

3110 continue

do 4000 iii = 1,numtrial

c
c   generate the first sample
c

do 3120 j = 1,num1

samp1(j) = rnor(0)

3120    continue

call rsort(samp1,sort1,iorigord,num1)

c
c   generate the second sample
c

do 3130 j = 1,num2

samp2(j) = rnor(0)

3130    continue

call rsort(samp2,sort2,iorigord,num2)

do 3137 i = 1,numcv

do 3135 j = 1,numred

numb(i,j) = 0

3135       continue

3137    continue

do 3900 i = 1,numcv

do 3140 j = 1,num2

3140       continue

do 3800 j = 1,numred

a = av(j)
b = bv(j)

do 3700 k = 1,num2

xmu = (strrat(k) - av(j))/bv(j)
c
c  The standard deviation of the ln(time to failure) distribution
c  must be based on data --- variability of lifetimes
c  given a fixed stress ratio.  From the plot of the Wood data on
c  page 3 of Chuck Gerhards' "Effect of Duration and Rate of Loading
c  on Strength of Wood and Wood-based Materials", I obtained an estimate
c  of sigma (from ranges) as 10.4 mm where a log base 10 (hours) unit is 15.2 mm.
c  Thus the sd in log base 10 (hours) units was 10.4/15.2 = .684 which leads
c  to a sd of 1.57 in log base e (hours) units.
c

xlnt = xmu + 1.57*rnor(0)

if (xlnt .lt. testlen) then

numb(i,j) = numb(i,j) + 1

endif

3700          continue

if (numb(i,j) .gt. numbok) numf(i,j) =
x                                    numf(i,j) + 1

3800       continue

3900    continue

4000 continue

xnumt = numtrial

do 4007 i = 1,numcv

do 4005 j = 1,numred

probf(i,j) = numf(i,j)/xnumt

4005    continue

4007 continue

write(6,4010)
write(7,4012) (red(i), i = 1,numred)
4010 format('Lognormal',
x'Nonparametric Approach',
x'Reductions are on top.  Coefficients',
x' of variation are along the left side.  Probabilities',
x' of failing the test are in the table.')
4012 format(///,1x,'Lognormal',//,1x,'Nonparametric Approach',
x//,1x,'Reductions are on top.  Coefficients of',
x' variation',/,1x,'are along the left side.  ',
x'Probabilities of failing the test',/,1x,'are in the table.',
x//,
x5x,5(3x,f4.2),/)

do 4016 i = 1,numcv
write(7,4014) cv(i),(probf(i,j),j = 1,numred)
4014    format(1x,f4.2,5(2x,f5.3))
4016 continue

write(6,4020)
4020 format("")

write(6,4030) numred
4030 format("",
x"")
do 4034 i = 1,numred
write(6,4032) red(i)
4032    format("")
4034 continue
write(6,4036)
4036 format("")

do 4050 i = 1,numcv

write(6,4042) cv(i),(probf(i,j), j = 1,numred)
4042    format("",(""))

write(6,4044)
4044    format("")

4050 continue

write(6,4052)
4052 format("Coefficientofvariation Reduction  ",f4.2,"",f5.3,"")

write(6,6060)
6060 format('For further information,',
x'',
x'sverrill@fs.fed.us',

stop

end

REAL FUNCTION RNOR(JD)
C***BEGIN PROLOGUE  RNOR
C***DATE WRITTEN   810915
C***REVISION DATE  830805
C***CATEGORY NO.  L6A14
C***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
C***AUTHOR    KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
C             MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV
C
C***PURPOSE  GENERATES QUASI NORMAL RANDOM NUMBERS, WITH MEAN ZERO AND
C             UNIT STANDARD DEVIATION, AND CAN BE USED WITH ANY COMPUTER
C             WITH INTEGERS AT LEAST AS LARGE AS 32767.
C***DESCRIPTION
C
C       RNOR generates quasi normal random numbers with zero mean and
C       unit standard deviation.
C       It can be used with any computer with integers at least as
C       large as 32767.
C
C
C   Use
C       First time....
C                   Z = RNOR(JD)
C                     Here JD is any  n o n - z e r o  integer.
C                     This causes initialization of the program
C                     and the first random number to be returned as Z.
C       Subsequent times...
C                   Z = RNOR(0)
C                     Causes the next random number to be returned as Z.
C
C.....................................................................
C
C    Note: Users who wish to transport this program to other
C           computers should read the following ....
C
C   Machine dependencies...
C      MDIG = A lower bound on the number of binary digits available
C              for representing integers, including the sign bit.
C              This must be at least 16, but can be increased in
C              line with remark A below.
C
C   Remarks...
C     A. This program can be used in two ways:
C        (1) To obtain repeatable results on different computers,
C            set 'MDIG' to the smallest of its values on each, or,
C        (2) To allow the longest sequence of random numbers to be
C            generated without cycling (repeating) set 'MDIG' to the
C            largest possible value.
C     B. The sequence of numbers generated depends on the initial
C          input 'JD' as well as the value of 'MDIG'.
C          If MDIG=16 one should find that
C            the first evaluation
C              Z=RNOR(87) gives  Z=-.40079207...
C            The second evaluation
C              Z=RNOR(0) gives   Z=-1.8728870...
C            The third evaluation
C              Z=RNOR(0) gives   Z=1.8216004...
C            The fourth evaluation
C              Z=RNOR(0) gives   Z=.69410355...
C            The thousandth evaluation
C              Z=RNOR(0) gives   Z=.96782424...
C
C***REFERENCES  MARSAGLIA & TSANG, "A FAST, EASILY IMPLEMENTED
C                 METHOD FOR SAMPLING FROM DECREASING OR
C                 SYMMETRIC UNIMODAL DENSITY FUNCTIONS", TO BE
C                 PUBLISHED IN SIAM J SISC 1983.
C***ROUTINES CALLED  I1MACH,XERROR
C***END PROLOGUE  RNOR
REAL V(65),W(65)
INTEGER M(17)
SAVE I1,J1,M,M1,M2,RMAX
DATA AA,B,C,RMAX/12.37586,.4878992,12.67706,3.0518509E-5/
DATA C1,C2,PC,XN/.9689279,1.301198,.1958303E-1,2.776994/
DATA V/ .3409450, .4573146, .5397793, .6062427, .6631691
+, .7136975, .7596125, .8020356, .8417227, .8792102, .9148948
+, .9490791, .9820005, 1.0138492, 1.0447810, 1.0749254, 1.1043917
+,1.1332738, 1.1616530, 1.1896010, 1.2171815, 1.2444516, 1.2714635
+,1.2982650, 1.3249008, 1.3514125, 1.3778399, 1.4042211, 1.4305929
+,1.4569915, 1.4834526, 1.5100121, 1.5367061, 1.5635712, 1.5906454
+,1.6179680, 1.6455802, 1.6735255, 1.7018503, 1.7306045, 1.7598422
+,1.7896223, 1.8200099, 1.8510770, 1.8829044, 1.9155830, 1.9492166
+,1.9839239, 2.0198430, 2.0571356, 2.0959930, 2.1366450, 2.1793713
+,2.2245175, 2.2725185, 2.3239338, 2.3795007, 2.4402218, 2.5075117
+,2.5834658, 2.6713916, 2.7769943, 2.7769943, 2.7769943, 2.7769943/
DATA W/   .10405134E-04, .13956560E-04, .16473259E-04,
+ .18501623E-04, .20238931E-04, .21780983E-04, .23182241E-04,
+ .24476931E-04, .25688121E-04, .26832186E-04, .27921226E-04,
+ .28964480E-04, .29969191E-04, .30941168E-04, .31885160E-04,
+ .32805121E-04, .33704388E-04, .34585827E-04, .35451919E-04,
+ .36304851E-04, .37146564E-04, .37978808E-04, .38803170E-04,
+ .39621114E-04, .40433997E-04, .41243096E-04, .42049621E-04,
+ .42854734E-04, .43659562E-04, .44465208E-04, .45272764E-04,
+ .46083321E-04, .46897980E-04, .47717864E-04, .48544128E-04,
+ .49377973E-04, .50220656E-04, .51073504E-04, .51937936E-04,
+ .52815471E-04, .53707761E-04, .54616606E-04, .55543990E-04,
+ .56492112E-04, .57463436E-04, .58460740E-04, .59487185E-04,
+ .60546402E-04, .61642600E-04, .62780711E-04, .63966581E-04,
+ .65207221E-04, .66511165E-04, .67888959E-04, .69353880E-04,
+ .70922996E-04, .72618816E-04, .74471933E-04, .76525519E-04,
+ .78843526E-04, .81526890E-04, .84749727E-04,
+ .84749727E-04, .84749727E-04, .84749727E-04/
DATA M(1),M(2),M(3),M(4),M(5),M(6),M(7),M(8),M(9),M(10),M(11),
1     M(12),M(13),M(14),M(15),M(16),M(17)
2                   / 30788,23052,2053,19346,10646,19427,23975,
3                     19049,10949,19693,29746,26748,2796,23890,
4                     29168,31924,16499 /
DATA M1,M2,I1,J1 / 32767,256,5,17 /
C Fast part...
C
C
C***FIRST EXECUTABLE STATEMENT  RNOR
IF(JD.NE.0)GO TO 27
10 CONTINUE
I=M(I1)-M(J1)
IF(I .LT. 0) I=I+M1
M(J1)=I
I1=I1-1
IF(I1 .EQ. 0) I1=17
J1=J1-1
IF(J1 .EQ. 0) J1=17
J=MOD(I,64)+1
RNOR=I*W(J+1)
IF( ( (I/M2)/2 )*2.EQ.(I/M2))RNOR=-RNOR
IF(ABS(RNOR).LE.V(J))RETURN
C Slow part; AA is a*f(0)
X=(ABS(RNOR)-V(J))/(V(J+1)-V(J))
Y=UNI(0)
S=X+Y
IF(S.GT.C2)GO TO 11
IF(S.LE.C1)RETURN
IF(Y.GT.C-AA*EXP(-.5*(B-B*X)**2))GO TO 11
IF(EXP(-.5*V(J+1)**2)+Y*PC/V(J+1).LE.EXP(-.5*RNOR**2))RETURN
C Tail part; 3.855849 is .5*XN**2
c
c   the three lines after the next empty comment line were added by
c   Steve Verrill on 10/23/89
c
22    tttttt = uni(0)
if (tttttt .eq. 0.0) go to 22
s = xn - alog(tttttt)/xn
c   22 S=XN-ALOG(UNI(0))/XN
c
c   the three lines after the next empty comment line were
c   added by Steve Verrill on 10/23/89
c
23    tttttt = uni(0)
if (tttttt .eq. 0.0) go to 23
if (3.855849 + alog(tttttt) - xn*s .gt. -.5*s**2) go to 22
c      IF(3.855849+ALOG(UNI(0))-XN*S.GT.-.5*S**2)GO TO 22
RNOR=SIGN(S,RNOR)
RETURN
11 RNOR=SIGN(B-B*X,RNOR)
RETURN
C  FILL
27 CONTINUE
MDIG=32
C          BE SURE THAT MDIG AT LEAST 16...
ccc      IF(MDIG.LT.16)CALL XERROR('RNOR--MDIG LESS THAN 16',23,1,2)
M1 = 2**(MDIG-2) + (2**(MDIG-2)-1)
M2 = 2**(MDIG/2)
JSEED = MIN0(IABS(JD),M1)
IF( MOD(JSEED,2).EQ.0 ) JSEED=JSEED-1
K0 =MOD(9069,M2)
K1 = 9069/M2
J0 = MOD(JSEED,M2)
J1 = JSEED/M2
DO 2 I=1,17
JSEED = J0*K0
J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2)
J0 = MOD(JSEED,M2)
2   M(I) = J0+M2*J1
J1=17
I1=5
RMAX = 1./FLOAT(M1)
C        Seed uniform (0,1) generator.  (Just a dummy call)
RNOR=UNI(JD)
DO 28 I=1,65
28  W(I)=RMAX*V(I)
GO TO 10
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
REAL FUNCTION UNI(JD)
C***BEGIN PROLOGUE  UNI
C***DATE WRITTEN   810915
C***REVISION DATE  830805
C***CATEGORY NO.  L6A21
C***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
C***AUTHOR    BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS
C             KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
C             MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV
C
C***PURPOSE  THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1)
C             AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS
C             AT LEAST AS LARGE AS 32767.
C***DESCRIPTION
C
C       THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTERVAL
C       [0,1).  IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS
C       INTEGERS AT LEAST AS LARGE AS 32767.
C
C
C   USE
C       FIRST TIME....
C                   Z = UNI(JD)
C                     HERE JD IS ANY  N O N - Z E R O  INTEGER.
C                     THIS CAUSES INITIALIZATION OF THE PROGRAM
C                     AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z.
C       SUBSEQUENT TIMES...
C                   Z = UNI(0)
C                     CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z.
C
C
C..................................................................
C   NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER
C         TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION.....
C
C   MACHINE DEPENDENCIES...
C      MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE
C              FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT.
C              THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED
C              IN LINE WITH REMARK A BELOW.
C
C   REMARKS...
C     A. THIS PROGRAM CAN BE USED IN TWO WAYS:
C        (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS,
C            SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR,
C        (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE
C            GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE
C            LARGEST POSSIBLE VALUE.
C     B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL
C          INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'.
C          IF MDIG=16 ONE SHOULD FIND THAT
C            THE FIRST EVALUATION
C              Z=UNI(305) GIVES Z=.027832881...
C            THE SECOND EVALUATION
C              Z=UNI(0) GIVES   Z=.56102176...
C            THE THIRD EVALUATION
C              Z=UNI(0) GIVES   Z=.41456343...
C            THE THOUSANDTH EVALUATION
C              Z=UNI(0) GIVES   Z=.19797357...
C
C***REFERENCES  MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM
C                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
C***ROUTINES CALLED  I1MACH,XERROR
C***END PROLOGUE  UNI
INTEGER M(17)
C
SAVE I,J,M,M1,M2
C
DATA M(1),M(2),M(3),M(4),M(5),M(6),M(7),M(8),M(9),M(10),M(11),
1     M(12),M(13),M(14),M(15),M(16),M(17)
2                   / 30788,23052,2053,19346,10646,19427,23975,
3                     19049,10949,19693,29746,26748,2796,23890,
4                     29168,31924,16499 /
DATA M1,M2,I,J / 32767,256,5,17 /
C***FIRST EXECUTABLE STATEMENT  UNI
IF(JD .EQ. 0) GO TO 3
C  FILL
MDIG=32
C          BE SURE THAT MDIG AT LEAST 16...
c      IF(MDIG.LT.16)CALL XERROR('UNI--MDIG LESS THAN 16',22,1,2)
M1= 2**(MDIG-2) + (2**(MDIG-2)-1)
M2 = 2**(MDIG/2)
JSEED = MIN0(IABS(JD),M1)
IF( MOD(JSEED,2).EQ.0 ) JSEED=JSEED-1
K0 =MOD(9069,M2)
K1 = 9069/M2
J0 = MOD(JSEED,M2)
J1 = JSEED/M2
DO 2 I=1,17
JSEED = J0*K0
J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2)
J0 = MOD(JSEED,M2)
2   M(I) = J0+M2*J1
I=5
J=17
C  BEGIN MAIN LOOP HERE
3 K=M(I)-M(J)
IF(K .LT. 0) K=K+M1
M(J)=K
I=I-1
IF(I .EQ. 0) I=17
J=J-1
IF(J .EQ. 0) J=17
UNI=FLOAT(K)/FLOAT(M1)
RETURN
END

c   This subroutine attempts to perform an n*log(n) sort
c   rather than an n*n sort.
c   It sorts reals.
c   It also keeps track of the original order of the x's.
c
c
c   Steve Verrill, 5/2/89
c
c
c
subroutine rsort(x,xsort,iorigord,n)
c
dimension x(n),xsort(n),iorigord(n)
c
do 10 i = 1,n
iorigord(i) = i
10    continue
c
num = 2
c
if (x(1) .lt. x(2)) then
xsort(1) = x(1)
iorigord(1) = 1
xsort(2) = x(2)
iorigord(2) = 2
else
xsort(1) = x(2)
iorigord(1) = 2
xsort(2) = x(1)
iorigord(2) = 1
endif
c
do 1000 j = 3,n
c
c   insert the next value
c
if (x(j) .ge. xsort(num)) go to 200
if (x(j) .le. xsort(1)) go to 160
il = 1
iu = num
145   numb = iu - il
if (numb .eq. 1) go to 170
numbd2 = numb/2
itest = il + numbd2
if (x(j) .le. xsort(itest)) go to 150
il = itest
go to 145
150   iu = itest
go to 145
160   il = 0
170   num = num + 1
do 180 k = num,il + 2,-1
km1 = k - 1
xsort(k) = xsort(km1)
iorigord(k) = iorigord(km1)
180   continue
190   ilp1 = il + 1
xsort(ilp1) = x(j)
iorigord(ilp1) = j
go to 1000
c
200   num = num + 1
xsort(num) = x(j)
iorigord(num) = j
c
1000  continue
c
return
end
```