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
        
    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
            x0! = branches(branch).nodes(node).x
            y0! = branches(branch).nodes(node).y
            z0! = branches(branch).nodes(node).z
            
            if node > 0 then
                x1! = branches(branch).nodes(node - 1).x
                y1! = branches(branch).nodes(node - 1).y
                z1! = branches(branch).nodes(node - 1).z
            else
                x1! = branches(branch).root.x
                y1! = branches(branch).root.y
                z1! = branches(branch).root.z
            end if
            
            getXY x0!, y0!, z0!, screenX0, screenY0
            getXY x1!, y1!, z1!, screenX1, screenY1
            
            line (screenX0, screenY0) - (screenX1, screenY1), 6
            
        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