clear;
format compact;

%--minimally constrained array-----

%----the first value---

M=6;
signum=2;
derta=0.5;
N=1024;

theta1=30*pi/180;   theta111=30;
theta2=60*pi/180;   theta222=60;
fai1=20*pi/180;     fai111=20;
fai2=40*pi/180;     fai222=40;

gama1=45*pi/180;    gama111=45;
gama2=25*pi/180;    gama222=25;
jita1=5*pi/180;    jita111=5;
jita2=25*pi/180;    jita222=25;

f1=0.50;
f2=1.00;



%----s---
for SSS=1:11, 
   snr=5+(SSS-1)*2;
   Amp=sqrt(2*10^(snr/10));
 
for MM=1:50, 
    SSS
    MM
s1=(exp(j*2*pi*(f1*0.1*[0:N-1]+rand)))*Amp;   % fc/fs=0.7*0.1
s2=(exp(j*2*pi*(f2*0.1*[0:N-1]+rand)))*Amp;   % fc/fs=1.0*0.1

s=[s1(1:N);s2(1:N)];
%----a---

u1x=sin(gama1)*cos(theta1)*cos(fai1)*exp(j*jita1)-cos(gama1)*sin(fai1);  % x polarization element
u1y=sin(gama1)*cos(theta1)*sin(fai1)*exp(j*jita1)+cos(gama1)*cos(fai1);  % y polarization element
p1=exp(j*2*pi*derta*sin(theta1)*cos(fai1)*f1); % x delay
q1=exp(j*2*pi*derta*sin(theta1)*sin(fai1)*f1); % y delay
a1=[u1x;u1y;u1x*p1;u1x*q1;p1^2;q1^2];


u2x=sin(gama2)*cos(theta2)*cos(fai2)*exp(j*jita2)-cos(gama2)*sin(fai2);
u2y=sin(gama2)*cos(theta2)*sin(fai2)*exp(j*jita2)+cos(gama2)*cos(fai2);
p2=exp(j*2*pi*derta*sin(theta2)*cos(fai2)*f2);
q2=exp(j*2*pi*derta*sin(theta2)*sin(fai2)*f2);
a2=[u2x;u2y;u2x*p2;u2x*q2;p2^2;q2^2];

a=[a1,a2];
%----noise---

noise=randn(6,N)+j*randn(6,N);

%----z---

z=a*s+noise;

%----caculate c0,c1,c2,c3 c4----

