DECLARE SUB getXY (x0!, y0!, z0!, sx%, sy%)
DECLARE SUB mRotate (an!, plane%)
DECLARE SUB mMult ()
DECLARE SUB mIdentity ()
declare sub mScale (factor!)
declare SUB rotpnt (x0!, y0!, z0!, xr!, yr!, zr!, nx!, ny!, nz!)

DEFINT A-Z
randomize timer

const xplane = 0, yplane = 1, zplane = 2

DIM SHARED apage, bpage, zback, rad as single
DIM SHARED matrix(8) AS SINGLE, tempmatrix(2, 8) AS SINGLE, pnti(4,1)

CONST MAX_NODES = 13, MAX_BRANCHES = 13, MAX_LEAVES = 5, MAX_SEGS = 3

type pointType
    x   as single
    y   as single
    z   as single
end type

TYPE nodeType
  x AS SINGLE
  y AS SINGLE
  z AS SINGLE
  globalNode    as integer
  plane(MAX_SEGS - 1) as pointType
END TYPE

TYPE nodeAngleType
  x AS SINGLE
  y AS SINGLE
  z AS SINGLE
end type

TYPE branchType
  nNodes  AS INTEGER
  nodes(MAX_NODES) as nodeType
  nodeAngle(MAX_NODES) as nodeAngleType
  root  as nodeType
END TYPE

DIM branches(MAX_BRANCHES) AS branchType


screen 13,,2

apage = 0
bpage = 1
rad = 3.141592 / 180
zback = 30
matrixScale! = .8
xan! = -20
rotSpd! = .05
scaleSpd! = .001


sag! = 10
growth! = 0
scale! = .5

top:
t# = timer + 1

twistyness = 30
nleaves = 0
nBranches = MAX_BRANCHES
totalNodes = MAX_NODES
baseBranches = 1

trunkLength = 4
twigLength = 1.5
branchLengthInter! = (twigLength - trunkLength) / totalNodes


trunkWidth = 1.4
twigWidth = 0
branchWidthInter! = (twigWidth - trunkWidth) / totalNodes

for branch = 0 to nBranches - 1
    for i = 0 to totalNodes - 1
        branches(branch).nodes(i).globalNode = 0
        branches(branch).nodeAngle(i).x = 0
        branches(branch).nodeAngle(i).y = 0
        branches(branch).nodeAngle(i).z = 0
        branches(branch).nodes(i).x = 0
        branches(branch).nodes(i).y = 0
        branches(branch).nodes(i).z = 0
    next
    do
        parentBranch = int(rnd * branch)
        if baseBranches > 0 then parentBranch = -1
        parentNode = -1
        if parentBranch > -1 then
            parentNode = int(rnd * branches(parentBranch).nNodes)
            branches(branch).root.x = branches(parentBranch).nodes(parentNode).x
            branches(branch).root.y = branches(parentBranch).nodes(parentNode).y
            branches(branch).root.z = branches(parentBranch).nodes(parentNode).z
            branches(branch).nodeAngle(0).x = branches(parentBranch).nodeAngle(parentNode).x
            branches(branch).nodeAngle(0).y = branches(parentBranch).nodeAngle(parentNode).y + rnd * 180 - 90
            branches(branch).nodeAngle(0).z = branches(parentBranch).nodeAngle(parentNode).z
            globalNode = branches(parentBranch).nodes(parentNode).globalNode
        else
            branches(branch).root.x = 0
            branches(branch).root.y = 0
            branches(branch).root.z = 0
            branches(branch).nodeAngle(0).x = 0
            branches(branch).nodeAngle(0).y = rnd * 360
            branches(branch).nodeAngle(0).z = 0
            globalNode = 0
            baseBranches = baseBranches - 1
        end if
        nNodes = totalNodes - globalNode
        branches(branch).nNodes = nNodes
        
        if nNodes > 0 and parentNode <> 0 then exit do
    loop
    
    ox! = branches(branch).root.x
    oy! = branches(branch).root.y
    oz! = branches(branch).root.z
    for node = 0 to nNodes - 1
        globalNode = globalNode + 1
        rx! = branches(branch).nodeAngle(node).x + rnd * twistyness - twistyness / 2
        ry! = branches(branch).nodeAngle(node).y + rnd * twistyness - twistyness / 2
        rz! = branches(branch).nodeAngle(node).z + rnd * twistyness - twistyness / 2
        rz! = rz! + rnd * sag! * rz! / 90 + sag! / 5
        rz! = rz! - rnd * growth! * rz! / 90 - growth! / 5
        rotpnt 0, (trunkLength + globalNode * branchLengthInter!) * scale!, 0, 0, 0, rz!, nx!, ny!, nz!
        rotpnt nx!, ny!, nz!, rx!, 0, 0, nx!, ny!, nz!
        rotpnt nx!, ny!, nz!, 0, ry!, 0, nx!, ny!, nz!
        branches(branch).nodes(node).x = ox! + nx!
        branches(branch).nodes(node).y = oy! + ny!
        branches(branch).nodes(node).z = oz! + nz!
        branches(branch).nodeAngle(node+1).x = rx!
        branches(branch).nodeAngle(node+1).y = ry!
        branches(branch).nodeAngle(node+1).z = rz!
        branches(branch).nodes(node).globalNode = globalNode
        ox! = branches(branch).nodes(node).x
        oy! = branches(branch).nodes(node).y
        oz! = branches(branch).nodes(node).z
        
        for p = 0 to MAX_SEGS - 1
            rotpnt 0, 0, (trunkWidth * cos(3.14159 / 2 * globalNode / totalNodes)) * scale!, 0, p * 360/MAX_SEGS, 0, nx!, ny!, nz!
            rotpnt nx!, ny!, nz!, 0, 0, rz!, nx!, ny!, nz!
            rotpnt nx!, ny!, nz!, rx!, 0, 0, nx!, ny!, nz!
            rotpnt nx!, ny!, nz!, 0, ry!, 0, nx!, ny!, nz!
        
            branches(branch).nodes(node).plane(p).x = ox! + nx!
            branches(branch).nodes(node).plane(p).y = oy! + ny!
            branches(branch).nodes(node).plane(p).z = oz! + nz!
        next
        
    next
