 |
Well, here is the Maple code for the basics on the gasket for your viewing pleasure:
>
#this procedure generates all words of length m in L letters...
p:=proc(L,m)
option remember;
local L1;
if m = 1 then L
else
L1 := p(L,m-1);
[op(map(add_one,L1,0)),op(map(add_one,L1,1)),op(map(add_one,L1,2))]
fi;
end;
>
add_one:=(u,d)-> [op(u),d];
>
#This function returns all vertices in Vm in a big list
vertices:= proc(m)
local i,outlist;
option remember;
if m = 1 then RETURN([[0],[1],[2],[0,1],[0,2],[1,2]]) fi;
if m=2 then
[[0],[1],[2],[0,1],[0,2],[1,2],[0,0,1],[0,0,2],[0,1,2],[1,0,1],[1,0,2],[1,1,2],[2,0,1],[2,0,2],[2,1,2]]; else
outlist:=[[0] , [1] , [2] , [0,1], [0,2], [1,2]];
for i from 2 to m do
outlist:=[op(outlist),op(map(add_pair, p([[0],[1],[2]], (i-1)), [0,1])), op(map(add_pair, p([[0],[1],[2]], (i-1)), [0,2])), op(map(add_pair, p([[0],[1],[2]], (i-1)), [1,2]))];
od;
outlist;
#[[0] , [1] , [2] , [0,1], [0,2], [1,2], op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,1])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,2])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [1,2]))];
fi;
end;
#vertices_m creates just the vertices in Vm, excluding Vk, k<m...
vertices_m:=proc(m)
option remember;
if m = 0 then RETURN([[0],[1],[2]]);fi;
if m = 1 then RETURN([[0,1],[0,2],[1,2]]);fi;
if m > 1 then [op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,1])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,2])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [1,2]))]; fi;
end;
>
add_pair := (u,v) -> [op(u), op(v)];
>
#This function returns a table of the gasket out to level m...
gasket:=proc(m)
local i, vm, T;
vm:= vertices(m);
T:=table();
for i from 1 to nops(vm) do
T[vm[i], 0] := vm[i] od;
T;
end;
>
#This function makes stores the neighbors of each vertex out to level m in
#the gasket table T. The neighbors are ordered clockwise; i.e., the 1st neighbor
#of 0,1 in V3 is 1,0,0,1, etc...
make_neighbors:=proc(T,m,vm)
local i,j,depth,nb1,nb2,nb3,nb4,iverts,last_verts;
#first, we have to deal with the special case of V1...
if m=1 then
print(`Yo, V1 is STOOOOOOPID!`);
T[[0],1]:=[0,1];
T[[0],2]:=[0,2];
T[[1],1]:=[1,2];
T[[1],2]:=[0,1];
T[[2],1]:=[0,2];
T[[2],2]:=[1,2];
T[[0,1],3]:=[0,2];
T[[0,1],4]:=[0];
T[[0,1],1]:=[1];
T[[0,1],2]:=[1,2];
T[[1,2],3]:=[0,1];
T[[1,2],4]:=[1];
T[[1,2],2]:=[0,2];
T[[1,2],1]:=[2];
T[[0,2],1]:=[0];
T[[0,2],2]:=[0,1];
T[[0,2],3]:=[1,2];
T[[0,2],4]:=[2];
#hard coding reverse neighbors...
T[T[[0],1],4] := [0];
T[T[[0],2],1] := [0];
T[T[[1],1],4] := [1];
T[T[[1],2],1] := [1];
T[T[[2],1],4] := [2];
T[T[[2],2],1] := [2];
#done with hard coding (yay!)
fi;
if m <>1 then #but if we are not in V1...
#hard coding the neighbors for V0,V1
T[[0],1]:=[seq(0,i=1..m),1];
T[[0],2]:=[seq(0,i=1..m),2];
T[[1],1]:=[seq(1,i=1..m),2];
T[[1],2]:=[seq(1,i=1..m-1),0,1];
T[[2],1]:=[seq(2,i=1..m-1),0,2];
T[[2],2]:=[seq(2,i=1..m-1),1,2];
T[[0,1],3]:=[0,seq(1,i=1..m-2),1,2];
T[[0,1],4]:=[0,seq(1,i=1..m-2),0,1];
T[[0,1],1]:=[1,seq(0,i=1..m-2),0,1];
T[[0,1],2]:=[1,seq(0,i=1..m-2),0,2];
T[[1,2],3]:=[1,seq(2,i=1..m-2),0,2];
T[[1,2],4]:=[1,seq(2,i=1..m-2),1,2];
T[[1,2],2]:=[2,seq(1,i=1..m-2),0,1];
T[[1,2],1]:=[2,seq(1,i=1..m-2),1,2];
T[[0,2],1]:=[0,seq(2,i=1..m-2),0,2];
T[[0,2],2]:=[0,seq(2,i=1..m-2),1,2];
T[[0,2],3]:=[2,seq(0,i=1..m-2),0,1];
T[[0,2],4]:=[2,seq(0,i=1..m-2),0,2];
#hard coding reverse neighbors...
T[T[[0],1],4] := [0];
T[T[[0],2],1] := [0];
T[T[[1],1],4] := [1];
T[T[[1],2],1] := [1];
T[T[[2],1],4] := [2];
T[T[[2],2],1] := [2];
T[T[[0,1],1],4] := [0,1];
T[T[[0,1],2],1] := [0,1];
T[T[[0,1],3],4] := [0,1];
T[T[[0,1],4],1] := [0,1];
T[T[[1,2],1],4] := [1,2];
T[T[[1,2],2],1] := [1,2];
T[T[[1,2],3],4] := [1,2];
T[T[[1,2],4],1] := [1,2];
T[T[[0,2],1],4] := [0,2];
T[T[[0,2],2],1] := [0,2];
T[T[[0,2],3],4] := [0,2];
T[T[[0,2],4],1] := [0,2];
#done with hard coding (yay!)
# First, we compute the neighbors of all pts in Vk, 1<k<m. Notice that all of the neighbors
# will be in Vm...
for i from 2 to m-1 do
#doing the 01 endings first...
iverts:=pick_vertices01(vm,i);
for j from 1 to nops(iverts) do
nb1:=[op(iverts[j][1..i-1]),op(T[[0,1],1][1..m-i]),0,1];
nb2:=[op(iverts[j][1..i-1]),op(T[[0,1],2][1..m-i]),0,2];
nb3:=[op(iverts[j][1..i-1]),op(T[[0,1],3][1..m-i]),1,2];
nb4:=[op(iverts[j][1..i-1]),op(T[[0,1],4][1..m-i]),0,1];
T[iverts[j],1]:=nb1: T[iverts[j],2]:=nb2: T[iverts[j],3]:=nb3: T[iverts[j],4]:=nb4;
T[nb1,4]:=iverts[j]: T[nb2,1]:=iverts[j]: T[nb3,4]:=iverts[j]: T[nb4,1]:=iverts[j];
od;
#doing the 12 endings second...
iverts:=pick_vertices12(vm,i);
for j from 1 to nops(iverts) do
nb1:=[op(iverts[j][1..i-1]),op(T[[1,2],1][1..m-i]),1,2];
nb2:=[op(iverts[j][1..i-1]),op(T[[1,2],2][1..m-i]),0,1];
nb3:=[op(iverts[j][1..i-1]),op(T[[1,2],3][1..m-i]),0,2];
nb4:=[op(iverts[j][1..i-1]),op(T[[1,2],4][1..m-i]),1,2];
T[iverts[j],1]:=nb1: T[iverts[j],2]:=nb2: T[iverts[j],3]:=nb3: T[iverts[j],4]:=nb4;
T[nb1,4]:=iverts[j]: T[nb2,1]:=iverts[j]: T[nb3,4]:=iverts[j]: T[nb4,1]:=iverts[j];
od;
#doing the 02 endings third...
iverts:=pick_vertices02(vm,i);
for j from 1 to nops(iverts) do
nb1:=[op(iverts[j][1..i-1]),op(T[[0,2],1][1..m-i]),0,2];
nb2:=[op(iverts[j][1..i-1]),op(T[[0,2],2][1..m-i]),1,2];
nb3:=[op(iverts[j][1..i-1]),op(T[[0,2],3][1..m-i]),0,1];
nb4:=[op(iverts[j][1..i-1]),op(T[[0,2],4][1..m-i]),0,2];
T[iverts[j],1]:=nb1: T[iverts[j],2]:=nb2: T[iverts[j],3]:=nb3: T[iverts[j],4]:=nb4;
T[nb1,4]:=iverts[j]: T[nb2,1]:=iverts[j]: T[nb3,4]:=iverts[j]: T[nb4,1]:=iverts[j];
od;
od;
#Here we do the neighbors of the vertices in Vm...
#Notice that we only have to assign 2 neighbors because each pt in
#Vm has 2 neighbors in Vm and 2 neighbors in Vk, k<m.
last_verts:=pick_vertices01(vm,m);
for j from 1 to nops(last_verts) do
nb2:=[op(last_verts[j][1..m-1]),1,2];
nb3:=[op(last_verts[j][1..m-1]),0,2];
T[last_verts[j],2]:=nb2: T[last_verts[j],3]:=nb3;
od;
last_verts:=pick_vertices12(vm,m);
for j from 1 to nops(last_verts) do
nb2:=[op(last_verts[j][1..m-1]),0,2];
nb3:=[op(last_verts[j][1..m-1]),0,1];
T[last_verts[j],2]:=nb2: T[last_verts[j],3]:=nb3;
od;
last_verts:=pick_vertices02(vm,m);
for j from 1 to nops(last_verts) do
nb2:=[op(last_verts[j][1..m-1]),0,1];
nb3:=[op(last_verts[j][1..m-1]),1,2];
T[last_verts[j],2]:=nb2: T[last_verts[j],3]:=nb3;
od;
fi;
T;
end;
>
#picks vertices out of vm of length i with ending [0,1]
pick_vertices01:=proc(vm,i)
select(len_list_eq, vm,i,[0,1]);
end;
>
len_list_eq:=proc(l,i,ending)
evalb(nops(l) = i+1 and l[nops(l)-1..nops(l)] = ending);
end;
>
#like pick_vertices01...
pick_vertices12:=proc(vm,i)
select(len_list_eq, vm,i,[1,2]);
end;
>
#like pick_vertices01...
pick_vertices02:=proc(vm,i)
select(len_list_eq, vm,i,[0,2]);
end;
>
#This function stores the x and y coordinates of each vertex
#in the 5 and 6 slots of T, respectively
make_points:=proc(T,m)
local i;
#func:=select(the_points,[indices(T)],0,3);
# first, we have to manually set the values of x and y out to V1...
T[[0],5]:=0;
T[[1],5]:=.5;
T[[2],5]:=-.5;
T[[0],6]:=sqrt(3)/2;
T[[1],6]:=0;
T[[2],6]:=0;
T[[0,1],5]:=.25;
T[[0,1],6]:=sqrt(3)/4;
T[[1,2],5]:=0;
T[[1,2],6]:=0;
T[[0,2],5]:=-.25;
T[[0,2],6]:=sqrt(3)/4;
#gen_points recursively creates the rest of the coordinates out to Vm...
gen_points(T,m-1,[0],[0,1],[0,2]);
gen_points(T,m-1,[1],[1,2],[0,1]);
gen_points(T,m-1,[2],[0,2],[1,2]);
end;
>
the_points:=proc(L,funcno,minlen,maxlen)
evalb(L[2] = funcno and nops(L[1]) >= minlen and nops(L[1]) <= maxlen);
end;
>
#This function is called by make_points. It calls itself recursively 'n stuff.
gen_points:=proc(T,m,add0,add1,add2)
if m <> 0 then
#So the general strategy is as follows: We can't assign values to the next level of
#vertices without knowing the orientation of the triangle (because we need to know whether we need to assign
#to a 0,1 point, 1,2 point, etc.). So what we do is call the function with 3 addresses
#and make it so that the 0th address (add0) is one of the vertices of the outside triangle
#(i.e., going from the V0 triangle, the triangles would be T1=([0],[0,1],[0,2]),T2=([1],[1,2],[0,1])
#and so forth). Thus, by looking at add1 of this triangle, we can determine the orientation of
#the triangle by checking whether add1 ends in 0,1, 0,2 or 1,2. Then we assign the values
#accordingly, after which we call gen_points recursively, making sure to keep the orientation
#of the 3 addresses in the recursive calls. We use this same method to generate harmonic functions
#and n-harmonic functions, as well as several other things. Notice that this recursive process
#has it's base case as V1, therefore requiring 3 seperate calls to gen_points from make_points.
#This is due to the idiosyncratic addressing of the vertices of V0. Oh well....
if add1[nops(add1)-1..nops(add1)] = [0,1] then
T[[op(add1[1..nops(add1)-2]),0,0,1],5] := .5 * T[add0,5] + .5 * T[add1,5];
T[[op(add1[1..nops(add1)-2]),0,1,2],5] := .5 * T[add1,5] + .5 * T[add2,5];
T[[op(add1[1..nops(add1)-2]),0,0,2],5] := .5 * T[add0,5] + .5 * T[add2,5];
T[[op(add1[1..nops(add1)-2]),0,0,1],6] := .5 * T[add0,6] + .5 * T[add1,6];
T[[op(add1[1..nops(add1)-2]),0,1,2],6] := .5 * T[add1,6] + .5 * T[add2,6];
T[[op(add1[1..nops(add1)-2]),0,0,2],6] := .5 * T[add0,6] + .5 * T[add2,6];
gen_points(T,m-1,add0,[op(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2]);
gen_points(T,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2]),0,0,1]);
gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2]),0,0,2],[op(add1[1..nops(add1)-2]),0,1,2]);
fi;
if add1[nops(add1)-1..nops(add1)] = [1,2] then
T[[op(add1[1..nops(add1)-2]),1,1,2],5] := .5 * T[add0,5] + .5 * T[add1,5];
T[[op(add1[1..nops(add1)-2]),1,0,2],5] := .5 * T[add1,5] + .5 * T[add2,5];
T[[op(add1[1..nops(add1)-2]),1,0,1],5] := .5 * T[add0,5] + .5 * T[add2,5];
T[[op(add1[1..nops(add1)-2]),1,1,2],6] := .5 * T[add0,6] + .5 * T[add1,6];
T[[op(add1[1..nops(add1)-2]),1,0,2],6] := .5 * T[add1,6] + .5 * T[add2,6];
T[[op(add1[1..nops(add1)-2]),1,0,1],6] := .5 * T[add0,6] + .5 * T[add2,6];
gen_points(T,m-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1]);
gen_points(T,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2]);
gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2]),1,0,1],[op(add1[1..nops(add1)-2]),1,0,2]);
fi;
if add1[nops(add1)-1..nops(add1)] = [0,2] then
T[[op(add1[1..nops(add1)-2]),2,0,2],5] := .5 * T[add0,5] + .5 * T[add1,5];
T[[op(add1[1..nops(add1)-2]),2,0,1],5] := .5 * T[add1,5] + .5 * T[add2,5];
T[[op(add1[1..nops(add1)-2]),2,1,2],5] := .5 * T[add0,5] + .5 * T[add2,5];
T[[op(add1[1..nops(add1)-2]),2,0,2],6] := .5 * T[add0,6] + .5 * T[add1,6];
T[[op(add1[1..nops(add1)-2]),2,0,1],6] := .5 * T[add1,6] + .5 * T[add2,6];
T[[op(add1[1..nops(add1)-2]),2,1,2],6] := .5 * T[add0,6] + .5 * T[add2,6];
gen_points(T,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..nops(add1)-2]),2,1,2]);
gen_points(T,m-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2]);
gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add1[1..nops(add1)-2]),2,0,1]);
fi;
fi;
end;
>
#This function returns a list of x,y,z coordinates from a function
#stored in T. This list can be sent directly into pointplot3d for plotting
SG_plot:=proc(T,m,funcno)
local i,plot_pts,func;
#func:=select(the_points,[indices(T)],5,0,m+1);
func:=vertices(m);
plot_pts:=[seq([T[func[i],5],T[func[i],6],T[func[i],funcno]],i=1..nops(func))];
plot_pts;
end;
>
#This is used to animate a time series of basis values.
#It takes the parameters in vals, which is a list of lists,
#with each sublist corresponding to a particular spline function
#representation. Funcno is just some place in which to store the
#actual values of the function. The output is a list of plot
#structures, which can then be sent directly to display, and if you
#use insequence = true, it will animate for you...
make_animation:=proc(T,m_basis,m_func,vals,funcno,harm_funcno)
local S,vm,SG_plots,final_points,inds,i;
SG_plots:=[];
S:=gasket(m_basis);
vm:=vertices(m_basis);
make_neighbors(S,m_basis,vm);
make_points(T,m_func);
for i from 1 to nops(vals) do
print(`computing plot number `,i);
make_pretty_picture(T,S,m_basis,m_func,vm,funcno,harm_funcno,vals[i]);
final_points:=SG_plot(T,m_func,funcno);
SG_plots:=[op(SG_plots),pointplot3d(final_points,color=BLUE,symbol=POINT)];
od;
SG_plots;
end;
>
>
>
#This function will take a line of the sierpinski gasket, and
#returns, in order, a list of the vertices along that line. It is
#used in conjunction with SG_plot_line in order to create line
#plots of a function along a line of the SG. The line is defined
#by starting at startpt. From there, we go to first_neighbor. After
#that, we continue following rest_neighbor, until the ending (ie, last
#2 digits) change, thereby indicating we have gotten to the end of the line.
#for example, the line from [0,2] to [1,2] would be gotten by
#line_vertices(4,[0,2],3,1), since we would start by taking the 3rd neighbor,
#and then continue along the 1st neighbor (the endings would all be 0,1)
#until we got to 1,2, at which point, the ending of rest_neighbor would be
#1,2, indicating that we have gotten to the end of the line.
line_vertices:=proc(m_func,startpt,first_neighbor,rest_neighbor)
local lineT,vm,outlist,ending,last_pt;
lineT:=gasket(m_func);
vm:=vertices(m_func);
make_neighbors(lineT,m_func,vm);
outlist:=[startpt,lineT[startpt,first_neighbor]];
ending:=last_digits(lineT[startpt,first_neighbor]);
last_pt:=lineT[startpt,first_neighbor];
while last_digits(lineT[last_pt,rest_neighbor])=ending do
outlist:=[op(outlist),lineT[last_pt,rest_neighbor]];
last_pt:=lineT[last_pt,rest_neighbor];
od;
outlist:=[op(outlist),lineT[last_pt,rest_neighbor]];
outlist;
end;
>
>
>
#This function is just like make_animation, except that it makes 2d line plots
#of the function along a line defined by startpt, first_neighbor, rest_neighbor.
#You also need to input linelength to tell the program how long the line should
#be. For instance, the line from [0,2] to [1,2] would require linelength to be
#1/2. For more information about how startpt, first_neighbor, and rest_neighbor
#define a line, look at line_vertices above.
make_animation_line:=proc(T,m_basis,m_func,vals,funcno,harm_funcno,startpt,first_neighbor,rest_neighbor,linelength)
local linepts,S,vm,SG_plots,final_points,inds,i;
SG_plots:=[];
S:=gasket(m_basis);
vm:=vertices(m_basis);
make_neighbors(S,m_basis,vm);
make_points(T,m_func);
linepts:=line_vertices(m_func,startpt,first_neighbor,rest_neighbor);
for i from 1 to nops(vals) do
print(`computing plot number `,i);
make_pretty_picture(T,S,m_basis,m_func,vm,funcno,harm_funcno,vals[i]);
final_points:=SG_plot_line(T,m_func,funcno,linepts,linelength);
SG_plots:=[op(SG_plots),plot(final_points,color=BLUE,symbol=POINT)];
od;
SG_plots;
end;
>
>
#This function is just like SG_plot, except that it creates a 2d
#line plot along a line which is defined by linepts, and is of length
#linelength.
SG_plot_line:=proc(T,m,funcno,linepts,linelength)
local i,plot_pts,numpts;
numpts:=nops(linepts);
plot_pts:=[seq([(i-1)*linelength/(numpts-1),T[linepts[i],funcno]],i=1..numpts)];
plot_pts;
end;
|
