#This script approximates test functions using the PSI functions.
#One of MAKEBIHARMONIC MAKESQUEAREBIHARMONIC MAKESQUAREEIGEN or MAKEEIGEN should be set to true. 
#Then in order to approximate that function, you set APPFUNCTION to true 

#The script plots the function, the aprroximate, and their difference (along woth the min/max of each


read "itall.txt":

#This is the function number to store the test function in
TESTFN:=12: 

MAKEBIHARMONIC:=false:
	#these are the values to set the normal derivatives to at the boundary vertices
	BLAP1:=0:	
	BLAP2:=50:
	BLAP3:=-100:
MAKESQUAREBIHARMONIC:=true: #makes the square of a biharmonic function
MAKESQUAREEIGEN:=false:
	#level of eigenfunction
	SELVL:=4: 
	#function number
	SENO:=95: 
MAKEEIGEN:=false:
	#level of eigenfunction
	TELVL:=5: 
	#function number
	TENO:=349: 
APPFUNCTION:=true:
	#function number to store the approximate in
	APPFN:=13: 
	#temporary storage space (function number)
	TMPFN:=14: 
#selects the level of the Psi functions to sample with
lvl:=4: 
phatlvl:=8:
path:="../data/":

T:=gasket(phatlvl);

verts:=subsop(1=NULL,2=NULL,3=NULL,ourVertices(lvl)):

if (MAKEBIHARMONIC) then
	T[[0],8]:=0:		
	T[[1],8]:=0:
	T[[2],8]:=0:
	T[[0],9]:=BLAP1;
	T[[1],9]:=BLAP2;
	T[[2],9]:=BLAP3;
	
	n_harm(T,R,2,phatlvl,[0],[1],[2],8):
	copyFunc(T,phatlvl,8,TESTFN):

	print("Done Making Bi-Harmonic Function");
end if;

if (MAKESQUAREBIHARMONIC) then
	T[[0],8]:=0:		
	T[[1],8]:=0:
	T[[2],8]:=0:
	T[[0],9]:=BLAP1;
	T[[1],9]:=BLAP2;
	T[[2],9]:=BLAP3;
	
	n_harm(T,R,2,phatlvl,[0],[1],[2],8):
	copyFunc(T,phatlvl,8,TESTFN):
	funcProduct(T,phatlvl,TESTFN,TESTFN,TESTFN):	

	print("Done Making Square Bi-Harmonic Function");
end if;

if (MAKESQUAREEIGEN) then
	readFunc(T,TESTFN,cat(path,"level",SELVL,"/",phatlvl,"phat",SELVL,"eigen.",SENO));
	funcProduct(T,phatlvl,TESTFN,TESTFN,TESTFN);
end if:

if (MAKEEIGEN) then
	readFunc(T,TESTFN,cat(path,"level",TELVL,"/",phatlvl,"phat",TELVL,"eigen.",TENO));
end if:

if (APPFUNCTION) then
	clearFunc(T,phatlvl,APPFN);
	for i from 1 to nops(verts) do
		co:=T[verts[i],TESTFN];
		readFunc(T,TMPFN,cat(path,"level",lvl,"/",phatlvl,"phat",lvl,"psi.",i));
		funcSumWithConstant(T,phatlvl,co,APPFN,TMPFN,APPFN);
	end do;
end if:

ourPlot(T,phatlvl,TESTFN,"Test Function");
print(cat("Max=",findFuncMax(T,phatlvl,TESTFN),", Min=",findFuncMin(T,phatlvl,TESTFN)));
ourPlot(T,phatlvl,APPFN,"Approximation of Test Function");
print(cat("Max=",findFuncMax(T,phatlvl,APPFN),", Min=",findFuncMin(T,phatlvl,APPFN)));
funcDifference(T,phatlvl,APPFN,TESTFN,TMPFN);
ourPlot(T,phatlvl,TMPFN,"Approximation - Test Function");
print(cat("Max=",findFuncMax(T,phatlvl,TMPFN),", Min=",findFuncMin(T,phatlvl,TMPFN)));

print(cat("level= ",lvl));