The following text is the program. The executable files are also available in a zipped file. Copy to your hard drive, unzip, and execute helix.exe.
' In order of appearance in program:
Dim trial, scaled, CH
Dim ELEV, BWANG, ANG, XMIN
Dim GRAD, NG, WHORLS, StartSeed
Dim Yr, YrD, Zr, ZrD ' For Sub Convert
Dim sigma, GROW, RAD, aXMIN
Dim seed, Bseed, Nit
Dim Pscale
Dim Hmin, Hmax, Vmin, Vmax
Dim pi, x, y, z
Dim j As Single
Dim jr ' To carry angle "j" to Sub Convert
Dim jj, k, kk
Dim tip(2000, 4), Ntip(2000, 4), Ttip(2000, 4)
' 4 tip subscripts = x, y, z, rotation angle "j"
Dim xp, yp, zp ' Plotting coordinates from Sub Convert
Dim g
Dim X1, Y1, z1
Dim xx, yy, zz
Dim A1, A2
Dim tip1x, tip1y, tip2x, tip2y
Private Sub Command1_Click()
trial = 0
Cls
scaled = 0
Do
WindowState = 2
Label1.Visible = False
If CH <> 1 Then ' CH = 1 indicates default values changed by user input
'----------Set Default Parameter and Variable Values----------
ELEV = 0.17 ' MODEL PARAMETERS of McKinney & Raup (1982)
BWANG = 60
ANG = 45
XMIN = 20
' OTHER VARIABLES:
GRAD = 1 ' Growth gradient at distal tip of helix
NG = 25 ' Number of growth increments
WHORLS = 2.8 ' Number of whorls in the helix
StartSeed = 54321 ' Initial random number seed
YrD = 0 ' Initial rotation around the y-axis
ZrD = 0 ' Initial rotation around the z-axis
ElseIf CH = 1 Then End If
'----------Constants (can be changed in the source code)----------
sigma = 0.05 ' Limit for random number generator
GROW = 0.5 ' Linear growth per increment
RAD = 2 ' Radius of central helix
aXMIN = XMIN / 2 ' Scaled to McKinney & Raup (1982) simulations
'----------Increment random number seed ("Try again" button)----------
Bseed = StartSeed + Nit
seed = Bseed
Nit = Nit + 1
Randomize (seed)
'----------Scaling Trial Run (without plotting; checks screen scale)----------
trial = trial + 1
If scaled = 0 Then
Pscale = 1
Hmin = 10000: Hmax = -10000: Vmin = 10000: Vmax = -10000
ElseIf scaled <> 0 Then End If
pi = 4 * Atn(1) ' Establish value of pi
Offset = pi / 4 ' Default orientation, McKinney & Raup (1982, fig. 2)
Command1.Caption = "Please wait . . . "
'----------Make Central Helix----------
x = RAD * Cos(Offset)
y = -RAD * Sin(Offset)
z = -ELEV * 360 / 5 * Offset / (2 * pi)
jr = Offset: Call convert
If trial > 1 Then Line (xp, zp) - (xp, zp)
For j = Offset To WHORLS * 2 * pi + 2 * pi / 36 + 0.01 Step 2 * pi / 36
' Line segments 10 degrees each
x = RAD * Cos(j)
y = -RAD * Sin(j)
z = -ELEV * 360 / 5 * j / (2 * pi)
jr = j: Call convert
If trial > 1 Then Line - (xp, zp)
jj = j
Next j
'----------Make Initial Growing Tips ANG degrees apart----------
'----------(jj is total rotation through all whorls in radians)----------
k = 1
For j = Offset to jj + 0.01 Step 2 * pi / (360 / ANG)
k = k + 1
x = (GROW + RAD) * Cos(j) ' New tip
y = -(GROW + RAD) * Sin(j)
jr = j: Call convert: If trial > 1 Then Line (xp, zp) - (xp, zp)
If trial <= 1 Then
If zp < Vmin Then Vmin = zp: If zp > Vmax Then Vmax = zp
If xp > Hmax Then Hmax = xp: If xp < Hmin Then Hmin = xp
ElseIf trial > 1 Then End If
tip(k, 1) = x ' Record growing tip
tip(k, 2) = y
tip(k, 4) = j ' Rotation angle for this tip
x = RAD * Cos(j) ' Point on helix
y = -RAD * Sin(j)
jr = j: Call convert: If trial > 1 Then Line -(xp, zp)
If trial <= 1 Then
If zp < Vmin Then Vmin = zp: If zp > Vmax Then Vmax = zp
If xp > Hmax Then Hmax = xp: If xp < Hmin Then Hmin = xp
ElseIf trial > 1 Then End If
kk = j
Next j
'----------(invisible boundaries to avoid edge effects at termini of
helix)----------
' start boundary
x = (GROW + RAD) * Cos(Offset - pi / (360 / ANG))
y = -(GROW + RAD) * Sin(Offset - pi / (360 / ANG))
' record starting boundary
tip(1, 1) = x: tip(1, 2) = y: tip(1, 4) = -pi / (360 / ANG) + Offset
' end boundary
x = (GROW + RAD) * Cos(kk + pi / (360 / ANG))
y = -(GROW + RAD) * Sin(kk + pi / (360 / ANG))
' record ending boundary
tip(k + 1, 1) = x: tip(k + 1, 2) = y: tip(k + 1, 4) = kk + pi / (360 / ANG)
Ntips = k + 1
'----------BEGIN MAIN PROGRAM (Grow All Tips in Colony)----------
For k = 2 To NG ' Calculate each growth increment "k"
For j = 1 To Ntips ' For each growing tip at "j"
If tip(j, 4) < jj - 2 * pi Or jj < 2 * pi Then
' Normal growth if WHORLS <= 1
Ntip(j, 1) = (k * GROW + RAD) * Cos(tip(j, 4))
Ntip(j, 2) = -(k * GROW + RAD) * Sin(tip(j, 4))
ElseIf tip(j, 4) >= jj - 2 * pi And GRAD > 0 Then
' Reduced growth for final whorl (using GRAD)
Ntip(j, 1) = (k * GROW * (1 - GRAD * (tip(j, 4) - _
(jj - 2 * pi)) / (2 * pi)) + RAD) * Cos(tip(j, 4))
Ntip(j, 2) = -(k * GROW * (1 - GRAD * (tip(j, 4) - _
(jj - 2 * pi)) / (2 * pi)) + RAD) * Sin(tip(j, 4))
End If
Ntip(j, 4) = tip(j, 4)
' Plot new tip, end of old tip
x = tip(j, 1): y = tip(j, 2): jr = tip(j, 4): Call convert
If j > 1 And j < Ntips And trial > 1 Then Line (xp, zp) - (xp, zp)
If trial <= 1 Then
If zp < Vmin Then Vmin = zp
If zp > Vmax Then Vmax = zp
If xp > Hmax Then Hmax = xp
If xp < Hmin Then Hmin = xp
ElseIf trial > 1 Then End If
' Start of new tip
x = Ntip(j, 1): y = Ntip(j, 2): jr = tip(j, 4): Call convert
If j > 1 And j < Ntips And trial > 1 Then Line -(xp, zp)
If trial <= 1 Then
If zp < Vmin Then Vmin = zp: If zp > Vmax Then Vmax = zp
If xp > Hmax Then Hmax = xp: If xp < Hmin Then Hmin = xp
ElseIf trial > 1 Then End If
Next j
' Rewrite tip(j, 1), tip(j, 2), & tip(j, 4)
For j = 1 To Ntips
tip(j, 1) = Ntip(j, 1)
tip(j, 2) = Ntip(j, 2)
tip(j, 4) = Ntip(j, 4)
Next j
'----------Loop for Bifurcations due to XMIN----------
LastSplit = 0
Do ' Census all growing tips for those that need to bifurcate
For j = LastSplit + 2 To Ntips - 1
x = tip(j - 1, 1): y = tip(j - 1, 2): jr = tip(j - 1, 4): Call convert
' x-y-z position of point at angle = j - 1
X1 = xp: Y1 = yp: z1 = zp
' x-y-z position of point at angle = j + 1
x = tip(j + 1, 1): y = tip(j + 1, 2): jr = tip(j + 1, 4): Call convert
xx = xp - X1: yy = yp - Y1: zz = zp - z1
' Distance between points in 3D
dist = Sqr(xx * xx + yy * yy + zz * zz)
' Get a random number from normal distribution
Call Rng
' Limit deviation from XMIN to 0.5
If Abs(g) > 0.5 Or g = 0 Then Call Rng
' Bifurcate if the following condition holds:
If dist / Pscale > (1 - g) * (aXMIN / 5) Then Exit For
Next j
If j >= Ntips - 1 Then Exit Do
' Calculate bifurcate at "j":
A1 = tip(j, 4) - (1 / 6) * (tip(j + 1, 4) - tip(j - 1, 4))
A2 = tip(j, 4) + (1 / 6) * (tip(j + 1, 4) - tip(j - 1, 4))
If tip(j, 4) < jj - 2 * pi Or jj < 2 * pi Then
' Normal growth if WHORLS <= 1
tip1x = (k * GROW + RAD) * Cos(A1)
tip1y = -(k * GROW + RAD) * Sin(A1)
tip2x = (k * GROW + RAD) * Cos(A2)
tip2y = -(k * GROW + RAD) * Sin(A2)
ElseIf tip(j, 4) >= jj - 2 * pi And GRAD > 0 Then
' Reduced growth for final whorl (using GRAD)
tip1x = (k * GROW * (1 - GRAD * (tip(j, 4) - _
(jj - 2 * pi)) / (2 * pi)) + RAD) * Cos(A1)
tip1y = -(k * GROW * (1 - GRAD * (tip(j, 4) - _
(jj - 2 * pi)) / (2 * pi)) + RAD) * Sin(A1)
tip2x = (k * GROW * (1 - GRAD * (tip(j, 4) - _
(jj - 2 * pi)) / (2 * pi)) + RAD) * Cos(A2)
tip2y = -(k * GROW * ( 1 - GRAD * (tip(j, 4) - _
(jj - 2 * pi)) / (2 * pi)) + RAD) * Sin(A2)
End If
' Draw crossbar
x = tip1x: y = tip1y: jr = A1: Call convert
If k < NG And trial > 1 Then Line (xp, zp) - (xp, zp)
x = tip2x: y = tip2y: jr = A2: Call convert
If k < NG And trial > 1 Then Line -(xp, zp)
' Save tips in Ttip( )
For m = 1 To Ntips
Ttip(m, 1) = tip(m, 1)
Ttip(m, 2) = tip(m, 2)
Ttip(m, 4) = tip(m, 4)
Next m
' Rewrite tip( ), inserting new tips
tip(j, 1) = tip1x: tip(j, 2) = tip1y: tip(j, 4) = A1
tip(j + 1, 1) = tip2x: tip(j + 1, 2) = tip2y: tip(j + 1, 4) = A2
For n = j + 1 To Ntips
tip(n + 1, 1) = Ttip(n, 1)
tip(n + 1, 2) = Ttip(n, 2)
tip(n + 1, 4) = Ttip(n, 4)
Next n
Ntips = Ntips + 1
LastSplit = j + 1
' Return to bifurcation census:
Loop Until j = Ntips - 1
' Return for next growth increment:
Next k
If trial = 1 Then scaled = 1
' Scaling trial complete:
Loop While trial = 1
'----------END MAIN PROGRAM----------
'----------Screen display of parameter values----------
CurrentX = 100: CurrentY = 8000
Print " ANG = "; ANG
Print " XMIN = "; XMIN
Print " BWANG = "; BWANG
Print " ELEV = "; ELEV
Print " "
Print " Whorls: "; WHORLS
Print " Number of growth increments: "; NG
Print " Growth gradient: "; GRAD
Print " Rotation: "; ZrD
Print " Tilt: "; YrD
Print " Starting Seed = "; StartSeed
CurrentY = 7500
Beep
Command1.Caption = "Run Again"
Command2.Visible = True
Command3.Visible = True
Command4.Visible = True
Command5.Visible = True
Command6.Visible = True
Command7.Visible = True
Command8.Visible = True
Command9.Visible = True
Command10.Visible = True
Command11.Visible = True
Command12.Visible = True
Command13.Visible = True
scaled = 0
End Sub
'----------COMMAND BUTTONS----------
Private Sub Command10_Click()
'----------Random Number Seed Button----------
Lt = InputBox ("Random number seed", "HIT ENTER for no change", StartSeed, 4000,
4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of seed
StartSeed = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command11_Click()
'----------Rotation Angle (degrees) Button----------
Lt = InputBox("Rotation in degrees about vertical axis (negative for
clockwise)", _
"HIT ENTER for no change", ZrD, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of rotation
ZrD = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command12_Click()
'----------Tilt Angle (degrees) Button----------
Lt = InputBox("Tilt in degrees (negative for clockwise)", _
"HIT ENTER for no change", YrD, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of tilt
YrD = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command13_Click()
'----------Growth Gradient Button----------
' GRAD = 0 for no growth gradient
' GRAD = 1 for full gradient (i.e., zero growth at distal tip)
Lt = InputBox("Growth gradient", "HIT ENTER for no change", GRAD, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of gradient
GRAD = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command2_Click()
End
End Sub
Private Sub Command3_Click()
Command1.Visible = False
Command2.Visible = False
Command3.Visible = False
Command4.Visible = False
Command5.Visible = False
Command6.Visible = False
Command7.Visible = False
Command8.Visible = False
Command9.Visible = False
Command10.Visible = False
Command11.Visible = False
Command12.Visible = False
Command13.Visible = False
Form1.PrintForm
Printer.EndDoc
Command1.Visible = True
Command2.Visible = True
Command3.Visible = True
Command4.Visible = True
Command5.Visible = True
Command6.Visible = True
Command7.Visible = True
Command8.Visible = True
Command9.Visible = True
Command10.Visible = True
Command11.Visible = True
Command12.Visible = True
Command13.Visible = True
End Sub
Private Sub Command4_Click()
'----------Model Parameter ANG Button----------
Lt = InputBox("ANG in degrees", "HIT ENTER for no change", ANG, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of ANG
ANG = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command5_Click()
'----------Model Parameter XMIN Button----------
Lt = InputBox("XMIN", "HIT ENTER for no change", XMIN, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of XMIN
XMIN = Lt
aXMIN = XMIN / 2
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command6_Click()
'----------Model Parameter BWANG Button----------
Lt = InputBox("BWANG in degrees", "HIT ENTER for no change", BWANG, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of BWANG
BWANG = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command7_Click()
'----------Model Parameter ELEV Button----------
Lt = InputBox("ELEV", "HIT ENTER for no change", ELEV, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of ELEV
ELEV = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command8_Click()
'----------Number of Whorls Button----------
Lt = InputBox("Number of Whorls", "HIT ENTER for no change", WHORLS, 4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of WHORLS
WHORLS = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub Command9_Click()
'----------Number of Growth Increments----------
Lt = InputBox("Number of growth increments", "HIT ENTER for no change", NG,
4000, 4000)
If Len(Lt) = 0 Then Exit Sub
CH = 1 ' bypass default value of increments
NG = Lt
Command1.SetFocus
Label1.Visible = True
End Sub
Private Sub convert()
'----------Convert Subroutine----------
' Given x, y, BWANG, and angle in radians (jr), returns plotting coordinates: xp
& zp
' plus yp for evaluating distance between branches (and for rotation and tilt)
xh = RAD * Cos(jr) ' Point on helix at this rotation angle (jr)
yh = -RAD * Sin(jr)
zh = -ELEV * 360 / 5 * jr / (2 * pi) ' Elevation of point on helix
d = Sqr((y - yh) ^ 2 + (x - xh) ^ 2) ' Distance from x, y to point on helix
dp = d * Cos(pi * (90 - BWANG) / 180) ' Distance from x, y to helix after
elevation
zp = zh - d * Sin(pi * (90 - BWANG) / 180) ' z-value for plotting
xp = xh + dp * Cos(jr) ' x-value for plotting
yp = yh - dp * Sin(jr) ' y-value for computing 'dist' and plotting (IF NEEDED)
' Rotation around y-axis (Yr) and/or z-axis (Zr):
Yr = YrD ' Negative for clockwise
Zr = ZrD ' Negative for clockwise
Yr = Yr - 90: Yr = pi * Yr / 180: Zr = Zr + 180: Zr = pi * Zr / 180
yyy = xp * Sin(Zr) + yp * Cos(Zr)
xxx = xp * Cos(Zr) * Sin(Yr) - yp * Sin(Zr) * Sin(Yr) + zp * Cos(Yr)
zzz = xp * Cos(Zr) * Cos(Yr) - yp * Sin(Zr) * Cos(Yr) - zp * Sin(Yr)
xp = xxx: yp = yyy: zp = zzz
If scaled = 0 Then Exit Sub
' Scale xp, zp for plotting:
If Abs(Hmax - Hmin) > Abs(Vmax - Vmin) Then _
Pscale = 7000 / (Hmax - Hmin) Else Pscale = 7000 / (Vmax - Vmin)
xp = 4000 + (xp - Hmin) * Pscale
zp = 8000 - (Vmax - zp) * Pscale
End Sub
Private Sub Form_Load()
WindowState = 2
Command2.Visible = False
Command3.Visible = False
Command4.Visible = False
Command5.Visible = False
Command6.Visible = False
Command7.Visible = False
Command8.Visible = False
Command9.Visible = False
Command10.Visible = False
Command11.Visible = False
Command12.Visible = False
Command13.Visible = False
Label1.Visible = False
End Sub
Private Sub Rng()
'----------Random Number Subroutine----------
' Random number generator (crude but works).
' Yields normally distributed values ('g') with mean of zero and
' standard deviation = 'sigma'.
r3 = 2
Do Until r3 < 1
r1 = 2 * Rnd - 1
r2 = 2 * Rnd - 1
r3 = r1 * r1 + r2 * r2
Loop
r3 = Sqr((-2 * Log(r3)) / r3)
g = sigma * (r1 * r3)
End Sub