read "itall.txt":
with(stats):
with(describe):


lvl:=5:
phatlvl:=8:
numpsis:=(3^(lvl+1)-3)/2:
path:=cat("../data2/","level",lvl,"/"):
T:=gasket(phatlvl):

fd:=fopen(cat(path,phatlvl,"phat",lvl,"dirintegrals.txt"),WRITE);
fprintf(fd,"l1, l2, and integral data for the %d phat %d dirichelt\n",phatlvl,lvl);
fclose(fd);

l1list:=[];
l2list:=[];
intlist:=[];

for i from 1 to numpsis do
	readFunc(T,12,cat(path,phatlvl,"phat",lvl,"dirichlet.",i));
	funcProduct(T,phatlvl,12,12,13);
	l2:=simp_rule(T,phatlvl,13);
	funcAbs(T,phatlvl,12,13);
	l1:=simp_rule(T,phatlvl,13);
	integ:=simp_rule(T,phatlvl,12);
	fd:=fopen(cat(path,phatlvl,"phat",lvl,"dirintegrals.txt"),APPEND);
	fprintf(fd,"%d:L1=%20.20g\tL2=%20.20g\tintegral=%20.20g\n",i,l1,l2,integ);
	fclose(fd);
	intlist:=[op(intlist),integ];
	l1list:=[op(l1list),l1];
	l2list:=[op(l2list),l2];
end do;

fd:=fopen(cat(path,phatlvl,"phat",lvl,"dirintegrals.txt"),APPEND);
fprintf(fd,"min:L1=%20.20g\tL2=%20.20g\tintegral=%20.20g\n",min(op(l1list)),min(op(l2list)),min(op(intlist)));
fprintf(fd,"max:L1=%20.20g\tL2=%20.20g\tintegral=%20.20g\n",max(op(l1list)),max(op(l2list)),max(op(intlist)));
fprintf(fd,"average:L1=%20.20g\tL2=%20.20g\tintegral=%20.20g\n",mean(l1list),mean(l2list),mean(intlist));
fclose(fd);