#algo24: function to do the eigenfunction extension algorithm at a specific triangle.
#	T is a Sierpinski Gasket Table , L is the triangle the extension algorithm is to be performed at (i.e. [1,0,0,2])
#	fnum is the function number we are working on, and lambda is the discrete eigenvalue of the extended function.
#	note: lambda should never be 2 or 5

algo24:=proc(T,L,fnum,lambda)		 
	local c1,c2,Lb,L0,L1,L2,lvl,vm;
	
	c1:=(4-lambda)/((2-lambda)*(5-lambda));
	c2:=2/((2-lambda)*(5-lambda));

	#if (lambda=2 or lambda=5) then print(lambda); fi;

	Lb:=add_one(add_one(L,0),1);
	L0:=T[Lb,4];
	L1:=T[Lb,1];
	L2:=T[T[Lb,2],1];
	#print(Lb);
	T[Lb,fnum]:=c1*(T[L0,fnum]+T[L1,fnum])+c2*(T[L2,fnum]);

	Lb:=add_one(add_one(L,0),2);
	L0:=T[Lb,4];
	L1:=T[Lb,1];
	L2:=T[T[Lb,2],1];
	T[Lb,fnum]:=c1*(T[L0,fnum]+T[L1,fnum])+c2*(T[L2,fnum]);

	Lb:=add_one(add_one(L,1),2);
	L0:=T[Lb,4];
	L1:=T[Lb,1];
	L2:=T[T[Lb,2],1];
	T[Lb,fnum]:=c1*(T[L0,fnum]+T[L1,fnum])+c2*(T[L2,fnum]);

	#T[[0,2],fnum]:=c1*(T[[0],fnum]+T[[2],fnum])+c2*(T[[1],fnum]);
	#T[[1,2],fnum]:=c1*(T[[1],fnum]+T[[2],fnum])+c2*(T[[0],fnum]);
end;

#doAlgo24: executes algo24 at all triangles at a specific level. This is the level "level" refers to. lambda, fnum, and T
#	are the same as in algo24

doAlgo24:=proc(T,fnum,level ,lambda)
	local tris,t,vm;
	tris:=p([0,1,2],level);
	vm:=vertices(level+1);
	make_neighbors(T,level+1,vm);
	for t in tris do
		#print(t);
		algo24(T,t,fnum,lambda);
	end do;
end;

#algo24Wrapper: feed it startlevel, startlambda, and seqEp (of the form [1,-1,1,1]) and it runs doAlgo24 to a   
#	certain depth (determined by seqEp). The lambdas are generated from startlambda, and seqEp.

algo24Wrapper:=proc(T,fnum,startlevel,startlambda,seqEp)
	local newlambda,level,epsilon;
	doAlgo24(T,fnum,startlevel,startlambda);
	newlambda:=startlambda;
	level:=startlevel;
	for epsilon in seqEp do
	#	print(level);
		level:=level+1;
		newlambda:=(5+epsilon*sqrt(25-4*newlambda))/2;
		doAlgo24(T,fnum,level,newlambda);
	end do;
end;

#ourPlot: feed this a function (which must be defined to level "level" and in T at fnum) and a string plotname
#	and you get a nice, blue, plot of a function on the SG

ourPlot:=proc(T,level,fnum,plotname)
	local final_points;
	make_points(T,level);
	final_points:=SG_plot(T,level,fnum);
	#print(final_points);
	pointplot3d(final_points,color=blue,symbol=POINT,title=plotname);
end;

#ourPlotFancy: basically the same as ourplot, with color not specified and font changed.

ourPlotFancy:=proc(T,level,fnum,plotname)
	local final_points;
	make_points(T,level);
	final_points:=SG_plot(T,level,fnum);
	#print(final_points);
	pointplot3d(final_points,symbol=POINT,title=plotname,titlefont=[SYMBOL,20]);
	#pointplot3d(final_points,symbol=POINT,title='blah',titlefont=[SYMBOL,20]);
end;

#ourPlot2: same as ourPlot, but only plots on the two cells of level 1 at the bottom of the SG 