for hang=1:6,
   for lie=1:6,
      zmat = z(1,:).* conj(z(1,:)).* z(hang,:);   
      R4 = (zmat * z(lie,:)') /N;      
      r12= (z(1,:) * z(1,:)') / N;
      r34= (z(hang,:)*z(lie,:)')/N;
      r13= (z(1,:) * z(hang,:).') / N;
      r24= (conj(z(1,:))*z(lie,:)')/N;
      r14= (z(1,:) * z(lie,:)') / N;
      r23= (conj(z(1,:))*z(hang,:).')/N;
      C0(hang,lie) = R4 - r12*r34-r13*r24-r14*r23; 
      
      zmat = z(1,:).* conj(z(3,:)).* z(hang,:);   
      R4 = (zmat * z(lie,:)') /N;      
      r12= (z(1,:) * z(3,:)') / N;
      r34= (z(hang,:)*z(lie,:)')/N;
      r13= (z(1,:) * z(hang,:).') / N;
      r24= (conj(z(3,:))*z(lie,:)')/N;
      r14= (z(1,:) * z(lie,:)') / N;
      r23= (conj(z(3,:))*z(hang,:).')/N;
      C1(hang,lie) = R4 - r12*r34-r13*r24-r14*r23; 

      zmat = z(1,:).* conj(z(4,:)).* z(hang,:);   
      R4 = (zmat * z(lie,:)') /N;      
      r12= (z(1,:) * z(4,:)') / N;
      r34= (z(hang,:)*z(lie,:)')/N;
      r13= (z(1,:) * z(hang,:).') / N;
      r24= (conj(z(4,:))*z(lie,:)')/N;
      r14= (z(1,:) * z(lie,:)') / N;
      r23= (conj(z(4,:))*z(hang,:).')/N;
      C2(hang,lie) = R4 - r12*r34-r13*r24-r14*r23; 
      
      zmat = z(1,:).* conj(z(2,:)).* z(hang,:);   
      R4 = (zmat * z(lie,:)') /N;      
      r12= (z(1,:) * z(2,:)') / N;
      r34= (z(hang,:)*z(lie,:)')/N;
      r13= (z(1,:) * z(hang,:).') / N;
      r24= (conj(z(2,:))*z(lie,:)')/N;
      r14= (z(1,:) * z(lie,:)') / N;
      r23= (conj(z(2,:))*z(hang,:).')/N;
      C3(hang,lie) = R4 - r12*r34-r13*r24-r14*r23;
      
      zmat = z(1,1:(N-1)).* conj(z(1,2:N)).* z(hang,1:(N-1));   %f
      R4 = (zmat * z(lie,1:(N-1))') /(N-1);      
      r12= (z(1,1:(N-1)) * z(1,2:N)') / (N-1);
      r34= (z(hang,1:(N-1))*z(lie,1:(N-1))')/(N-1);
      r13= (z(1,1:(N-1)) * z(hang,1:(N-1)).') / (N-1);
      r24= (conj(z(1,2:N))*z(lie,1:(N-1))')/(N-1);
      r14= (z(1,1:(N-1)) * z(lie,1:(N-1))') / (N-1);
      r23= (conj(z(1,2:N))*z(hang,1:(N-1)).')/(N-1);
      C4(hang,lie) = R4 - r12*r34-r13*r24-r14*r23; 
   end
end

%----estimate f, DOA and polarization----
[u s v]=svd(C0);
u1=u(:,1:signum);ss1=s(1:signum,1:signum);v1=v(:,1:signum);

Tc1=u1'*C1*v1;
[up faip]=eig(ss1,Tc1);

Tc2=u1'*C2*v1;
[uq faiq]=eig(ss1,Tc2);

Tc3=u1'*C3*v1;
[ur fair]=eig(ss1,Tc3);

Tc4=u1'*C4*v1;
[uf faif]=eig(ss1,Tc4);

c=[C0;C1;C2;C3;C4];
%----grouping----

up1=up(:,1)/norm(up(:,1));
up2=up(:,2)/norm(up(:,2));
uq1=uq(:,1)/norm(uq(:,1));
uq2=uq(:,2)/norm(uq(:,2));
ur1=ur(:,1)/norm(ur(:,1));
ur2=ur(:,2)/norm(ur(:,2));
uf1=uf(:,1)/norm(uf(:,1));
uf2=uf(:,2)/norm(uf(:,2));


faip1=faip(1,1);
faip2=faip(2,2);
if abs(up1'*uq1)<abs(up1'*uq2)
    faiq1=faiq(2,2);
    faiq2=faiq(1,1);
else
    faiq1=faiq(1,1);
    faiq2=faiq(2,2);
end

if abs(up1'*ur1)<abs(up1'*ur2)
    fair1=conj(fair(2,2));
    fair2=conj(fair(1,1));
else
    fair1=conj(fair(1,1));
    fair2=conj(fair(2,2));
end

if abs(up1'*uf1)<abs(up1'*uf2)
    faif1=faif(2,2);
    faif2=faif(1,1);
else
    faif1=faif(1,1);
    faif2=faif(2,2);
end
    
% if abs(faip1-p1)>abs(faip2-p1)
%   temp=faip2;
%   faip2=faip1;
%   faip1=temp;
%   temp=faiq2;
%   faiq2=faiq1;
%   faiq1=temp;
%   temp=fair2;
%   fair2=fair1;
%   fair1=temp;
%   temp=faif2;
%   faif2=faif1;
%   faif1=temp; 
% end

%----caculate the f, DOA and p----
f11(MM)=angle(faif1)*10/2/pi;
f22(MM)=angle(faif2)*10/2/pi;

theta11(MM)=asin((sqrt((angle(faip1)*angle(faip1))+(angle(faiq1)*angle(faiq1))))/pi/2/derta/f11(MM));
theta22(MM)=asin((sqrt((angle(faip2)*angle(faip2))+(angle(faiq2)*angle(faiq2))))/pi/2/derta/f22(MM));

fai11(MM)=atan(angle(faiq1)/angle(faip1));
fai22(MM)=atan(angle(faiq2)/angle(faip2));

ipsilon1=(fair1*cos(fai11(MM))+sin(fai11(MM)))/(cos(theta11(MM))*((-1)*fair1*sin(fai11(MM))+cos(fai11(MM))));
ipsilon2=(fair2*cos(fai22(MM))+sin(fai22(MM)))/(cos(theta22(MM))*((-1)*fair2*sin(fai22(MM))+cos(fai22(MM))));

jita11(MM)=angle(ipsilon1);
jita22(MM)=angle(ipsilon2);
gama11(MM)=atan(abs(ipsilon1));
gama22(MM)=atan(abs(ipsilon2));


theta11(MM)=theta11(MM)*180/pi;
theta22(MM)=theta22(MM)*180/pi;
fai11(MM)=fai11(MM)*180/pi;
fai22(MM)=fai22(MM)*180/pi;
gama11(MM)=gama11(MM)*180/pi;
gama22(MM)=gama22(MM)*180/pi;
jita11(MM)=jita11(MM)*180/pi;
jita22(MM)=jita22(MM)*180/pi;
   
if abs(theta11(MM)-theta111)>abs(theta22(MM)-theta111)
  temp=theta22(MM);
  theta22(MM)=theta11(MM);
  theta11(MM)=temp;
  temp=fai22(MM);
  fai22(MM)=fai11(MM);
  fai11(MM)=temp;
  temp=gama22(MM);
  gama22(MM)=gama11(MM);
  gama11(MM)=temp;
  temp=jita22(MM);
  jita22(MM)=jita11(MM);
  jita11(MM)=temp;
  temp=f22(MM);
  f22(MM)=f11(MM);
  f11(MM)=temp;
end
 
end

%----caculate the mean----
   mf1(SSS)=mean(f11);
   mf2(SSS)=mean(f22);
   mtheta1(SSS)=mean(theta11);
   mtheta2(SSS)=mean(theta22);
   mfai1(SSS)=mean(fai11);
   mfai2(SSS)=mean(fai22);
   mgama1(SSS)=mean(gama11);
   mgama2(SSS)=mean(gama22);
   mjita1(SSS)=mean(jita11);
   mjita2(SSS)=mean(jita22);
 
%----caculate the var---
   varf1(SSS)=std(f11);
   varf2(SSS)=std(f22);
   vart1(SSS)=std(theta11);
   vart2(SSS)=std(theta22);
   varfa1(SSS)=std(fai11);
   varfa2(SSS)=std(fai22);
   varg1(SSS)=std(gama11);
   varg2(SSS)=std(gama22);
   varj1(SSS)=std(jita11);
   varj2(SSS)=std(jita22);

end





%---crossed dipoles array-----

%----the first value---

L=3;
signum=2;
derta=0.5;
N=1024;

theta1=30*pi/180;
theta2=60*pi/180;
fai1=20*pi/180;
fai2=40*pi/180;

gama1=45*pi/180;
gama2=25*pi/180;
jita1=5*pi/180;
jita2=25*pi/180;

f1=0.50;
f2=1.00;

%----s---
for SSS=1:11, 
   snr=5+(SSS-1)*2;
   Amp=sqrt(2*10^(snr/10));
 
for MM=1:50,
    SSS
    MM
s1=(exp(j*2*pi*(f1*0.1*[0:N-1]+rand)))*Amp;
s2=(exp(j*2*pi*(f2*0.1*[0:N-1]+rand)))*Amp;

s=[s1(1:N);s2(1:N)];
%----a---

u1=[sin(gama1)*cos(theta1)*cos(fai1)*exp(j*jita1)-cos(gama1)*sin(fai1);
    sin(gama1)*cos(theta1)*sin(fai1)*exp(j*jita1)+cos(gama1)*cos(fai1)];
p1=exp(j*2*pi*derta*sin(theta1)*cos(fai1)*f1);   % *f1 
q1=exp(j*2*pi*derta*sin(theta1)*sin(fai1)*f1);   % *f1
a1=[u1;u1*p1;u1*(p1^2);u1*q1;u1*(q1^2)];


u2=[sin(gama2)*cos(theta2)*cos(fai2)*exp(j*jita2)-cos(gama2)*sin(fai2);
    sin(gama2)*cos(theta2)*sin(fai2)*exp(j*jita2)+cos(gama2)*cos(fai2)];
p2=exp(j*2*pi*derta*sin(theta2)*cos(fai2)*f2);
q2=exp(j*2*pi*derta*sin(theta2)*sin(fai2)*f2);
a2=[u2;u2*p2;u2*(p2^2);u2*q2;u2*(q2^2)];