function [Eout,Gout] = Focus_beam_with_telescope(Ein,array_L_D)
% Focus_beam_with_lens() use this function for lens with short focal length
%  This function re-adapt the grid size since the focusing may lead to a
%  small beam size
%  array_L_D is a vector representing the focal length of the lens and the
%  distance of propagation, several lens and distance can be simulated. For
%  example if array_L_D = [100 5 -10 4] means that the beam first cross a
%  lens of focal length 100m then propagate 5m, then meet another lens of
%  focal length -10m and finally propagate 4m.

if ~isa(Ein, 'E_Field')
    error('Focus_beam_with_lens2(): the first argument must be an instance of the class E_Field')
end

if rem(length(array_L_D),2)
    error('Focus_beam_with_lens2(): the input array must have a odd  number of element')
end

% Nomber of iteration for the focusing and propagation:
num_iter = length(array_L_D) / 2;

% first derive the parameters of the input beam and remove the wavefront
[~,R_in] = Fit_TEM00(Ein);
WF_change = exp(1i * Ein.k_prop *  Ein.Grid.D2_square * (1/(2*R_in)) );
Ein.Field = Ein.Field .* WF_change;

for pp=1:num_iter
    % Change the sign of the lens to be consitent with the following
    f_lens = -array_L_D(2*(pp-1)+1);
    
    if (f_lens + R_in) ~= 0
        f_new = -(f_lens * R_in)/(f_lens + R_in);
    else
        f_new = 1E99; % So the lens will cancel the incident wavefront of curvature
    end
    
    % Corrected distance to propagate:
    dist_prop = - array_L_D(2*(pp-1)+2) * f_new / (array_L_D(2*(pp-1)+2) - f_new);
    
    E_prop = Propagate_E(Ein,dist_prop);
    
    Scaling_factor = (f_new - array_L_D(2*(pp-1)+2)) / f_new;
    
    Gout = Grid(Ein.Grid.Num_point,Ein.Grid.Length*abs(Scaling_factor));
    
    Eout = Ein;
    Eout.Grid = Gout;
    Eout.Field = E_prop.Field / Scaling_factor;
    
     % Check if no power fall outside the grid
    Check_Grid_size(E_prop,0.10)    
    New_wavefront = -(array_L_D(2*(pp-1)+2) - f_new);
    
    % Remove the wavefront curvature from the beam
    [~,R_in2] = Fit_TEM00(Eout);
    WF_change = exp(1i * Ein.k_prop *  Eout.Grid.D2_square * (1/(2*R_in2)) );
    Eout.Field = Eout.Field .* WF_change;
    
    % Calculate the new wavefront
    R_in = 1/(1/R_in2 - 1/New_wavefront);
    Ein = Eout;
    
end

% At the end, add the final wavefront
Eout.Field = Eout.Field .* exp(1i * Eout.k_prop *  Eout.Grid.D2_square * (-1/(2*R_in)) );

end

