clear all
close all
clc
name2='profilotdsimplesso';
name2pdf='profilotdsimplesso.pdf';
name2eps='profilotdsimplesso.eps';


AccMeth='column';
col=3;
polygon=[-1,-1; 1,-1; 1,1 ;-1,1];
rect=[-1,1,-1,1];
degmin=1;
degmax=13;
degstep=1;
degvalues=degmin:degstep:degmax;

for deg=degvalues;
    degmesh=floor(deg);
    x=cos(pi.*(1:-1/(2*deg):0));
    [X,Y]=meshgrid(x,x);
    wam_pts=[X(:),Y(:)];
%[fek,R1,R2,leb] = approxfek(deg,wam_pts,0,wam_pts,rect);
%[fek,nweights] = wamcomprex2(deg,wam_pts,0);
    fek=wam_pts;
    cardmesh=size(fek,1);
    s(deg)=(deg+1)*(deg+2)/2;
    V = chebvand(deg,fek,rect)./sqrt(cardmesh);
    [U,S]=svd(V);
    %diametro transfinito approx del quadrato in base di Cheb
    d(deg)=prod((diag(S)).^(3/(4*deg*s(deg))));
end
figure(1),plot(2.*d,'.-'),hold on;
title('Ratio transfinite diameter of the square/ same computation in the Cheb basis')

h=0.6; l=sqrt(2)/2;
if h+l>2
    error('h must be smaller')
end
for deg=degvalues;
    wam=polygon_wam([[-l;l;-0],[1-h-l;1-h-l;1-h]],deg);
    cardwam=size(wam,1);
    s(deg)=(deg+1)*(deg+2)/2;
    V = chebvand(deg,wam,rect)./sqrt(cardwam);
    [U,S]=svd(V);
    
    dd(deg)=prod((diag(S).^2).^(3/(4*deg*s(deg))));
end
tdS_approx(degvalues)=dd(degvalues)./(2.*d(degvalues));
tdS=1./(2*exp(1));

%%
[Table_tdS_approx_acc,TI]=...
    scaledvectorrho(tdS_approx(degvalues),degvalues,'originalsequence',10^-16,[]);
switch AccMeth
    case 'diagonal'
degvaluesAcc=diag(TI);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    tdS_acc(j)=Table_tdS_approx_acc(j,:,j);
% compute errors 
    ErrAccS(j)=abs(tdS_acc(end)-tdS);
    ErrAccRelS(j)=ErrAccS(j)/tdS;
end
    case 'column'
degvaluesAcc=TI(:,col);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    tdS_acc(j)=Table_tdS_approx_acc(j,:,col);
% compute errors 
    ErrAccS(j)=abs(tdS_acc(end)-tdS);
    ErrAccRelS(j)=ErrAccS(j)/tdS;
end
    case 'row'
degvaluesAcc=TI(col,:);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    tdS_acc(j)=Table_tdS_approx_acc(col,:,j);
% compute errors 
    ErrAccS(j)=abs(tdS_acc(end)-tdS);
    ErrAccRelS(j)=ErrAccS(j)/tdS;
    
end
end
%%
f3=figure(3);
semilogy(degvalues,abs(tdS-tdS_approx(degvalues)),'-')
hold on
semilogy(degvalues,abs(tdS-tdS_approx(degvalues))/tdS,'d-')
semilogy(degvaluesAcc,ErrAccS,'.-.r');
semilogy(degvaluesAcc,ErrAccRelS,'d-.r');
hold off
legend('abs err','rel err','abs err acc','rel err acc','Location','SouthWest')
%saveas(gcf,name2)
print(f3,name2pdf)
print(f3,name2eps)

%%
figure(4);
plot(degvalues,tdS_approx(degvalues));
hold on
plot(degvaluesAcc,tdS_acc,'r')
plot([degmin,degmax],[tdS,tdS])