! 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