clear all;
tic
N=51;  F=N*(N+1)/2; %U=4.5;  
LL3(1:F)=0; LL2(1:F)=0; 
PF(1:F)=0;  c1=zeros(N); SSS(1:F,1:5)=0;

for ii=1:1, %be(ii)=2*pi*(ii-1)/100; 
        U=0.5+(ii-1)/10; UU(ii)=U;
R1=zeros(F); P(1:F)=0;
n=(N+1)/2;  m=n-1; f=m+n*(n-1)/2;  P(f)=1;
%n=(N-1)/2;  m=n-1; f=m+n*(n-1)/2;  P(f)=1;
for n=1:N, ee(n)=2.05*sin(1+2*pi*(n-(N+1)/2)*(sqrt(5)-1)/2); end


for n=1:N, 
    for m=1:n, f=m+n*(n-1)/2; 
R1(f,f)=ee(n)+ee(m);
    end
   m=n; f=m+n*(n-1)/2; R1(f,f)=ee(n)+ee(n)+U;
end



for n=4:N-1,
    for m=2:n-2, 
f=m+n*(n-1)/2; m1=m+1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m-1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n+1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n-1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
    end
end

for n=3:N-1, m=n-1; 
f=m+n*(n-1)/2; m1=m+1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);
f=m+n*(n-1)/2; m1=m-1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n+1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n-1; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);
end

for n=2:N-1, m=n; 
f=m+n*(n-1)/2; m1=m-1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);
f=m+n*(n-1)/2; m1=m; n1=n+1; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);
end

for n=3:N-1, m=1; 
f=m+n*(n-1)/2; m1=m+1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n+1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n-1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
end

n=2; m=1;
f=m+n*(n-1)/2; m1=m+1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);
f=m+n*(n-1)/2; m1=m; n1=n+1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n-1; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);

n=1; m=1;
f=m+n*(n-1)/2; m1=m; n1=n+1; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);

for m=2:N-2, n=N; 
f=m+n*(n-1)/2; m1=m+1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m-1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n-1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
end

n=N; m=N-1;
f=m+n*(n-1)/2; m1=m+1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);
f=m+n*(n-1)/2; m1=m-1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n-1; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);

n=N; m=N;
f=m+n*(n-1)/2; m1=m-1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=sqrt(2);

n=N; m=1;
f=m+n*(n-1)/2; m1=m+1; n1=n; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;
f=m+n*(n-1)/2; m1=m; n1=n-1; f1=m1+n1*(n1-1)/2; R1(f,f1)=1;


[V,s]=eig(R1);
E1=real(diag(s));
[E10,d]=sort(E1);
for p1=1:F, VV(:,p1)=V(:,d(p1)); end


for p1=1:F, LLm(p1)=0; LLn(p1)=0;
for n=1:N, 
for m=1:n
        f=m+n*(n-1)/2;  c1(n,m)=VV(f,p1); cc1(n,m,p1)=VV(f,p1); LLm(p1)=LLm(p1)+m*(c1(n,m))^2; LLn(p1)=LLn(p1)+n*(c1(n,m))^2;
end
end
%  LLn=0;
% for n=1:N, 
% for m=1:n
%         f=m+n*(n-1)/2;  LLn=LLn+n*(c1(n,m))^2;
% end
% end
%  LL2(p1)=0; LL3(p1)=0;
% for n=1:N, 
% for m=1:n
%         f=m+n*(n-1)/2;  LL2(p1)=LL2(p1)+(m+n-LLm-LLn)^2*(c1(n,m))^2;
%         LL3(p1)=LL3(p1)+((m-LLm)^2+(n-LLn)^2)*(c1(n,m))^2;
% end
% end
% 
% 
% 



for n=1:N, c11(n)=0; 
for m=1:n 
        c11(n)=c11(n)+(c1(n,m))^2; 
end
end

for m=1:N, c22(m)=0; 
for n=m:N 
        c22(m)=c22(m)+(c1(n,m))^2; 
end
end

c33=(c11+c22)/2;
 PPR=0;
for m=1:N, 
        PPR=PPR+c33(m)^2;
end


%c2=(c1+c1')/2;
%PPR=sum(sum(c2.^4))^0.5;

PF(p1)=1/PPR; SSS(p1,ii)=1/PPR;
end
ii
RR(ii,1)=U;
[RR(ii,2) max_index]=max(PF);
RR(ii,3)=LLn(max_index);
RR(ii,4)=LLm(max_index);
RR(ii,5)=E10(max_index);
% FF(ii)=max(LL2); 
% CC(ii)=max(LL3);
end
PFout(:,1)=PF';
PFout(:,2)=E10;
fid=fopen('www.dat','w');
fprintf(fid,'%8.4f\n',PFout); 
fclose(fid);

toc

% PFout(:,1)=PF';
% PFout(:,2)=E10;
% save RR301_2_05_0.txt RR -ascii
% save RF301_2_05_0.txt PFout -ascii