function Cout = Calculate_RT_mat(Cin,varargin)
%Calculate_RT_mat() Calculate the kernel for one round trip in the cavity
% C1 = Calculate_RT_mat(C1), this function calculate the kernel for one
% round trip of the light in the cavity. From this kernel, one can derive
% the eigen modes and eigen vectors of the cavity.
% !! only use with small size of grid 64X64  with 4GB RAM, or 128X128 on
% more powerful machine

if ~isempty(Cin.Cavity_EM_mat)
    disp('Calculate_RT_mat(): Cavity kernel has already been calculated  ')
end

resampled_grid = false;

if nargin == 2
    if ~isa(varargin{1}, 'Grid')
        error(' Calculate_RT_mat(): the second argument must be a grid')
    else
        resampled_grid = true;
        Gr = varargin{1};
    end
end

Cout = Cin;

if ~resampled_grid
    Num_point = Cin.Laser_in.Grid.Num_point;
else
    Num_point = Gr.Num_point;
end

tmp_mat_EM =  complex(zeros(Num_point^2));

for pp=1:Cin.Nb_mirror
    Cin.Propagation_mat_array(pp).Use_DI = true;
end

disp('Calculating the kernel:       ')

if license('test','distrib_computing_toolbox')          % check if the Parallel Computing Toolbox exists
    
    matlabpool open 4
    
    for mx=1:Num_point
        
        Tmp_RT_mat2 = complex(zeros(Num_point,Num_point^2));
        
        parfor py=1:Num_point
            
            if ~resampled_grid
                E_in = Cin.Laser_in;
                E_in.Field = complex(zeros(Num_point));
                E_in.Field(mx,py) = 1;
                Field_Circ = E_in;
            else
                E_in = E_Field(Gr,'w0',450E-6); % dummy value
                E_in.Field = complex(zeros(Num_point));
                E_in.Field(mx,py) = 1;
                Field_Circ = Resample_E(E_in,Cin.Laser_in.Grid);
            end
            
            for pp=1:Cin.Nb_mirror
                if pp ~= Cin.Nb_mirror % check we are not at the last iteration
                    Field_Circ = Propagate_E(Field_Circ,Cin.Propagation_mat_array(pp));
                    Field_Circ = Reflect_mirror(Field_Circ,Cin.I_array(pp+1),'Ref',1);
                else
                    Field_Circ = Propagate_E(Field_Circ,Cin.Propagation_mat_array(pp));
                    Field_Circ = Reflect_mirror(Field_Circ,Cin.I_array(1),'Ref',1);
                end
            end
            
            if ~resampled_grid
                Tmp_RT_mat2(py,:)  =  Field_Circ.Field(:);
            else
                Field_tmp = Resample_E(Field_Circ,Gr);
                Tmp_RT_mat2(py,:) = Field_tmp.Field(:);
            end
        end
        
        tmp_mat_EM((mx-1)*Num_point+1:(mx-1)*Num_point+Num_point,:) = Tmp_RT_mat2;
    end
    
    matlabpool close
    
else % if the PCT is not installed
    
    for mx=1:Num_point
        for py=1:Num_point
            
            if ~resampled_grid
                E_in = Cin.Laser_in;
                E_in.Field = complex(zeros(Num_point));
                E_in.Field(mx,py) = 1;
                Field_Circ = E_in;
            else
                E_in = E_Field(Gr,'w0',450E-6); % dummy value
                E_in.Field = complex(zeros(Num_point));
                E_in.Field(mx,py) = 1;
                Field_Circ = Resample_E(E_in,Cin.Laser_in.Grid);
            end
            
            for pp=1:Cin.Nb_mirror
                if pp ~= Cin.Nb_mirror % check we are not at the last iteration
                    Field_Circ = Propagate_E(Field_Circ,Cin.Propagation_mat_array(pp));
                    Field_Circ = Reflect_mirror(Field_Circ,Cin.I_array(pp+1),'Ref',1);
                else
                    Field_Circ = Propagate_E(Field_Circ,Cin.Propagation_mat_array(pp));
                    Field_Circ = Reflect_mirror(Field_Circ,Cin.I_array(1),'Ref',1);
                end
            end
            
            if ~resampled_grid
                tmp_mat_EM((mx-1)*Num_point+py,:) = Field_Circ.Field(:);
            else
                Field_tmp = Resample_E(Field_Circ,Gr);
                tmp_mat_EM((mx-1)*Num_point+py,:) = Field_tmp.Field(:);
                
            end
        end
        fprintf('\b\b\b\b\b\b\b %-3.0i %% ',round(100 * mx / Num_point))
    end
 fprintf('\b\b\b\b\b\b\b done! \n')            
    
end

Cout.Cavity_EM_mat = tmp_mat_EM;
















