#########################################################################
#
#                   Peano's Space Filling Curve
#
#                    (Miguel A. Lerma, 7/8/95)
#                    (Modified for Maple V R.5 - 5/4/00)
#
# A sequence of curves is defined by a recursive procedure, starting 
# with a circle and using a set or rules to transform arcs of each 
# curve into list of arcs of the next curve. The arcs are of two types: 
# left (L) and right (R). The rules are the following:
#
# Start -> L L L L   (a circle)
# L -> R L L L R
# R -> R L R
#
#########################################################################


# Use just 3 digts in floating point.

Digits := 3:


# 'rule' gives how each arc is transformed from one iteration 
# to the next one. 

rule := x -> if x='L' then 'R','L','L','L','R' 
             elif x='R' then 'R','L','R'
             else ERROR 
             fi:


# 'list_of arcs' computes the list of arcs of the curve by iterating 'rule'.

list_of_arcs := proc(n)
                   if type(n,integer) then 
                      if n<1 then ERROR(`Argument shoud be positive`) fi:
                      if n=1 then ['L','L','L','L'] 
                      else map(rule,list_of_arcs(n-1))
                      fi
                   else ERROR(`Argument should be an integer`)
                   fi:
                end:


# 'left_arc' computes m points of a left arc of center c, radius r 
# and phase a.

left_arc := proc (c,r,a,m)
            local k,re,t:
               re:=[]: 
               for k from 1 to m do
                  t:=(k/m+a)*Pi/2:
               re:=[op(re),[evalf(r*cos(t)+c[1]),evalf(r*sin(t)+c[2])]]
               od:
               re:
            end:


# 'right_arc' computes m points of a right arc of center c, radius r 
# and phase a.

right_arc := proc (c,r,a,m)
             local k,re,t:
                re:=[]: 
                for k from 1 to m do
                   t:=(2-k/m+a)*Pi/2:
                   re:=[op(re),[evalf(r*cos(t)+c[1]),evalf(r*sin(t)+c[2])]]
                od:
                re:
             end:


# 'move_center' computes the center of the next arc

move_center := 
           proc(c,r,a,s):
              [c[1]+s*r*cos(a*Pi/2),c[2]+s*r*sin(a*Pi/2)]:
           end:  


# 'lopcurve' computes de list of points of the curve.
           
lopcurve := proc(n) 
            local r,c,a,t,k,N_P,loarcs,len,list_of_points:
  
              if n<1 or n>5 then 
                 ERROR(`Argument should be between 1 and 5.`)
              fi:
  
              loarcs:=list_of_arcs(n):
              len:=nops(loarcs):
  
                                   # number of points to graph of each arc
                                   # (it is greater for the first curves)
              N_P:=floor(10/n): 
  
              r:=1/2^n:            # radius
              if n<=1 then c:=[0,0]  
              else c:=[2*r,0]      # starting center
              fi:
              a:=0:                # starting phase
  
                              # Initialization of the list of points
  
              if loarcs[1]='L' then
                 list_of_points:=[[r,0],op(left_arc(c,r,a,N_P))]:
                 a:=a+1:
              elif loarcs[1]='R' then 
                 list_of_points:=[[r,0],op(right_arc(c,r,a,N_P))]:
                 a:=a-1:
              else ERROR
              fi:
  
              for k from 2 to len do
  
                              # If there is a change in the type of arc, 
                              # move the center to its new position.
  
                 if loarcs[k-1]='L' and loarcs[k]='R' then 
                   c:=move_center(c,r,a,2.) 
                 fi:
                 if loarcs[k-1]='R' and loarcs[k]='L' then 
                    c:=move_center(c,r,a,-2.) 
                 fi:
  
                              # Compute the next arc.            
  
                 if loarcs[k]='L' then 
                    list_of_points:=
                              [op(list_of_points),op(left_arc(c,r,a,N_P))]:
                    a:=a+1:
                 fi:
                 if loarcs[k]='R' then 
                    list_of_points:=
                              [op(list_of_points),op(right_arc(c,r,a,N_P))]:
                    a:=a-1:
                 fi:
  
              od:
              list_of_points:   # Return the list of points.
           end:



# 'plot_curve' plots the curve

plot_curve := n -> plot(lopcurve(n),-1..1,-1..1,style=LINE,axes=BOXED):