algo24:=proc(T,L,fnum,lambda)
#option remember;
local c1,c2,Lb,L0,L1,L2,lvl,vm;

#lvl:=nops(L)+1;
#vm:=vertices(lvl);
#make_neighbors(T,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:=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:=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:=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:=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:=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;

ourPlotBlack:=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=black,symbol=POINT);
end;

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:=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:=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:=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:=proc(fname)
	plotsetup(gif,plotoutput=fname);
end;

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

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

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:=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:=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:=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:=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;

#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:=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:=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:=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:=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;