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 pushMatrix ()
declare sub popMatrix ()
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), pnt(4,7) as single

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
  connected as integer
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

type leafType
    x   as single
    y   as single
    z   as single
    size    as single
end type

DIM branches(MAX_BRANCHES) AS branchType, leaves(MAX_LEAVES) AS leafType


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
leafscale! = 3.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
        branches(branch).nodes(node).connected = branches(branch).nodes(node).connected + 1
        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    


for leaf = 0 to MAX_LEAVES
    try = 0
    do
        try = try + 1
        if try = 100 then exit do
        branch = int(rnd * nBranches)
        node = int(rnd * branches(branch).nNodes)
        globalNode = branches(branch).nodes(node).globalNode
        pass = 0
        if globalNode > MAX_NODES / 2 and globalNode < MAX_NODES - 1 then pass = 1
        for i = 0 to leaf - 1
            ff = 0
            if leaves(i).x = branches(branch).nodes(node).x then ff = ff + 1
            if leaves(i).y = branches(branch).nodes(node).y then ff = ff + 1
            if leaves(i).z = branches(branch).nodes(node).z then ff = ff + 1
            if ff = 3 then pass = 0
        next
        if pass = 1 then exit do
        
    loop
    
    leaves(leaf).x = branches(branch).nodes(node).x
    leaves(leaf).y = branches(branch).nodes(node).y
    leaves(leaf).z = branches(branch).nodes(node).z
    leaves(leaf).size = branches(branch).nodes(node).connected
    if leaves(leaf).size > 2 then leaves(leaf).size = 2
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
    
    pnt(0,0) = -.5
    pnt(0,1) = -.5
    pnt(1,0) = -.5
    pnt(1,1) = .5
    pnt(2,0) = .5
    pnt(2,1) = .5
    pnt(3,0) = .5
    pnt(3,1) = -.5
    
    for leaf = 0 to MAX_LEAVES
        pushMatrix
        mIdentity
        mRotate -xan!, xplane
        mRotate -yan!, yplane
        
        for p = 0 to 3
            x0! = pnt(p,0) * leaves(leaf).size * leafscale!
            y0! = pnt(p,1) * leaves(leaf).size * leafscale!
            z0! = 0
            pnt(p,2) = x0! * matrix(0) + y0! * matrix(3) + z0! * matrix(6) + leaves(leaf).x
            pnt(p,3) = x0! * matrix(1) + y0! * matrix(4) + z0! * matrix(7) + leaves(leaf).y
            pnt(p,4) = x0! * matrix(2) + y0! * matrix(5) + z0! * matrix(8) + leaves(leaf).z
        next
        popMatrix
        
        for p = 0 to 3
            getXY pnt(p,2), pnt(p,3), pnt(p,4), sx, sy
            pnt(p,5) = sx
            pnt(p,6) = sy
        next
        
        
        line (pnt(0,5), pnt(0,6)) - (pnt(1,5), pnt(1,6)), 2
        line (pnt(0,5), pnt(0,6)) - (pnt(3,5), pnt(3,6)), 2
        line (pnt(2,5), pnt(2,6)) - (pnt(1,5), pnt(1,6)), 2
        line (pnt(2,5), pnt(2,6)) - (pnt(3,5), pnt(3,6)), 2
    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 pushMatrix
    for i = 0 to 8
        tempmatrix(2,i) = matrix(i)
    next
end sub

sub popMatrix
    for i = 0 to 8
        matrix(i) = tempmatrix(2,i)
    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