theta0=30; % угол места
phi0=30;  % азимут

N=32;
Nxr=N;
Nyr=N;
dolx=0.5;
doly=0.5;

eps = 0.001;
nx = 0:Nxr-1;
ny = 0:Nyr-1;
i = sqrt(-1);

% choose proper size fft
Nrx = 10 * Nxr;
Nry = 10 * Nyr;
nfftx = 2^(ceil(log(Nrx)/log(2)));
nffty = 2^(ceil(log(Nry)/log(2)));

% generate array of elements
array = ones(Nxr,Nyr);

%======================================================================

% convert steering angles (theta0, phi0) to radians
theta0 = theta0 * pi / 180;
phi0 = phi0 * pi / 180;
% convert steering angles (theta0, phi0) to U-V sine-space
u0 = sin(theta0) * cos(phi0);
v0 = sin(theta0) * sin(phi0);

% Use formula thetal = (2*pi*n*dol) * sin(theta0)
    phasem = exp(i*2*pi*dolx*u0 .* nx );
    phasen = exp(i*2*pi*doly*v0 .* ny );

% add the phase shift terms
array = array .* (transpose(phasem) * phasen);
w = array;  

%==========================================================================

% Compute array pattern
arrayfft = abs(fftshift(fft2(w,nfftx,nffty))).^2;

%compute [su,sv] matrix
U = [-nfftx/2:(nfftx/2)-1] ./(dolx*nfftx);
indexx = find(abs(U) <= 1);
U = U(indexx);

V = [-nffty/2:(nffty/2)-1] ./(doly*nffty);
indexy = find(abs(V) <= 1);
V = V(indexy);

%Normalize to generate gain patern
rbar=sum(sum(arrayfft(indexx,indexy))) / dolx/doly/4./nfftx/nffty;
arrayfft = arrayfft(indexx,indexy) ./rbar;

V = asin(V);
U = asin(U);
V = V .* (180.0 / pi);
U = U .* (180.0 / pi);

pattern = 10*log10(arrayfft +eps);
figure(1)
mesh(V,U,pattern);