global U V k1 k2 k3 k4 k21 k22 k23 k24 a t0 m1 m2 m3 m4 mu d0 d2;
t0=1;U=1.0;V=0.0;n = 1;k1 = 0.1;k2 = 0.1;k3 = 0.1;k4 = 0.1;k21=0.0;k22=0.0;k23=0.0;k24=0.0;
a=0.001;m1=2;m2=2;m3=2;m4=2;mu=0.0;nd0=4;nd1=0.00;nd2=0.0011;nd3=0.0001;N0=2^12;d0=0;d2=0;
for l =1:1
% t0=0.01.*(l-1)
%mu=0.5*(l-1)
% U=0.01.*(l-1)
% a=0.01*l*0.01;
% d0=2.*pi.*0.01.*(l-1)
% d2=-2.*pi.*0.01*(l-1)
% d=0.01*0.5*l
for k = 1:n
r=0.0.*exp(i.*rand(1));
g=[-0.5,0.5,0.5,-0.5,-0.5,0.5,0.5,-0.5,0.5,0.5,0.5,-0.5,-0.5,0.5,0.5,-0.5,0.0+r,0.0+r,0.0+r,0.0+r,0.0+r,0.0+r,0.0+r,0.0+r];
%[t,u]=ode45('loreq12',[N0],[gm',0.01+r,0.01+r,0.01+r,0.01+r,0.01+r,0.01+r,0.01+r,0.01+r]');
[t,u]=ode45('loreq17',[N0],g');
%tau31/(real(u(1,1)-u(1,4)+u(1,5)-u(1,8)).*10);
%tau=3000;
xi = 0.0:1:N0;
y(:,1,k) = interp1(t,u(:,1),xi);
y(:,2,k) = interp1(t,u(:,2),xi);
y(:,3,k) = interp1(t,u(:,3),xi);
y(:,4,k) = interp1(t,u(:,4),xi);
y(:,5,k) = interp1(t,u(:,5),xi);
y(:,6,k) = interp1(t,u(:,6),xi);
y(:,7,k) = interp1(t,u(:,7),xi);
y(:,8,k) = interp1(t,u(:,8),xi);
y(:,9,k) = interp1(t,u(:,9),xi);
y(:,10,k) = interp1(t,u(:,10),xi);
y(:,11,k) = interp1(t,u(:,11),xi);
y(:,12,k) = interp1(t,u(:,12),xi);
y(:,13,k) = interp1(t,u(:,13),xi);
y(:,14,k) = interp1(t,u(:,14),xi);
y(:,15,k) = interp1(t,u(:,15),xi);
y(:,16,k) = interp1(t,u(:,16),xi);
y(:,17,k) = interp1(t,u(:,17),xi);
y(:,18,k) = interp1(t,u(:,18),xi);
y(:,19,k) = interp1(t,u(:,19),xi);
y(:,20,k) = interp1(t,u(:,20),xi);
%y(:,13,k) = 8000000000000./((tau.*tau.*(xi - 7.*10.^(-2)).*(xi - 7.*10.^(-2))+1).*tau);
%end
%for j=1:N0/0.001;
% h(1)=y(j,1);
% h(2)=y(j,2);
% h(3)=y(j,3);
% h(4)=y(j,4);
% h(5)=y(j,5);
% h(6)=y(j,6);
% h(7)=y(j,7);
% h(8)=y(j,8);
% h(9)=y(j,9);
% h(10)=y(j,10);
% h(11)=y(j,11);
% h(12)=y(j,12);
% h(13)=y(j,13);
% h(14)=y(j,14);
% h(15)=y(j,15);
% h(16)=y(j,16);
%X=[h(1) h(9); h(13) h(5)];
%b(j)=-trace(X.*logm(X));
% c(j)=y(j,1)-y(j,4)+y(j,5)-y(j,8);
%% Y=-log(sum(eigs(X)));
%% b(j)=Y;
%end
%for k =1:n
%% for j =0:400/0.0025
% h(:,1)=+y(:,1,k)/n;
% h(:,2)=+y(:,2,k)/n;
% h(:,3)=+y(:,3,k)/n;
% h(:,4)=+y(:,4,k)/n;
% h(:,5)=+y(:,5,k)/n;
% h(:,6)=+y(:,6,k)/n;
% h(:,7)=+y(:,7,k)/n;
%% h(:,8)=+y(:,8,k)/n;
% %end
end
w1=fft(y(:,1)+y(:,4) +y(:,5) +y(:,8) ).*conj(fft(y(:,1) +y(:,4) +y(:,5) +y(:,8)));
w2=fft(y(:,1)-y(:,4) +y(:,5) -y(:,8) ).*conj(fft(y(:,1) -y(:,4) +y(:,5) -y(:,8)));
w3=fft(y(:,1)+y(:,4) -y(:,5) -y(:,8) ).*conj(fft(y(:,1) +y(:,4) -y(:,5) -y(:,8)));
w4=fft(y(:,1)-y(:,4) -y(:,5) +y(:,8) ).*conj(fft(y(:,1) -y(:,4) -y(:,5) +y(:,8)));
w5=fft(y(:,9)+y(:,12)+y(:,13)+y(:,15)).*conj(fft(y(:,9)+y(:,12)+y(:,13)+y(:,15)));
w6=fft(y(:,9)-y(:,12)+y(:,13)-y(:,15)).*conj(fft(y(:,9)-y(:,12)+y(:,13)-y(:,15)));
w7=fft(y(:,9)+y(:,12)-y(:,13)-y(:,15)).*conj(fft(y(:,9)+y(:,12)-y(:,13)-y(:,15)));
w8=fft(y(:,9)-y(:,12)-y(:,13)+y(:,15)).*conj(fft(y(:,9)-y(:,12)-y(:,13)+y(:,15)));
w9=y(:,1)+y(:,4)+y(:,5)+y(:,8);
w10=y(:,9)+y(:,12)+y(:,13)+y(:,15);
w11=fft(y(:,3)+y(:,2)+y(:,4)+y(:,1)+y(:,7)+y(:,6)+y(:,8)+y(:,5)).*conj(fft(y(:,3)+y(:,2)+y(:,4)+y(:,1)+y(:,7)+y(:,6)+y(:,8)+y(:,5)));
w12=fft(y(:,3)+y(:,2)+y(:,4)+y(:,1)-(y(:,7)+y(:,6)+y(:,8)+y(:,5))).*conj(fft(y(:,3)+y(:,2)+y(:,4)+y(:,1)-(y(:,7)+y(:,6)+y(:,8)+y(:,5))));
w13=fft(y(:,11)+y(:,10)+y(:,12)+y(:,9)+y(:,15)+y(:,14)+y(:,16)+y(:,13)).*conj(fft(y(:,11)+y(:,10)+y(:,12)+y(:,9)+y(:,15)+y(:,14)+y(:,16)+y(:,13)));
w14=fft(y(:,11)+y(:,10)+y(:,12)+y(:,9)-(y(:,15)+y(:,14)+y(:,16)+y(:,13))).*conj(fft(y(:,11)+y(:,10)+y(:,12)+y(:,9)-(y(:,15)+y(:,14)+y(:,16)+y(:,13))));
w15=fft(y(:,17)).*conj(fft(y(:,17)));
w16=fft(y(:,18)).*conj(fft(y(:,18)));
w17=fft(y(:,19)).*conj(fft(y(:,19)));
w18=fft(y(:,20)).*conj(fft(y(:,20)));
%plot(w1);
%w1=fft(y(:,10)).*conj(fft(y(:,10)));
%plot(log10(xi),-log10(w2));
%plot(log(abs(fft(y(:,13)))));
%b(:)=trace(X(;) .* logm(X(:)))
%plot(log10(real(y(:,1)-y(:,4)+y(:,5)-y(:,8))),log10(imag(y(:,1)-y(:,4)+y(:,5)-y(:,8))));
%plot(log(abs(fft(y(:,13)))));
%plot(real(y(:,1)),imag(y(:,1)-y(:,4)+y(:,5)-y(:,8)));
%plot(abs(b(:)));
%plot3(y(:,9),real(y(:,1)+y(:,4)+y(:,5)+y(:,8)),imag(y(:,1)+y(:,4)+y(:,5)+y(:,8)));
N=N0./0.01
%for j=1:N
% wn1(j)=w1(j);
% wn2(j)=w2(j);
% wn3(j)=w3(j);
% wn4(j)=w4(j);
%end
n1=abs(sum(w1));
n2=abs(sum(w2));
n3=abs(sum(w3));
n4=abs(sum(w4));
n5=abs(sum(w5));
n6=abs(sum(w6));
n7=abs(sum(w7));
n8=abs(sum(w8));
n9=imag(sum(w9));
n10=imag(sum(w10));
n11=abs(sum(w11));
n12=abs(sum(w12));
n13=abs(sum(w13));
n14=abs(sum(w14));
n15=abs(sum(w15));
n16=abs(sum(w16));
n17=abs(sum(w17));
n18=abs(sum(w18));
%n5=abs(sum(wn5));
%tau1(l)=(2.*N)/(3.*pi) - (-4.*N^2 - 3.*pi^2)/(3.*pi.*(8.*N^3 - 27.*n1.*pi^2 + 9.*N.*pi^2 +3.*sqrt(3).*(-16.*n1.*N^3.*pi^2+ 27.*n1^2.*pi^4 - 18.*N.*n1.*pi^4 - N^2.*pi^4-pi^6)^(1/3)))+1./((8.*N^3-27.*n1.*pi^2+3.*3^(1/3).*(-16.*n1.*N^3.*pi^2+27.*n1^2.*pi^4-N^2.*pi^4-pi^6)^(1/3)).*(3.*pi));
tau1(l)=N./(n1-pi);
tau2(l)=N./(n2-pi);
tau3(l)=N./(n3-pi);
tau4(l)=N./(n4-pi);
tau5(l)=N./(n5-pi);
tau6(l)=N./(n6-pi);
tau7(l)=N./(n7-pi);
tau8(l)=N./(n8-pi);
tau9(l)=abs(n9);
tau10(l)=abs(n10);
tau11(l)=N./(n11-pi);
tau12(l)=N./(n12-pi);
tau13(l)=N./(n13-pi);
tau14(l)=N./(n14-pi);
tau15(l)=N./(n15-pi);
tau16(l)=N./(n16-pi);
tau17(l)=N./(n17-pi);
tau18(l)=N./(n18-pi);
%total1(l)=imag(sum(y(:,1)+y(:,4)+y(:,5)+y(:,8)));
%total2(l)=imag(sum(y(:,1)-y(:,4)+y(:,5)-y(:,8)));
%total3(l)=imag(sum(y(:,1)+y(:,4)-y(:,5)-y(:,8)));
%total4(l)=imag(sum(y(:,1)-y(:,4)-y(:,5)+y(:,8)));
%w2=(y(:,1)-1/tau1(l)).*conj((y(:,1)-1/tau1));
%plot(w4);
end
%plot(real(y(:,1)-y(:,4)+y(:,5)-y(:,8)),abs(y(:,17)),'c.')
%plot(real(y(:,1)-y(:,4)-y(:,5)+y(:,8)))
plot(abs(y(:,1)),abs(y(:,2)),'c.');
%plot(w15)
%plot(tau2)
%plot3(real(y(:,1)-y(:,4)+y(:,5)-y(:,8)),imag(y(:,1)-y(:,2)+y(:,5)-y(:,8)),abs(y(:,17)))