function Cout = Calculate_RT_mat(Cin)
%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
% a more powerful machine

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

Cout = Cin;

Num_point = Cin.Laser_in.Grid.Num_point;
tmp_mat_EM =  complex(zeros(Num_point^2));

Cin.Propagation_mat.Use_DI = true;

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
            
            E_in = Cin.Laser_in;
            E_in.Field = complex(zeros(Num_point));
            E_in.Field(mx,py) = 1;
            
            Circ_field = Propagate_E(E_in,Cin.Propagation_mat);
            Circ_field = Reflect_mirror(Circ_field,Cin.I_end,'Ref',1);
            Circ_field = Propagate_E(Circ_field,Cin.Propagation_mat);
            Circ_field = Reflect_mirror(Circ_field,Cin.I_input,'Ref',1);
            
            Tmp_RT_mat2(py,:)  =  Circ_field.Field(:);
            
        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
            
            E_in = Cin.Laser_in;
            E_in.Field = complex(zeros(Num_point));
            E_in.Field(mx,py) = 1;
            
            Circ_field = Propagate_E(E_in,Cin.Propagation_mat);
            Circ_field = Reflect_mirror(Circ_field,Cin.I_end,'Ref',1);
            Circ_field = Propagate_E(Circ_field,Cin.Propagation_mat);
            Circ_field = Reflect_mirror(Circ_field,Cin.I_input,'Ref',1);
            
            tmp_mat_EM((mx-1)*Num_point+py,:) =  Circ_field.Field(:);            
        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;

end

