! SUNFLOWER
! size -- tells the size of the circle of influence of a
! leaf at the specified height
def size(t) = .3*exp(-t/10)
! flags
let plotting = 1
! arrays to hold information about the leaves
dim leaf_height(500),leaf_turn(500),leaf_size(500)
! 'tolerance' tells how fussy to be about leaf location
let tolerance = 1e-3
! 'pt_siz' tells the plotting sub 'point_me' how big pts are
let pt_siz = .1
! remember to declare functions early and ofter
declare function find_turn,find_x,bump
if plotting=1 then call plot_init
! start off with one lone leaf at the bottom of the stalk
call zap_leaves
call new_leaf(0,.25,size(0))
! make that two
!call new_leaf(0,.75,size(0))
! now move up the stalk, adding new leaves whenever we
! find room.
let height = 0
do
call next_place(height,next_height,next_turn)
call new_leaf(next_height,next_turn,size(next_height))
let height = next_height
loop
! next_place -- tell the next place at or above the
! specified height that a leaf can go.
sub next_place(h,h1,t1)
if find_turn(h) <= 1 then
let h1 = h
let t1 = find_turn(h)
exit sub
end if
let h0 = h
let h1 = h0 + 3*size(h0)
let t1 = find_turn(h1)
if t1 > 1 then
print "oops"
stop
end if
do until h1-h0 < tolerance
let h2 = (h0+h1)/2
let t2 = find_turn(h2)
if t2 > 1 then
let h0 = h2
let t0 = t2
else
let h1 = h2
let t1 = t2
end if
loop
end sub
! find_turn -- tries to find the first available place for
! a leaf at the specified height. Returns a value >1 if
! no leaf will fit at this height.
def find_turn(h)
let find_turn=find_x(h,size(h),leaves,leaf_turn,leaf_height,leaf_size)
end def
! find_x -- find the first x in [0,1] so that a circle of
! radius r about (x,y) on the cylinder [0,1]*R
! doesn't intersect any of the cnum
! circles whose centers and radii are specified in the
! arrays cx(), cy(), cr(). Return a value >1 if there
! is no such x.
! WARNING: This routine assumes that cy(i)+cr(i) increases
! with i.
def find_x(y,r,cnum,cx(),cy(),cr())
! start x off at 0 and keep trying to bump it to the right
let x_to_bump = 0
do
! compute an amount to bump by
for cindex = cnum to 1 step -1 ! (call it a hunch)
if y - r > cy(cindex)+cr(cindex) then exit do ! are we above it?
for cdisp = -1 to 1 ! mod 1, each circle is really 3 circles
let x_bumped = x_to_bump + bump(cx(cindex)+cdisp-x_to_bump,cy(cindex)-y,cr(cindex)+r)
if x_bumped > x_to_bump then let bumped = 1 else let bumped = 0
let x_to_bump = x_bumped
if bumped = 1 then exit for
next cdisp
if bumped = 1 then exit for
next cindex
loop until bumped = 0 or x_to_bump > 1
let find_x = x_to_bump
end def
! bump -- return the first x>=0 for which (x,0)
! is not inside the circle of radius r about (x0,y0).
def bump(x0,y0,r)
if x0^2 + y0^2 >= r^2 then
let bump = 0
else
let bump = x0 + sqr(r^2-y0^2)
end if
end def
! zap_leaves -- start with a clean stalk
sub zap_leaves
let leaves = 0
end sub
! new_leaf -- sprout a new leaf that the specified point,
! and having the specified effective size
sub new_leaf(new_height, new_turn, new_size)
let leaves = leaves+1
let leaf_height(leaves) = new_height
let leaf_turn(leaves) = new_turn
let leaf_size(leaves) = new_size
if plotting=1 then call plot_leaf(new_height,new_turn,new_size)
end sub
! print_leaves -- prints out the list of leaves
sub print_leaves
print "number","height","turn","size"
for isub = 1 to leaves
print isub, leaf_height(isub), leaf_turn(isub), leaf_size(isub)
next isub
print
end sub
! plot_init -- set up for plotting
sub plot_init
open #1:screen .4,.6,.1,.9
box lines 0,1,0,1
set window 0,1,0,10
end sub
! plot_leaf -- show the location of a leaf and its sphere
! of influence
sub plot_leaf(y,x,s)
for displace = -1 to 1
call circle_me(x+displace,y,s)
call point_me(x,y)
next displace
end sub
! circle_me, point_me
sub circle_me(x,y,s)
box circle x-s,x+s,y-s,y+s
end sub
sub point_me(x,y)
plot x,y
box area x-pt_siz/8,x+pt_siz/8,y-pt_siz/4,y+pt_siz/4
end sub
end