'This program produces density plots of hydrogen atom wavefunctions.
'It produces plots for the |100>, |200>, and |210> states.
'It produces projections of 3D wavefuctions.
'To produce cuts set y=0.

DIM a1(60, 60), a2(60, 60), a3(60, 60), psi(3), col(10)
'Here we choose colors to represent different densities.
col(1) = 8
col(2) = 1
col(3) = 9
col(4) = 5
col(5) = 13
col(6) = 4
col(7) = 12
col(8) = 6
col(9) = 14
col(10) = 15

'We choose a VGA screen.
SCREEN 12
WINDOW (0, 0)-(64, 48)
CLS

RANDOMIZE TIMER
begin:
'We pick a random point in space.
x = -10 + 20 * RND(1)
y = -10 + 20 * RND(1)
'Set y=0 if you want a cut instead of a projection.
z = -10 + 20 * RND(1)
r = SQR(x * x + y * y + z * z)

'We evaluate the square of the wavefunctions at this point.
'The wavefunctions are normalized so that the maximum value
'of the square is one.
IF r < 10 THEN
        psi(1) = EXP(-2 * r)
        psi(2) = 13.6 * (2 - r) ^ 2 * EXP(-r)
        psi(3) = 2.7 * z ^ 2 * EXP(-r)
       
'We pick a random number.  If the random number is smaller
'than the square of the wavefunction (proportional to the probability 
'of finding the electron at our point) we say that we found an electron
'at this point.  We increment an array element which holds the number of
'electrons detected at this point.  If the random number is larger, we do
'not increment. We then pick the color that corresponds to the number in the
'array element and plot a point.
'The |100> state:
        IF RND(1) < psi(1) THEN
                x1 = (x + 10) * 3
                z1 = (z + 10) * 3
                a1(x1, z1) = a1(x1, z1) + 1
                pt = SQR(a1(x1, z1))
                IF pt > 10 THEN GOTO fini
                CIRCLE (x + 16, z + 36), .1, col(pt)
        END IF
'The |200> state:
        IF RND(1) < psi(2) THEN
        x1 = (x + 10) * 3
        z1 = (z + 10) * 3
        a2(x1, z1) = a2(x1, z1) + 1
        pt = SQR(a2(x1, z1))
        IF pt > 10 THEN GOTO fini
        CIRCLE (x + 48, z + 36), .1, col(pt)
        END IF
'The |210> state
        IF RND(1) < psi(3) THEN
        x1 = (x + 10) * 3
        z1 = (z + 10) * 3
        a3(x1, z1) = a3(x1, z1) + 1
        pt = SQR(a3(x1, z1))
        IF pt > 10 THEN GOTO fini
        CIRCLE (x + 32, z + 12), .1, col(pt)
        END IF

END IF
GOTO begin
fini:
END