Making a Peano 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]]:

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:
The procedure decoration 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]]:

We use l1 for the first stage of the Peano curve. 

l1 := decorate(l0):
plot(l1, axes=none, scaling=constrained);

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

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);
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]]))));

It is a little difficult to SEE the way the curve winds its way through the square, so lets change decoration a little bit, so that the decorations dont quite meet up. Then we should be able to follow its path through the square. 

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]]))));

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....