ourPlot2:=proc(T,level,fnum,plotname)
        local final_points;
        make_points(T,level);
        final_points:=OurNewSG_plot(T,level,fnum);
        #print(final_points);
        pointplot3d(final_points,color=blue,symbol=POINT,title=plotname);
end;

#LocPlot: same as ourPlot, but additionally you feed it a list of vertices to plot at.

LocPlot:=proc(T,level,fnum,vertlist,plotname)
        local final_points;
        make_points(T,level);
        final_points:=LocSG_plot(T,fnum,vertlist);
        #print(final_points);
        pointplot3d(final_points,color=blue,symbol=POINT,title=plotname);
end;

#OurNewSG_plot: function needed to make ourPlot2

OurNewSG_plot:=proc(T,m,funcno)
   local i,plot_pts,func, newfunc;
   	
   #func:=select(the_points,[indices(T)],5,0,m+1);
   func:=vertices(m);
   newfunc:=[];
   for i from 1 to nops(func) do
	if (sameCell(2,[1,2],func[i])) then newfunc:=[op(newfunc),func[i]]; end if;
   end do;		
   
   func:=newfunc;
   plot_pts:=[seq([T[func[i],5],T[func[i],6],T[func[i],funcno]],i=1..nops(func))];
 
   plot_pts;
 end; 

#LocSG_plot: function needed to make LocPlot

LocSG_plot:=proc(T,funcno,vertlist)
   local i,plot_pts,func, newfunc;
        
   func:=vertlist;
   plot_pts:=[seq([T[func[i],5],T[func[i],6],T[func[i],funcno]],i=1..nops(func))];
   plot_pts;
 end;

#ourAnimate: makes animation of a list (L) of function numbers.

ourAnimate:=proc(T,lvl,L,plotname)
	local final_points,pl,fnum;
	make_points(T,lvl);
	pl:=[];
	for fnum in L do
		final_points:=SG_plot(T,lvl,fnum);
	       #pl:=[op(pl),PLOT3D(POINTS(final_points),color=blue,symbol=POINT,title=plotname)];
		pl:=[op(pl),pointplot3d(final_points,color=blue,symbol=POINT,title=plotname)];
	       #pl:=[op(pl),PLOT3D(POINTS(final_points))];
	end do;
	display3d(pl,insequence=true);

	#display3d(pl,insequence=false);
end;

#displayToGif: feed this a filename (string) and run it before you plot. Instead of plotting to screen it plots to a file. 
#	note: different fonts don't work.

displayToGif:=proc(fname)
	plotsetup(gif,plotoutput=fname);
end;

#displayToPs same as displayToGif except Postscript file instead of GIF (ugly sometimes) and fonts do work.

displayToPs:=proc(fname)
	plotsetup(ps,plotoutput=fname,plotoptions=`noborder`);
end;

#defaultDisplay: run this to make it start plotting to the screen again.

defaultDisplay:=proc()
	plotsetup(default);
end;

#newAlgoWrapper: similar to algoWrapper, but the levels you feed it are more obvious to choose.

newAlgoWrapper:=proc(T,startlevel,fnum,startlambda,seqEp)
	local newlambda,level,epsilon;
        newlambda:=startlambda;
        level:=startlevel;
        for epsilon in seqEp do
                newlambda:=evalf((5+epsilon*sqrt(25-4*newlambda))/2,20);
                doAlgo24(T,fnum,level,newlambda);
        	level:=(level+1);
	end do;
end;

#makeLambda: feed this an initial lambda, and a list of plus and minus ones, and it generates a new lambda 
#	and returns it.

makeLambda:=proc(lambda, epslist)
	local newlambda, e;
	newlambda:=lambda;
	for e in epslist do
		newlambda:=(5+e*sqrt(25-4*newlambda))/2;
	end do;
	return(newlambda);
end;

#doLocAlgo24: same as LocAlgo24, but only runs it a list of the triangles that you feed it.
#	note: you need to run T:=verGasket(MAXLEVELYOUAREWORKINGAT) before this function because it accesses
#	a list of vertices stored at T[level,-1]. 

doLocAlgo24:=proc(T,fnum,level,lambda,trilist)
        local t;
        make_neighbors(T,level+1,T[(level+1),-1]);
        for t in trilist[level] do   
                algo24(T,t,fnum,lambda);
        end do;
