REM *** Increasing differences model: Solar System
REM *** Use broken stick technique to generate distances,
REM *** then sort the distances to get random system with
REM *** increasing distances. Then figure out new A(I)'s.
REM *** Then calculate R-squared and add to histogram.

REM *** input the number of orbitals (planets)
LET ORBITALS = 9

REM *** input the semimajor axis of the innermost planet
REM *** note that INNERMOST = mercury's semimajor axis
REM *** less A(min) (which is 0.15 AU)
LET INNERMOST = 0.23710

REM *** input the semimajor axis of the outermost planet
LET OUTERMOST = 29.91107

REM *** input average x for least squares now, because it
REM *** never changes
LET AVGX = 5

REM *** number of trials between PRNG seed changing
LET MAXTRIALS = 100000

REM *** number of PRNG seed changes
REM *** MAXTRIALS x MAXCYCLES = total # systems generated
LET MAXCYCLES = 10000

REM *** A(i) is for the semimajor axes. A(1) and A(outermost)
REM *** never change, so we can set their values now
DIM A(ORBITALS)
LET A(1) = INNERMOST
LET A(ORBITALS) = OUTERMOST

REM *** we transform the semimajor axes into logs to compute
REM *** the R-squared goodness of fit statistic
DIM LOGA(ORBITALS)
LET LOGA(1) = LOG(INNERMOST)
LET LOGA(ORBITALS) = LOG(OUTERMOST)

REM *** D(i) the differences between consecutive orbitals
DIM D(ORBITALS - 1)

REM *** X(i) is for the x-axis in the least squares zone
REM *** X(i) aka the ordinal of the orbitals w/ planets

DIM X(ORBITALS)

REM *** load the x-axis that keeps track of 9 planets
REM *** note it goes {1,2,3,4,5,6,7,8,9} because the
REM *** scenario calls for 9 planets, no skipped orbitals,
FOR I = 1 TO ORBITALS
LET X(I) = I
NEXT I

REM *** HIST(i) is the histogram that defines the
REM *** probability distribution which is the final goal
DIM HIST(10001)

FOR CYCLE = 1 TO MAXCYCLES
REM *** we print the cycle # to keep track of progress of program
PRINT CYCLE
RANDOMIZE
FOR TRIAL = 1 TO MAXTRIALS
REM *** We are normalizing to the size of the whole system.
REM *** We generate random numbers from within the innermost
REM *** to the outermost orbitals
FOR I = 2 TO (ORBITALS - 1)
LET A(I) = INNERMOST + (OUTERMOST - INNERMOST) * RND
NEXT I

REM *** now sort in ascending order
LET J = 3
DO
LET I = J
DO
IF A(I - 1) > A(I) THEN
LET BUFF = A(I)
LET A(I) = A(I - 1)
LET A(I - 1) = BUFF
END IF
LET I = I - 1
LOOP UNTIL I = 2
LET J = J + 1
LOOP UNTIL J = ORBITALS

REM *** calculate the length of the sticks
FOR I = 1 TO ORBITALS - 1
LET D(I) = A(I+1) - A(I)
NEXT I

REM *** now sort the sticks
LET J = 2
DO
LET I = J
DO
IF D(I - 1) > D(I) THEN
LET BUFF = D(I)
LET D(I) = D(I - 1)
LET D(I - 1) = BUFF
END IF
LET I = I - 1
LOOP UNTIL I = 1
LET J = J + 1
LOOP UNTIL J = ORBITALS

REM *** now reconstruct the semimajor axes
REM *** you already know A(1) and A(ORBITALS)
FOR I = 2 TO ORBITALS - 1
LET A(I) = A(I - 1) + D(I - 1)
NEXT I

REM *** the following statements are for debugging purposes
REM FOR I = 1 TO ORBITALS
REM PRINT A(I)
REM NEXT I
REM PRINT "R-SQUARED = ";

REM *** time to figure out the r-squared score
REM *** start by tranforming data by taking logarithm
REM *** of the semimajor axes, A(i)
REM *** we already know the logs of the first and last A(i)
LET SUMY = LOGA(1) + LOGA(ORBITALS)
FOR I = 2 TO ORBITALS - 1
LET LOGA(I) = LOG(A(I))
LET SUMY = SUMY + LOGA(I)
NEXT I
LET AVGY = SUMY / ORBITALS

REM *** now figure the sums of the squares of the
REM *** deviations (x - avg(x)) and (y - avg(y))
REM *** note I don't do X^2 because X*X is faster
LET SUMXSQ = 0
LET SUMYSQ = 0
LET SUMXY = 0
FOR I = 1 TO ORBITALS - 1
LET SUMXSQ = SUMXSQ + ((X(I) - AVGX) * (X(I) - AVGX))
LET SUMYSQ = SUMYSQ + ((LOGA(I) - AVGY) * (LOGA(I) - AVGY))
LET SUMXY = SUMXY + ((X(I) - AVGX) * (LOGA(I) - AVGY))
NEXT I

REM *** B is the slope of the best-fit least squared
REM *** line of the form y = a + bx
LET B = SUMXY / SUMXSQ

REM *** SSRESID is the so-called residual sum of squares
LET SSRESID = SUMYSQ - (B * SUMXY)

REM *** now we calculate the R-squared score!
LET RSQ = 1 - (SSRESID / SUMYSQ)

REM *** next statement for debugging purposes
REM PRINT RSQ

REM *** this adds the R-squared to the histogram that will
REM *** eventually form the probability distribution for
REM *** determining the critical value that the observed
REM *** R-squared of 55Cnc must exceed in order to reject
REM *** the law of increasing differences null hypothesis
LET HIST(INT(RSQ * 10000)) = HIST(INT(RSQ * 10000)) + 1
NEXT TRIAL
NEXT CYCLE

REM *** time to output the probability distribution
REM *** the left hand tail is of no interest, so we skip that part
FOR I = 1 TO 4999
LET COUNTER = COUNTER + HIST(I)
NEXT I
FOR I = 5000 TO 9999
REM *** prints out the bins corresponding an R-squared score
PRINT USING "%.####": (I / 10000);

REM *** print out the number of systems for that bin
PRINT USING "##########": HIST(I);

REM *** the running total of bins from lowest to highest
LET COUNTER = COUNTER + HIST(I)
PRINT USING "###########": COUNTER
NEXT I

LET TOTALTRIALS = (TRIAL -1) * (CYCLE - 1)
PRINT "TOTAL TRIALS = "; TOTALTRIALS
PRINT "FOR TOTAL PLANETS = "; ORBITALS

PRINT "TOTAL LENGTH IS "; OUTERMOST - INNERMOST; "FROM "; INNERMOST; "TO "; OUTERMOST
END

__________________
"Ignorance more frequently begets confidence than does knowledge" -- Charles Darwin