REM written sept 7, 1997, gerald L Swanson
REM
REM     compute solid elipsiod  coordinates
REM
REM     X * X   Y * Y   Z * Z
REM     ----- + ----- + -----  <= 1
REM     A * A   B * B   C * C
REM
rem output is layer by layer plan consisting of
rem	row	(0-A)
rem	start column (0-b)
rem	end column (0-b)
rem	asterix's.
rem
rem	level 1
rem	0 1 4  ****
rem     1 0 3 ****
rem     2 0 2 ***
rem     3 0 0 *


REM coordinates
DIM X, Y, Z AS SINGLE

REM axes length. for a sphere, A = B = 6*C/5
DIM A, B, C AS SINGLE

REM polinomials
DIM xx, yy, zz AS SINGLE
DIM aa, bb, cc AS SINGLE

REM beginning and ending for x axis
DIM xb, xe AS SINGLE

REM current and next layer
DIM cz, nz  AS SINGLE

REM Thickness of walls
DIM S, SS AS SINGLE

REM Start of each string of bricks
DIM mm AS INTEGER

INPUT A
aa = A * A
INPUT B
bb = B * B
INPUT C
cc = C * C

INPUT S
SS = 2 * S

REM end of each string of bricks for the current layer and S layers above
REM and below.
DIM M(A + S, SS) AS INTEGER



REM initialize
cz = 0
nz = S
FOR Z = 0 TO SS
        FOR X = 0 TO A + S
                M(X, Z) = -1
NEXT X, Z

REM generate each layer
FOR Z = C + S TO -S STEP -1
        zz = (Z * Z) / cc
        PRINT "level ", Z + S

REM compute end of each row of the next layer

        FOR X = 0 TO A + S
                xx = (X * X) / aa
                yy = (1! - xx - zz) * bb
                IF yy >= 0 THEN
                        M(X, nz) = FIX(SQR(yy))
                ELSE
                        M(X, nz) = -1
                END IF
        NEXT X

REM get beginning of each row  of the current layer
REM and output
      
        FOR X = 0 TO A + S
                IF M(X, cz) > -1 THEN
                        mm = M(X, cz) - S
                        xb = X - 1
                        IF xb < 0 THEN xb = 0
                        xe = X + 1
                        IF xe > A THEN xe = A
                        FOR xx = xb TO xe
                                FOR zz = 0 TO SS
                                        IF mm > M(xx, zz) - S THEN mm = M(xx, zz) - S
                                NEXT zz
                        NEXT xx
                        IF mm < 0 THEN mm = 0
                        PRINT USING "## #### #### &"; X; mm; M(X, cz); STRING$(mm + 1, " ") + STRING$(M(X, cz) - mm + 1, "*")
                END IF
        NEXT X
        INPUT X

REM next layer
cz = cz + 1
IF cz > SS THEN cz = 0
nz = nz + 1
IF nz > SS THEN nz = 0
NEXT Z