function planearray_jitter

% Inspect the effects of the elements' amp and phase jitter to the
% array pattern. We inspect this only in the xz-plain only, and assume 
% no phase steeting in the y-direction. 
%
%    y  ^
%       |
%       *   *   *   *
%       |                   M_y elements in y-direction
%    0  *--*--*--*--*--> x
%       0
%                           M_x elements in x-direction
%

% 11-Oct-2005  Jm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(0,'defaultaxesfontsize',14);
rad = pi/180;
deg = 180/pi;

%=======================% CHANGE THESE INPUT PARAMETERS AS APPROPRIATE

%randn('state',0);      % Comment this out if new set of jitter is required
                        % for every run.

fradar = 225e6;         % Radar freq in Hz.  

Mx = 10;                % Number of elements x direction.
My = 10;                % Number of elements in y direction.
Dx = 1.50;              % Element distance divided by wavelenth.

sigma_tim = 500E-12;    % Sdev (seconds) of timing jitter at each element.
sigma_amp = 0.20;       % Sdev of amplitude jitter, assuming amplitude = 1.
                        % Either (or both) of these can also be zero.

delta = -90.0*rad;      % Phasing step between the elements in radians
                        % (in the x-direction only).
                        % Should be between -pi and pi.

Ntrials = 1000;         % Number of trials in the Monte-Carlo run.  
show_movie = 0;         % If set, show the AF's in a separate window              
                        % as they get computed. This slows things down a lot.


%====================== NORMALLY too not touch anything below this line


if Mx < 1 || My < 1 || sigma_amp < 0
    error ('Invalid parameter')
end

lightcyan = [0.8 1 1];
blue = [ 0 0 1];
lightgray = [1 1 1]*0.85;
lightblue = [0.8 0.9 1];

% The main beam lobe and nearest few sidelobes on both side if it, 
% when there is no phase errors
%
%   sin theta(+-) = -delta +-2*2pi/M

n = 2.5; 
x1 = (-delta - n*2*pi/Mx) / (2*pi*Dx );      
x2 = (-delta + n*2*pi/Mx) / (2*pi*Dx );       
x0 = -delta/(2*pi*Dx);

theta1 = asin(x1);                      % Direction of null
theta2 = asin(x2);                      % Direction of null
theta0 = asin(x0);                      % Direction of main maximum 

theta = linspace (theta1,theta2,2000);   % Zenith-angle grid
N = length(theta);

% Compute all the terms going into the array factor sum, for all m and theta.
% This is common in all trials, so we only compute it once.
%
%   Z0(m,n) = exp { (2piDsin(theta(n)) + delta)m }
%

Z0 = zeros(Mx,N);

Z0(1,:) = exp (i * (2*pi*Dx*sin(theta) + delta) );
for m = 2:Mx
    Z0(m,:) = Z0(m-1,:) .* Z0(1,:);
end

% The un-distorted array factor is

AF0 = abs(sum ( Z0, 1))/Mx;

% ===================================================================
% A Monte-Carlo with random phase error added to each element's phase
% ===================================================================

% Generate the timing jitter and compute the corresponding 
% phase jitter, for all trials. 
%
% We assume the y-direction to be trivial, so that we can sum
% the jitter factors in the y-direction, to reduce the situation to 
% the one-dimensional case (see linearray_jitter.m).
%
%                                                   
%   Jitter(trial,mx) = sum_my {amp(trial,mx,my) * exp[i phi(trial,mx,my)]}

Jitter = zeros (Ntrials,Mx);
Amp = zeros(Ntrials*Mx,1);
Pha = Amp;
k1 = 1;
k2 = My*Mx;
for trial = 1:Ntrials

    % Time and phase jitter at the elements.

    tim = sigma_tim*randn(My,Mx);
    pha = 2*pi*fradar*tim;

    % Amplitude jitter. Normalize to unit total power.

    amp = 1 + sigma_amp*randn(My,Mx);

    normalizer = 1/sum(sum (amp.^2, 1));
    amp = normalizer * amp;

    % Jitter summed over y-direction

    Jitter(trial,:) = sum (amp.*exp(i*pha), 1);

    % We want to plot the tim and amp jitter distributions also,
    % so we accumulate the distribution here.

    Amp(k1:k2) = amp;
    Pha(k1:k2) = pha;
    k1 = k1+Mx;
    k2 = k2+Mx;
end

% Initialize result arrays

max_AF = zeros (Ntrials,1);
dir_AF = zeros (Ntrials,1);
beamwidth_AF = zeros (Ntrials,1);

% Find the half power beam width of the undistorted beam.