end;

#LocAlgoWrapper: same as newAlgoWrapper, but only runs at triangles in the trilist you feed it. see note on doLocAlgo24
#	(must use verGasket first)

LocAlgoWrapper:=proc(T,startlevel,fnum,startlambda,seqEp,trilist)
        local newlambda,level,epsilon;
        newlambda:=startlambda;
        level:=startlevel;
        for epsilon in seqEp do
                newlambda:=evalf((5+epsilon*sqrt(25-4*newlambda))/2,20);
                doLocAlgo24(T,fnum,level,newlambda,trilist);
                level:=(level+1);
        end do; 
end;  

#doOneTriList: this is used in makeTriList

doOneTriList:=proc(tri,level)
	local t,tris,l;
        tris:=p([0,1,2],level-nops(tri));
	l:=[];
	for t in tris do
		l:=[op(l),[op(tri),op(t)]];
	end do;
	return(l);
end;


#makeTriList: feed this function a list of triangles and levels and it gives you (returns) all the triangles 
#	down to that level within those triangles. format of trilist is as follows.
#trilist=[[triangle,level],[triangle,level],...]

makeTriList:=proc(trilist)
	local tri,l,lst,i,maxlvl;
	lst:=[];
	maxlvl:=0;
	for tri in trilist do
		if(tri[2]>maxlvl) then
			maxlvl:=tri[2];
		end if;
	end do;

	for i from 1 to maxlvl do
		lst:=[op(lst),[]];
	end do;

	for tri in trilist do
		l:=tri[1];
		while nops(l)>0 do
			lst[nops(l)]:=[op(lst[nops(l)]),l];
			if nops(l)>1 then
				l:=[seq(l[i],i=1..nops(l)-1)];
			else
				l:=[];
			end if;
		end do;
		for i from nops(tri[1])+1 to tri[2] do
			lst[i]:=[op(lst[i]),op(doOneTriList(tri[1],i))];
		end do;
	end do;
	for i from 1 to nops(lst) do
		lst[i]:=removeDuplicates(lst[i]);
	end do;
	return (lst);
end;

#doOneVertList: used in makeVertList.

doOneVertList:=proc(T, tri,level)
	local v, list;
	list:=[];
	for v in T[level,-1] do
		if (sameCell(nops(tri),[op(tri),0,1],v)) then
			list:=[op(list),v];
		end if;
	end do;
	return(list);
end;

#makeVertList: feed this thing the same trilist you feed makeTriList, and it will give you (return) the vertices 
#	contained within those triangles. T is needed for the vertex list (must be initialized with verGasket() 

makeVertList:=proc(T,trilist)
	local t, list;
	list:=[];
	for t in trilist do
		list:=[op(list),op(doOneVertList(T,t[1],t[2]+1))];
	end do;
	list:=removeDuplicates(list);
	return(list);
end;

#lambdaLimit: this function is useful for finding the eigenvalues of the actual (non-discrete, that is) Laplacian on 
#	the SG. It's a short function, you can figure out what it does.

lambdaLimit:=proc(lambda,tol)
	local tmp, ev, lm, count;
	
	count:=1;
	lm:=makeLambda(lambda,[-1]);
	tmp:=lambda;
	ev:=5^(count)*lm;

	while(evalf(abs(tmp-ev))>evalf(tol)) do
		count:=count+1;
		tmp:=ev;
		lm:=makeLambda(lm,[-1]);
		ev:=5^(count)*lm;
	end do;

	return([ev,count]);
end;

#lambdaLimit2: this does the same thing as lambdaLimit, but slower. We used it to check lambdaLimit (laziness)

lambdaLimit2:=proc(lambda,tol)
	local count;
	count:=1;
	while evalf(abs( (5^count)*makeLambda(lambda,[seq(-1,i=1..count)]) - (5^(count-1))*makeLambda(lambda,[seq(-1,i=1..(count-1))]) ))>evalf(tol) do
		count:=count+1;
	end do;
	return [(5^count)*makeLambda(lambda,[seq(-1,i=1..count)]),count];
end;