# 
#           Making a Peano Space-Filling Curve                   
# 
# 
# 
# First, we write a "helper routine" to allow us to  add a vector to a point without explictly writing both terms.
# 
> addvec := (p,q) -> [ p[1]+q[1], p[2]+q[2]]:

                  addvec := (p,q) -> [p[1] + q[1], p[2] + q[2]]
# Now, we write the main guts, which turns a segment into 9 segments, arranged appropriately.
# Notice that 4 of the segments are perpendicular to the original, so we need to compute "perp", which is the vector perpendicular to the given 
# one, and 1/3 as long.  the vector pq3 is the vector parallel to the input segment, but 1/3 the length.
# 
# Note that decoration only returns the 8 points making up the boxlike structure in the middle,  NOT the endpoints p & q.
> decoration:=proc(p,q) \
  local a,b,c,d,e,f, pq3, perp, perpneg;\
  pq3 := [ (q[1]-p[1])/3, (q[2]-p[2])/3 ];\
  perp:= [ pq3[2], -pq3[1] ];\
  perpneg:= [ -pq3[2], pq3[1] ];\
 \
  a := addvec(p, pq3);\
  b := addvec(a, perp);\
  c := addvec(b, pq3); \
  d := addvec(a, pq3);\
  e := addvec(a, perpneg);\
  f  := addvec(e, pq3);\
  RETURN(  a, b, c, d, a, e, f, d );\
end:
# 
# decorate takes a curve (actually a list of points), and calls decorate on each segment,  inserting the result between each pair of points.
#  
> decorate := proc ( incurve )\
local i, outcurve;\
\
outcurve := [incurve[1]];\
for i from 2 to nops(incurve) do\
     outcurve := [ op(outcurve),  decoration(incurve[i-1], incurve[i]), incurve[i] ];\
od;\
RETURN(outcurve);\
end:
--------------------------------------------------------------------------------
# 
# Now we are ready to make the curve.  first, we write the initiator, the unit segment.
> l0 := [[0,0], [1,0]]:
# 
# l1 is the first stage. 
> l1 := decorate(l0):\
plot(l1, axes=none, scaling=constrained);

   ** Maple V Graphics **

> l2:=decorate(l1):\
plot(l1, axes=none, scaling=constrained);

   ** Maple V Graphics **

# after two steps, we already see it beginning to fill out the square, although lets do a few more, just to see how it goes.
> l3:=decorate(l2):\
l4:= decorate(l3):\
plot(l4,axes=none);

   ** Maple V Graphics **

> 
--------------------------------------------------------------------------------
# Notice that for the way we wrote this, the orientation of the initiator didnt really matter; we can fill the unit square (instead of the diamond)
# by building on the diagonal instead.
# 
> plot(decorate(decorate(decorate([[0,0],[1,1]]))));

   ** Maple V Graphics **

# 
--------------------------------------------------------------------------------
# its a little difficult to SEE the way the curve winds its way through the square, so lets change it a little bit, so that the decorations dont quite meet 
# up.  Then we should be able to follow its path through.
> decoration:=proc(p,q) \
  local a,b,c,d,e,f, pq3, perp, perpneg, pq30, a2, d2;\
  pq3 := [ (q[1]-p[1])/3, (q[2]-p[2])/3 ];\
  pq30 := [ (q[1]-p[1])/20, (q[2]-p[2])/20 ];\
  perp:= [  .9*pq3[2], -.9*pq3[1] ];\
  perpneg:= [ - .9*pq3[2],  .9*pq3[1] ];\
 \
  a := addvec(p, pq3);\
  b := addvec(a, perp);\
  c := addvec(b, pq3); \
  d := addvec(a, pq3);\
  a2 := addvec(a,pq30);\
  e := addvec(a2, perpneg);\
  f  := addvec(e, pq3);\
  d2:= addvec(f,perp);\
  RETURN(  a, b, c, d, a2, e, f, d2 );\
end:
--------------------------------------------------------------------------------
# Now we can follow the way the curve winds through the square.   Of course, it no longer will fill out the whole square, but you can see the trend 
# of the original.
> plot(decorate(decorate(decorate([[0,0],[1,0]]))));

   ** Maple V Graphics **

# This square is skewed because I was not careful enough about the way I avoided self-intersections (some of the segments are longer than 
# others, etc.)  Oh well.  Fix it if you feel like it....
# 
--------------------------------------------------------------------------------
>