next    

DO

    swap apage, bpage
    screenset apage, bpage
    cls

    mIdentity
    mRotate yan!, yplane
    mRotate xan!, xplane
    mScale matrixScale!
    
    FOR x! = -10 TO 10 STEP 2
        FOR z! = -10 TO 10 STEP 2
        
            getXY x!, 0, z!, sx, sy
            PSET (sx, sy), 7
          
        NEXT
    NEXT
    
    
    for branch = 0 to nBranches - 1
        nNodes = branches(branch).nNodes
        for node = 0 to nNodes - 1
            for p = 0 to MAX_SEGS - 1
                n = p + 1
                if n = MAX_SEGS then n = 0
                nx! = branches(branch).nodes(node).plane(n).x
                ny! = branches(branch).nodes(node).plane(n).y
                nz! = branches(branch).nodes(node).plane(n).z
                getXY nx!, ny!, nz!, pnti(0,0), pnti(0,1)
                n = p
                nx! = branches(branch).nodes(node).plane(n).x
                ny! = branches(branch).nodes(node).plane(n).y
                nz! = branches(branch).nodes(node).plane(n).z
                getXY nx!, ny!, nz!, pnti(3,0), pnti(3,1)
                
                nn = node - 1
                if nn > -1 then
                    n = p + 1
                    if n = MAX_SEGS then n = 0
                    nx! = branches(branch).nodes(nn).plane(n).x
                    ny! = branches(branch).nodes(nn).plane(n).y
                    nz! = branches(branch).nodes(nn).plane(n).z
                    getXY nx!, ny!, nz!, pnti(1,0), pnti(1,1)
                    n = p
                    nx! = branches(branch).nodes(nn).plane(n).x
                    ny! = branches(branch).nodes(nn).plane(n).y
                    nz! = branches(branch).nodes(nn).plane(n).z
                    getXY nx!, ny!, nz!, pnti(2,0), pnti(2,1)
                    refs = 4
                else
                    nx! = branches(branch).root.x
                    ny! = branches(branch).root.y
                    nz! = branches(branch).root.z
                    getXY nx!, ny!, nz!, pnti(1,0), pnti(1,1)
                    refs = 3
                end if
                
                c = 6
                if refs = 3 then
                    line (pnti(0,0), pnti(0,1)) - (pnti(1,0), pnti(1,1)), c
                    line (pnti(0,0), pnti(0,1)) - (pnti(3,0), pnti(3,1)), c
                    line (pnti(1,0), pnti(1,1)) - (pnti(3,0), pnti(3,1)), c
                else
                    line (pnti(0,0), pnti(0,1)) - (pnti(1,0), pnti(1,1)), c
                    line (pnti(2,0), pnti(2,1)) - (pnti(1,0), pnti(1,1)), c
                    line (pnti(2,0), pnti(2,1)) - (pnti(3,0), pnti(3,1)), c
                    line (pnti(0,0), pnti(0,1)) - (pnti(3,0), pnti(3,1)), c
                end if
            next
        next
    next
    
    if multikey(&h1) then exit do
    
    if multikey(&h48) then xan! = xan! + rotSpd!
    if multikey(&h50) then xan! = xan! - rotSpd!
    if multikey(&h4b) then yan! = yan! - rotSpd!
    if multikey(&h4d) then yan! = yan! + rotSpd!
    
    if multikey(&h4e) then matrixScale! = matrixScale! + scaleSpd!
    if multikey(&h4a) then matrixScale! = matrixScale! - scaleSpd!
    
    if t# < timer then
        if multikey(&h1c) then goto top
    end if
    
    while inkey$ <> "": wend

