% 2-D glint simulation 
clear
rad=pi/180;
k=2*pi;
% amplitudes and phases of the scatterers (rectangular form)
% A=[1 1+j 1+1j 1 2+j];
  A=[1 1 0 0 0];
a=abs(A); phi=angle(A);
Nscat=length(A);
W=2;  % wingspan
L=4;  % nose to tail
% coordinates of the scatterers
x=[-W/2 W/2 0 0 0];  % first two wingtips, second two nose to tail
y=[0 0 L/2 -L/2 0];
% x=[0]; y=[0];    % point source at origin
iplt=0;
if iplt==0;
figure(1),clf
plot(x(1:2),y(1:2),'-o'),hold on, plot(x(3:4),y(3:4),'-x')
D=max([max(x),max(y)]); D1=D+0.2*D;
axis([-D1,D1,-D1,D1])
axis square
xlabel('x')
ylabel('y')
title('scatterer locations in z=0 plane')
end  % end of plot

icent=0;
if icent==0
% generate grid of observation point coordinates in x-y plane
X=-30:.1:30; Y=X; Nx=length(X); Ny=length(Y);
[u,v]=meshgrid(X,Y);
% compute the total scattered field at each observation point
for n=1:Nx
    for m=1:Ny
        if n<4 & m<4, E(m,n)=NaN; P(m,n)=NaN; 
        else
        summ=0;
        for i=1:Nscat
            R=sqrt((u(m,n)-x(i))^2+(v(m,n)-y(i))^2);
            summ=summ+a(i)*exp(j*phi(i))*exp(-j*k*R)/R; 
        end
        E(m,n)=abs(summ);
        P(m,n)=angle(summ);
    end  % end of if
end   % end of m loop
end  % end of n loop
Edb=20*log10(E);
Q=P; %unwrap(P);
figure(2)
contour(u,v,Edb)
axis square
title('Amplitude contours, dB')
xlabel('x, wavelengths'),ylabel('y, wavelengths')
figure(3)
contour(u,v,Q)
axis square
title('Phase contours, radians')
xlabel('x, wavelengths'),ylabel('y, wavelengths')
end % end of icent
ifar=1;
if ifar==0;
% compute contours in a box at large range
X=10000:.05:10020; Y=10000:.05:10020; Nx=length(X); Ny=length(Y);
% X=-10:.05:10; Y=-10:.05:10; Nx=length(X); Ny=length(Y);
[u,v]=meshgrid(X,Y);
% compute the total scattered field at each observation point
for n=1:Nx
    for m=1:Ny
    summ=0;
    for i=1:Nscat
        R=sqrt((u(m,n)-x(i))^2+(v(m,n)-y(i))^2);
        summ=summ+a(i)*exp(j*phi(i))*exp(-j*k*R)/R; 
    end
    E(m,n)=abs(summ);
    P(m,n)=unwrap(angle(summ));
    end
end
Edb=20*log10(E);
% Q=unwrap(P);
% take out phase to the center of the box
%R0=sqrt((10150)^2+(10150)^2);
%P0=-k*R0;
%Q=Q-P0;
Q=P;
figure(4)
contour(u,v,Edb)
axis square
title('Amplitude contours, dB')
xlabel('x, wavelengths'),ylabel('y, wavelengths')
figure(5)
contour(u,v,Q,1)
axis square
title('Phase contours, rad')
xlabel('x, wavelengths'),ylabel('y, wavelengths')
end
figure(6)
clf
% compute a polar plot at max range
R=max(X);
dt=0.5;
it=0;
for ang=0:dt:360
    U=cos(ang*rad);
    V=sin(ang*rad);
    it=it+1;
    summ=0;
        for i=1:Nscat
            summ=summ+a(i)*exp(j*k*(x(i)*U+y(i)*V)); 
        end
        Eff(it)=abs(summ);
        Ang(it)=ang;
end  % end of n loop
Edb=20*log10(Eff/max(Eff));
polardb(Ang,Edb,'k-')