% SCRIPT:  Beamforming with an antenna array of omnidirectional elements.
%
% An increasing number of elements is placed randomly in the x-y plane.
%
% The response is normalized to 1 at direction t0=0,
% and its max. magnitude minimized outside the beam.
%
% The optimization is done for a beamwidth that is inversely proportional
% to the number of elements.
%
% Plots the achieved sidelobe level as a fct of the number of elements.
%

Nu=20;		% SOCP parameter

ma=20;		% n. of angle samples per array element
mb=0;		% constant n. of angle samples

S=1;		% std of element x-y positions (note that lambda=2*pi)

n0=4;		% run from n0 elements (n0>=2) to nmax elements
nmax=32;	% XXX for nmax larger than 15, seems to need larger ma factor
		% to stay well conditioned (lapack is complaining)
	% initial array positions
px=random('norm',0,S,n0-1,1);
py=random('norm',0,S,n0-1,1);
	% initialize lists to record results
l_SLL=[];
l_iter=[];

for n=n0:nmax,
	% beamwidth= 30;
	beamwidth= 2*360/n;

		% angles: center, left and right of beam
	t0=0;
	t1=beamwidth/2/180*pi;
	t2=(360-beamwidth/2)/180*pi;

		% add in another element to the array
	px=[px; random('norm',0,S,1,1)];
	py=[py; random('norm',0,S,1,1)];
		% draw array
	figure(1);
	plot(px,py,'o');
	xlabel('x');
	ylabel('y');
	title('position of array elements');

		% number of samples for angle discretization
	m=ma*n+mb;
	dt=(t2-t1)/(m-1);

		% compute response at discretization angles
		% (i.e. y=A*w, with complex values)
	A=[];
	C=[];
	for t=t1:dt:(t2+eps),
		a=[];
		for j=1:n,
			gamma = px(j)*cos(t)+py(j)*sin(t);
			a_re  = cos(gamma);
			a_im  = sin(gamma);
			a=[a, [a_re,-a_im; a_im,a_re]];
		end;
		A=[A; [a,[0;0]]];
		C=[C; [zeros(1,2*n),1]];
	end;
	b=zeros(2*m,1);
	d=zeros(m,1);
	N=2*ones(m,1);
	f=[zeros(2*n,1); 1];	% minimize SLL

		% find response at theta=0  (i.e. y0=a0*w, with complex vals.)
	t=t0;
	a0=[];
	for j=1:n,
		gamma = px(j)*cos(t)+py(j)*sin(t);
		a_re  = cos(gamma);
		a_im  = sin(gamma);
		a0=[a0, [a_re,-a_im; a_im,a_re]];
	end;
		% compute initial point satifying beam=1 at theta=0
	w0=a0\[1;0];
	y0=A*[w0;0];
	Ire=1:2:2*m;
	Iim=Ire+1;
	SLL0=sqrt(max(y0(Ire).^2+y0(Iim).^2));
	x0=[w0; 1.1*SLL0];
		% normalize problem to enforce beam=1 at theta=0
	a1=[a0,[0;0]];
	T=null(a1);
	v=a1\[1;0];
	f1=T'*f;
	A1=A*T;
	C1=C*T;
	b1=b+A*v;
	d1=d+C*v;
	x1=T\(x0-v);

		% solve  (can also use B0=2*m*SLL0)
	[x2,info,z2,w2,hist,time]= ...
			socp(f1,A1,b1,C1,d1,N,x1,[],[],0,1e-3,0,100,Nu);

		% undo normalization
	x=T*x2+v;

		% extract SLL and weights from problem variable
	SLL=x(2*n+1);
	w=x(1:2*n);

		% update lists
	l_iter=[l_iter; time(3)];
	l_SLL=[l_SLL; SLL];

		% plot number of iterations vs. n
	l_n=n0:n;
	figure(2);
	plot(l_n,l_iter,'o');
	axis([0 n 0 max(l_iter)+1]);
	xlabel('n');
	ylabel('i');
	title('SOCP iterations');

		% plot optimal sidelobe level vs. n
	figure(3);
	plot(l_n,20*log10(l_SLL));
	xlabel('n');
	ylabel('dB');
	title('SLL');

		% compute and plot angular response with twice the resolution
	G=[];
	for t=0:dt/2:2*pi,
		y_re=0;
		y_im=0;
		for j=1:n,
			gamma = px(j)*cos(t)+py(j)*sin(t);
			a_re  = cos(gamma);
			a_im  = sin(gamma);
			y_re=y_re+a_re*w(2*j-1)-a_im*w(2*j);
			y_im=y_im+a_im*w(2*j-1)+a_re*w(2*j);
		end;
 		G=[G; sqrt(y_re^2+y_im^2)];
	end;
	figure(4);
	t=(0:dt/2:2*pi)';
	G1=20*log10(G);
	SLL1=20*log10(SLL);
	dG1=2*SLL1;
	G1=max(0,G1-dG1);
	SLL1=SLL1-dG1;
	plot(cos(t).*G1,sin(t).*G1,'-');
	axis('equal');
	axis('off');
	hold on;
	t=(t1:dt/2:t2)';
	plot(cos(t)*SLL1,sin(t)*SLL1,'--');
	plot([cos(t1)*SLL1 cos(t1)*2*SLL1],[sin(t1)*SLL1 sin(t1)*2*SLL1],'--');
	plot([cos(t2)*SLL1 cos(t2)*2*SLL1],[sin(t2)*SLL1 sin(t2)*2*SLL1],'--');
	t=[t2:dt/2:2*pi 0:dt/2:t1]';
	plot(cos(t)*2*SLL1,sin(t)*2*SLL1,'--');
	hold off;

	%%% save array-xy.mat
end;
