close all
clear all
name='profilotdcerchio';
namepdf='profilotdcerchio.pdf';
nameeps='profilotdcerchio.eps';
name2='profilotdsimplesso';
name2pdf='profilotdsimplesso.pdf';
name2eps='profilotdsimplesso.eps';

clc
AccMeth='diagonal';
col=2;
polygon=[-1,-1; 1,-1; 1,1 ;-1,1];
rect=[-1,1,-1,1];
degmin=5;
degmax=25;
degstep=3;
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).^2).^(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')
% save('d','d');
% save('degmin','degmin');
% save('degmax','degmax');
% 
% load('d')

% calcolo td_approx del cerchio
for deg=degvalues;
    t=pi.*[0:1:2*deg]./(2*deg);
    r=cos(t);
    [T,R]=meshgrid(t,r');
    x=R(:).*cos(T(:));
    y=R(:).*sin(T(:));
    wam=[x,y];
    cardwam=length(x);
    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


td_approx(degvalues)=dd(degvalues)./(2.*d(degvalues));
td=sqrt(1./(2*exp(1)));
[Table_td_approx_acc,TI]=...
    scaledvectorrho(td_approx(degvalues),degvalues,'originalsequence',10^-16,[]);
%%
switch AccMeth
    case 'diagonal'
degvaluesAcc=diag(TI);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    td_acc(j)=Table_td_approx_acc(j,:,j);
% compute errors 
    ErrAcc(j)=abs(td_acc(end)-td);
    ErrAccRel(j)=ErrAcc(j)/td;
end
    case 'column'
degvaluesAcc=TI(:,col);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    td_acc(j)=Table_td_approx_acc(j,:,col);
% compute errors 
    ErrAcc(j)=abs(td_acc(end)-td);
    ErrAccRel(j)=ErrAcc(j)/td;
end
    case 'row'
degvaluesAcc=TI(col,:);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    td_acc(j)=Table_td_approx_acc(col,:,j);
% compute errors 
    ErrAcc(j)=abs(td_acc(end)-td);
    ErrAccRel(j)=ErrAcc(j)/td;
    
end
end
%%
f2=figure(2);
semilogy(degvalues,abs(td-td_approx(degvalues)),'-')
hold on
semilogy(degvalues,abs(td-td_approx(degvalues))/td,'d-')
semilogy(degvaluesAcc,ErrAcc,'.-r');
semilogy(degvaluesAcc,ErrAccRel,'d-.r');
hold off
legend('abs err','rel err','abs err acc','rel err acc','Location','SouthWest')
saveas(gcf,name)
print(f2,namepdf)
print(f2,nameeps)
%%
% calcolo td_approx del simplesso
for deg=degvalues;
    wam=polygon_wam([[-1;1;-1],[-1;-1;1]],deg);
    cardwam=length(x);
    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=2./(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);
semilogy(degvalues,tdS_approx(degvalues));
hold on
semilogy(degvaluesAcc,tdS_acc,'r')
semilogy([degmin,degmax],[tdS,tdS])