%% On the Use of Convex Optimization for Electromagnetic Near-field Shaping
%  AJ Malcolm
%% The Paper
% 
%  The paper discusses the use of convex optimization to synthesize an
%  electromagnetic field on an aperture (Eap) such that the near field (Enf) at
%  position x = 0, y = 0, z = z0 is maximized, and constrained around the
%  maximum value to a max sidelobe level of rho. The field is calculated with:
% 
% $$E_{NF}(x,y,z_0)= \frac{1}{4\pi^2} F^{-1}(gF(E_{AP}(x,y,z=0)))$$
% 
%  where 
%
% $F, F^{-1}, g = e^{-jk_zz}$ 
%
%  are the direct, inverse two-dimensional Fourier transform, and  propagation 
%  factor in z respectively. This paper builds upon a previous paper by the
%  same authors where they use an iterative algorithm to slowly optimize
%  the near-field shape based on some initial guess. They propose the
%  convex optimization approach as it does not require the initial guess.
%  
%  Based on suggested parameters in the paper, a simplification of the
%  aforementioned equation can be made. The height z0 is given as
%  10*lambda, while the wavenumber k_z is equal to (2*pi/lambda). Plugging
%  this into the progagation factor yields an exponent power of (j*20*pi),
%  making the term equal to 1. This means the inverse and direct Fourier
%  transforms cancel each other out, with the result that:
% 
% $$E_{NF}(x,y,z_0)= \frac{1}{4\pi^2} (E_{AP}(x,y,z=0))$$
% 
%  Thus, maximizing the near-field is equivalent to maximizing the aperture
%  field. The issue lies in this fact, as although there are constraints
%  specified on the near-field, there are none on the aperture field. Thus,
%  the solution to maximizing the field is a point-source of some arbitrary
%  value at x=0, y=0. Obviously, this is the trivial solution.
%  
%  A common aperture field shape would be that of a Bessel function, defined in 
%  polar coordinates. However, this is not a good constraint to use with the 
%  aperture for two reasons: one, CVX doesn't support the Matlab Bessel function, 
%  and two, the paper demonstrates that the aperture field is not constrained 
%  by a Bessel function, as the value at r is not constant around theta.



%% Defining the Problem
%The problem is to maximize the near field, Enf, at x=0, y=0, z=z0. The
%near field itself is a function of the aperture field, Eap. Thus, the
%method of approach would be to maximize over Eap as the variable, because
%f(Eap) = Enf. Then constraints would be placed on f to constrain Enf. Due
%to the linear dependence between them, maximizing over Eap would be
%equivalent to maximizing Enf.

%;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% problem specifications
%;;;;;;;;;;;;;;;;;;;;;;;;;;;;
n = 33;                          % Matrix size
lambda = 1;                     % Wavelength 
alpha = 10*lambda;              % Aperture diameter
BW = 1.12*lambda;               % Beam width
kz = 2*pi/lambda;               % Wavenumber of z

prop_factor = @(z) exp(-1j*kz*z); % Propagation factor calculation for z

r = round(n/2);                 % Radius of largest circle in Matrix

%;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Field and constraint domains
%;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Defining the domain of Eap as a circle with radius r
D_Eap = zeros(n);               
for i = 1:n
    for j = 1:n
        if norm([i j]-[r r]) <= r
            D_Eap(i,j) = 1;     % 1 if x,y in the domain 
        else
            D_Eap(i,j) = 0;     % 0 if x,y not in the domain
        end
    end
end

% Defining the domain of Omega_SLL as a plane with circular opening with 
% radius BW
D_Osll = zeros(n);               
for i = 1:n
    for j = 1:n
        if norm([i j]-[r r]) >= BW
            D_Osll(i,j) = 1;     % 1 if x,y in the domain 
        else
            D_Osll(i,j) = 0;     % 0 if x,y not in the domain
        end
    end
end

    

%% Different Attempts to Solve the Problem
%;;;;;;;;;;;;;;;;;;;
% Constraint Values
%;;;;;;;;;;;;;;;;;;;
rho = 15*D_Osll;
g = prop_factor(alpha);

%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Fast Fourier Transform
%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Initial attempt to solve the problem. This does not work as CVX does not
% support the FFT or IFFT.
%{
cvx_begin
    variable Eap(n,n) complex
    variable Enf(n,n) 
    maximize (1/4*pi^2)*ifft2(g*fft2(Eap(r,r)))
    subject to
        abs(Enf) <= rho
cvx_end
%}

%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Bessel Function Objective
%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Attempt to solve the problem by starting with optimizing over a Bessel
% function before applying it as a constraint to the aperture field. CVX
% does not support the Bessel function, so it cannot be used as a
% constraint.
%{
[X Y] = meshgrid(-1:0.01:1,-1:0.01:1);
sqrt(X.^2+Y.^2);
cvx_begin
    variable a
    maximize max((2*besselj(1,a*R(:))./R(:)).^2)
cvx_end
%}

%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Fourier independent along z
%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
%An attempt to model the field along z = 0 to z = z0, instead of just at
%those two plans. As the Fourier transforms are on (x,y), the z-dependent
%propagation factor is independent of them. Thus, it can be extracted and
%the inverse and discrete transforms cancel. 
[X Y] = meshgrid(-1:0.01:1,-1:0.01:1);
sqrt(X.^2+Y.^2);
cvx_begin
    variable z
    variable Eap(n,n)
    maximize abs(-1j*kz*z)*log(Eap(r,r))
    subject to
        prop_factor(z)*Eap <= rho
cvx_end

%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
% Discrete Fourier Transform
%;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
%CVX does not support the FFT or IFFT functions, thus it would be necessary
%to implement them in Matrix form to perform a DFT. Unfortunately, this
%method is very computationally intensive, being O(n^2) rather than
%O(nlogn) for FFT. This method is not practical.

W = exp(1j*2*pi*[0:n-1]'*[0:n-1]/n); %This defines the DFT matrix

%{
cvx_begin
    variable Eap(n,n) complex
    maximize max(abs(W*Eap*W.'))
    subject to
        %Eap == 1/(4*pi)*ifft2(g*fft2(Enf))
        -20*log(abs(W*Eap*W.')) <= 0
        abs(W*Eap*W.') <= rho
cvx_end
surf(Enf)
figure
surf(Eap)
%}

%% Analysis and Conclusion
%  No headway was made with any of the attempted methods at solving the
%  optimization problem. I was not able to develop constraints on the
%  aperture that yielded anything but the trivial solution, and as such was
%  not able to develop a method to optimize the near-field shaping that
%  yielded anything but the trivial solution as well. I thought that an
%  initial guess might be fruitful, similar to the author's original paper,
%  but according to Michael Grant, creator of CVX, inputting an initial
%  guess is not supported by the software
%  (http://ask.cvxr.com/t/is-there-any-way-to-specify-the-starting-point-in-cvx/151).
%  Thus, my problem implementing the method described in the paper was in
%  failing to synthesize a realistic/accurate/similar electric field shape 
%  across the aperture (and thus in the near-field). Without a good method
%  to create this, implementing the propagation to the near-field was
%  unable to be properly tested or worked on.