[maxAF0, k0, k1, k2] = beamparam (AF0,sqrt(0.5));
beamwidth_AF0 = theta(k2) - theta(k1);

% Prepare to plot only every NN'th AF so that the total number
% plottes is is about 20.

NN = ceil(Ntrials/20);

% Open the jitter figure and its panel for the AFs.

getfig ('plane-array jitter (1)', [0.25 2.5 7 8]);
subplot ('position',[0.1300 0.1100 0.7750 0.32]);

% The next loops computes and plots the phase and amp jitter-affected
% array functions around the beam main direction, 
% and collects AF shape jitter statistics.

for trial = 1:Ntrials

    % The array factor, as a function of theta. 
    % (Jitter(trial,:) is 1xM, Z0 is MxN, so AF is 1xN.

    AF = abs (Jitter(trial,:)*Z0);

    if rem(trial,NN) == 1
        plot (theta*deg,AF,'c-'); 
        hold on;
    end

    % Find the half power beam width etc of the distorted beam.

    [maxAF, k, k1, k2] = beamparam (AF,sqrt(0.5));
    max_AF(trial) = maxAF;
    dir_AF(trial) = theta(k);
    beamwidth_AF(trial) = theta(k2) - theta(k1);

    if show_movie

        % Show the distorted AF on a separate window,
        % together with the undistored one.

        getfig ('plane-array jitter (3)')
        plot (theta*deg,AF0,'r--', theta*deg, AF, 'b-')
        set (gca,'tickdir','out');
        set (gca,'ylim',[0 1.1]);
        set (gca,'xlim',[theta1, theta2]*deg);
        xlabel ('Zenith angle (\circ)');
        ylabel ('AF');
        pause(0.20)
%        pause
        getfig ('plane-array jitter (1)');
        subplot ('position',[0.1300 0.1100 0.7750 0.32]);
    end

end %trials

h = plot (theta*deg,AF0,'r-'); 
hold off;
set(h,'linewidth',1.5);
set (gca,'tickdir','out');
set (gca,'ylim',[0,1.1]);
set (gca,'layer','top');
xlabel ('Zenith angle \theta [\circ]');
ylabel ('AF');
h = title('Jitter of the antenna array factor');
set(h,'fontweight','bold');


% Find and plot the changes in the main beam's shape parameters
% due to timing jitter.

% Jitter in AFmax 
% ===============

err_maxAF = max_AF - 1;
mean_errmaxAF = mean (err_maxAF);
sdev_errmaxAF = std (err_maxAF);

Nbins = 15;
X = err_maxAF*100;
x_mean = mean_errmaxAF*100;
x_sigma = sdev_errmaxAF*100;

x1 = x_mean -4*x_sigma;
x2 = x_mean +3*x_sigma;

if x_sigma > 1
    Xc = linspace (x1,x2,Nbins);
else
    x1 = floor(x_mean)-3.5;
    x2 = ceil(x_mean+3.5);
    Xc = (x1-0.5):0.25:(x2+0.5); 
end
Nx = hist (X,Xc)/length(X)*100;

getfig ('plane-array jitter (2)',[0.25 2.5 7 6]);

subplot('position',[0.12 0.50, 0.33, 0.2]);     % <--- Fig (2), pane 1
    h = jstairs (Xc, Nx, lightcyan, blue);
    set (gca,'tickdir','out');
    set (gca,'ticklen',[0.02 0.05]);
    set (gca,'layer','top');
    set (gca,'xlim',[x1 x2]);
    set (gca,'ylim',[0 max(Nx)*1.1]);
    xlabel ('\Delta AF_{max} [%]')
    ylabel ('% / bin');
    h = title('Beam amplitude jitter');
    set(h,'fontweight','bold');

% Beam direction jitter
% =====================

err_dirAF = theta0 - dir_AF;
mean_err_dirAF = deg * mean (err_dirAF);
sdev_err_dirAF = deg * std (err_dirAF);

X = err_dirAF*deg;
x_mean = mean_err_dirAF;
x_sigma = sdev_err_dirAF;

Nbins = 15;
if x_sigma > 0.02
    x1 = -3*x_sigma;
    x2 = +3*x_sigma;
    Xc = linspace (x1,x2,Nbins);
else
    x1 = -0.055;
    x2 = +0.055;
    Xc = -0.05:0.01:0.05;
end
Nx = hist (X,Xc)/length(X)*100;

subplot('position',[0.57 0.50, 0.33, 0.2]);      % <--- Fig (2), pane 2
    h = jstairs (Xc, Nx, lightcyan, blue);
    set (gca,'tickdir','out');
    set(gca,'ticklen',[0.02 0.05]);
    set(gca,'layer','top');
    set (gca,'xlim',[x1 x2]);
    set (gca,'ylim',[0 1.1*max(Nx)]);
    xlabel ('\Delta Dir [\circ]')
    ylabel ('% / bin');
    h = title('Beam direction jitter');
    set(h,'fontweight','bold');

% Beam width jitter
% =================

err_bwidth = (beamwidth_AF0 - beamwidth_AF)*deg;
mean_err_bwidth = mean (err_bwidth);
sdev_err_bwidth = std (err_bwidth);

X = err_bwidth;
x_mean = mean_err_bwidth;
x_sigma = sdev_err_bwidth;

if x_sigma > 0.02
    Nbins = 15;
    x1 = -3*x_sigma;
    x2 = +3*x_sigma;
    Xc = linspace (x1,x2,Nbins);
else
    x1 = -0.055;
    x2 = +0.055;
    Xc = -0.05:0.01:0.05;
end
Nx = hist (X,Xc)/length(X)*100;

subplot('position',[0.12 0.12, 0.33, 0.2]);     % <--- Fig (2), pane 3
    h = jstairs (Xc, Nx, lightcyan, blue);
    set (gca,'tickdir','out');
    set (gca,'ticklen',[0.02 0.05]);
    set (gca,'layer','top');
    set (gca,'xlim',[x1 x2]);
    set (gca,'ylim',[0 1.05*max(Nx)]);
    xlabel ('\Delta BW [\circ]');
    ylabel ('% / bin');
    h = title('Beam width jitter');
    set(h,'fontweight','bold');

getfig ('plane-array jitter (1)');              % <=== Fig 1

% Actual phase jitter distribution
% ================================

Nbins = 21;
sigma_pha = 2*pi*fradar*sigma_tim;
pha1 = -3*sigma_pha*deg;
pha2 = +3*sigma_pha*deg;
sigma_pha_actual = std (Pha);
if sigma_pha_actual*deg > 1.5
    pha_c = linspace (pha1,pha2,Nbins);
else
    pha1 = -5;
    pha2 = 5;
    pha_c = -5:0.2:5;
end
N_pha = 100 * hist (Pha*deg,pha_c)/length(Pha);

subplot('position',[0.57 0.55 0.3347 0.23]);    % <--- Fig (1), pane 2
    h = jstairs (pha_c, N_pha, lightblue, blue);
    set (gca,'tickdir','out');
    set (gca,'ticklen',[0.02 0.05]);
    set (gca,'xlim',1.1*[pha1,pha2]);
    set (gca,'ylim',1.1*[0,max(N_pha)]);
    set (gca,'layer','top');
    xlabel ('{\Delta}\phi [\circ]');
    ylabel ('% / bin');
    h = title('Phase jitter');
    set(h,'fontweight','bold');


% Amplitude jitter
% ================

Amp = Mx*My*Amp;
 
sigma_amp = std (Amp);
ave_amp = mean (Amp);
amp1 = ave_amp - 3*sigma_amp;
amp2 = ave_amp + +3*sigma_amp;
sigma_amp_pcent = 100*sigma_amp/ave_amp;
if sigma_amp_pcent < 5
    amp_c = 0.92:0.02:1.08;
else
    amp_c = linspace (amp1,amp2,Nbins);
end
N_amp = 100 * hist (Amp, amp_c) / length (Amp);

subplot('position',[0.13 0.55 0.3347 0.23]);    % <--- Fig (1), pane 1
    h = jstairs (amp_c, N_amp, lightblue, blue);
    set (gca,'tickdir','out');
    set (gca,'ticklen',[0.02 0.05]);
    set (gca,'layer','top');
    set (gca,'xlim',[0,1.7]);
    set (gca,'ylim',1.1*[0,max(N_amp)]);
    xlabel ('A');
    ylabel ('% / bin');
    h = title('Amplitude jitter');
    set(h,'fontweight','bold');

% Figure (1) header
% =================

h1 = subplot('position',[0.05,0.85,0.8550,0.10]);   % <--- Fig (1), header

    fmt1 = 'Mx=%d My=%d d/\\lambda=%1.2f';  
      fmt2 = '  Steering \\delta=%1.1f\\circ'; 
        fmt3 = '  Dir=%1.2f\\circ  BW=%1.2f\\circ';
    hdr{1} = sprintf ([fmt1 fmt2 fmt3], ...
                Mx, My, Dx, ...
                    -delta*deg, ... 
                        theta0*deg, beamwidth_AF0*deg);

    fmt1 = '{\\Delta}t=%1.0f ps  \\Delta\\phi=%1.1f\\circ'; 
      fmt2 = '  {\\Delta}A=%1.1f%%   N=%d';
    hdr{2} = sprintf ([fmt1 fmt2], ...
                  sigma_tim*1E12, sigma_pha_actual*deg, ...
                     sigma_amp_pcent , Ntrials );

    fmt1 = '{\\Delta}AF=%1.1f%%{\\pm}%1.1f';
      fmt2 = '  {\\Delta}Dir=%1.3f\\circ{\\pm}%1.3f';
        fmt3 = '  {\\Delta}BW=%1.3f\\circ{\\pm}%1.3f';
    hdr{3} = sprintf ( [fmt1 fmt2 fmt3], ...
                  mean_errmaxAF*100, sdev_errmaxAF*100, ...
                     mean_err_dirAF, sdev_err_dirAF, ...
                        mean_err_bwidth, sdev_err_bwidth ...
            );

    header (h1, hdr);


% Copy Fig (1) header to Fig (2)
% ==============================

getfig('plane-array jitter (2)');    
    h1 = subplot('position',[0.05,0.80,0.8550,0.13]);
    header (h1, hdr);


%end linearray_jitter


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%               SUBFUNCTIONS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function hh = getfig (Figurename, Paperpos)
%----------------------------------------
% Create a figure window if one does not already exists,
% otherwise just select the named figure.
%----------------------------------------

Width = 500;

if nargin == 0
    error ('Wrong number of arguments')
end 

h = findobj('Tag',Figurename);
if isempty(h)
    h = figure;
    set(h,'Tag',Figurename);
    set(h,'Numbertitle','off');
    colordef(h,'white');
    set(h,'PaperType','A4');
    set(h,'Name',Figurename);
    if nargin == 2
        Height = (Paperpos(4)/Paperpos(3))*Width;
        pos = get(gcf,'position');
        set(h,'Position',[pos(1), pos(2), Width, Height]);
        set(gcf,'paperpos',Paperpos);
    end
else
    figure(h);
end

if nargout > 0
    hh = h;
end    


function header (h,txt)
%============================================================
% 
%============================================================
cla
Eventnum = [];
set(h,  'xticklabel',[],...
        'yticklabel',[],...
        'xtick',[],...
        'ytick',[],...
        'box','on',...
        'color',[1 1 1]);

        
N = length(txt);        
n = N:-1:1;
y = (2*n-1)/(2*N);
for k = 1:N
    hh = text(0.02,y(k),txt{k});
    set(hh,...
        'Color',[0 0 0],...
        'VerticalAlignment','middle'...
        );
        set(hh,'fontsize',14);

end    

if ~isempty(Eventnum)

    pos = get(h,'pos'); 

    wd = 0.07;
    hgt = 0.03;
    x1 = pos(1)+pos(3) - wd;
    y1 = pos(2)+pos(4) - hgt;

    h = axes('position',[x1,y1,wd,hgt]);
    set(h,  'xticklabel',[],...
            'yticklabel',[],...
            'xtick',[],...
            'ytick',[],...
            'box','on',...
            'color',[1 1 1]);
    h = text(0.5,0.5,num2str(Eventnum));
    set(h,'horizontalAlignment','center');
    set(h,'verticalAlignment','middle');
end

%header


function h = jstairs (xc, N, facecolor, edgecolor)
%-----------------------------------------------------
% Staircase plot which centers the steps at xc, and the
% value of the step level is N. Basically, this is meant
% to draw a "histogram without the vertical lines".
%-----------------------------------------------------

dx = xc(2)-xc(1);
xc2 = [xc(:);xc(end)+dx];

[xx,yy] = stairs ( xc2-dx/2,[N(:);N(end)] );
h = area(xx,yy);
set(h,'facecolor',facecolor,'edgecolor',edgecolor);


function jsubplot (n,m,p, hscale)
%-----------------------------------------------------
% Vertically squeezed subplot window
%-----------------------------------------------------
h = subplot(n,m,p);
x = get(h,'pos');
set(h,'pos',[x(1) x(2) x(3) hscale*x(4)]);


function [ymax, k, k1, k2] = beamparam (y, ref)
%---------------------------------------------------------------
% Maximum value, and the indexes k1 and k2 of the nearest the points
% P1 and P2 both sides of the maximum  nearest to the max where
% y gets lower than the given ref level. The ref level is relative
% to ymax, so that ref = 1.0 means ymax.
%---------------------------------------------------------------

[ymax, k] = max (y);
m1 = min (find (y(k:-1:1) <= ref*ymax));
m2 = min (find (y(k:end) <= ref*ymax));
k1 = k - m1 + 1;
k2 = k + m2 - 1;