LOOP





SUB getXY (x0!, y0!, z0!, sx, sy)

  nx! = x0! * matrix(0) + y0! * matrix(3) + z0! * matrix(6)
  ny! = x0! * matrix(1) + y0! * matrix(4) + z0! * matrix(7)
  nz! = x0! * matrix(2) + y0! * matrix(5) + z0! * matrix(8) - zback

  IF nz! = 0 THEN nz! = -.1
  sx = 160 + nx! / nz! * 256
  sy = 150 + ny! / nz! * 256

END SUB

SUB mIdentity

  for i = 0 to 8
    matrix(i) = 0
  NEXT

  matrix(0) = 1
  matrix(4) = 1
  matrix(8) = 1


END SUB

SUB mMult
 
  FOR i = 0 TO 8
    tempmatrix(0, i) = matrix(i)
  NEXT

  FOR y = 0 TO 2
    FOR x = 0 TO 2
      v! = 0
      FOR a = 0 TO 2
        i1 = y * 3 + a
        i2 = a * 3 + x
        v! = v! + tempmatrix(0, i1) * tempmatrix(1, i2)
      NEXT
      i3 = y * 3 + x
      matrix(i3) = v!
    NEXT
  NEXT

END SUB

SUB mRotate (an!, plane)

  SELECT CASE plane
    CASE xplane
      tempmatrix(1, 0) = 1
      tempmatrix(1, 1) = 0
      tempmatrix(1, 2) = 0
      tempmatrix(1, 3) = 0
      tempmatrix(1, 4) = COS(rad * an!)
      tempmatrix(1, 5) = -SIN(rad * an!)
      tempmatrix(1, 6) = 0
      tempmatrix(1, 7) = SIN(rad * an!)
      tempmatrix(1, 8) = COS(rad * an!)
      mMult
    CASE yplane
      tempmatrix(1, 0) = COS(rad * an!)
      tempmatrix(1, 1) = 0
      tempmatrix(1, 2) = SIN(rad * an!)
      tempmatrix(1, 3) = 0
      tempmatrix(1, 4) = 1
      tempmatrix(1, 5) = 0
      tempmatrix(1, 6) = -SIN(rad * an!)
      tempmatrix(1, 7) = 0
      tempmatrix(1, 8) = COS(rad * an!)
      mMult
    CASE zplane
      tempmatrix(1, 0) = COS(rad * an!)
      tempmatrix(1, 1) = SIN(rad * an!)
      tempmatrix(1, 2) = 0
      tempmatrix(1, 3) = -SIN(rad * an!)
      tempmatrix(1, 4) = COS(rad * an!)
      tempmatrix(1, 5) = 0
      tempmatrix(1, 6) = 0
      tempmatrix(1, 7) = 0
      tempmatrix(1, 8) = 1
      mMult
  END SELECT


END SUB


sub mScale (factor!)
    for i = 0 to 8
        matrix(i) = matrix(i) * factor!
    next
end sub

SUB rotpnt (x0!, y0!, z0!, xr!, yr!, zr!, nx!, ny!, nz!)


  'xplane
  x1! = x0!
  y1! = z0! * SIN(rad * xr!) + y0! * COS(rad * xr!)
  z1! = z0! * COS(rad * xr!) + y0! * -SIN(rad * xr!)

  'yplane
  x2! = x1! * COS(rad * yr!) + z1! * -SIN(rad * yr!)
  y2! = y1!
  z2! = x1! * SIN(rad * yr!) + z1! * COS(rad * yr!)

  'zplane
  x3! = x2! * COS(rad * zr!) + y2! * -SIN(rad * zr!)
  y3! = x2! * SIN(rad * zr!) + y2! * COS(rad * zr!)
  z3! = z2!

  nx! = x3!
  ny! = y3!
  nz! = z3!

END SUB