You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: CELES_MAIN.m
    path of file: ./repos/celes
    
    the code of the file until where you have to start completion: % =========================================================================
%> @brief Run this script to start the simulation
% =========================================================================

%% add all folders to the MATLAB search path
addpath(genpath('./src'))

%% import example file with sphere positions, radii and refractive indices
data = dlmread('examples/sphere_parameters.txt');

%% initialize the CELES class instances

% initialize particle class instance
%   - positionArray:        Nx3 array (float) in [x,y,z] format
%   - refractiveIndexArray: Nx1 array (complex) of complex refractive indices
%   - radiusArray:          Nx1 array (float) of sphere radii
particles = celes_particles('positionArray',        data(:,2:4), ...
                            'refractiveIndexArray', data(:,5)+1i*data(:,6), ...
                            'radiusArray',          data(:,1) ...
                            );

% initialize initial field class instance
%   - polarAngle:           scalar (float) polar angle of incoming beam/wave,
%                           in radians. for Gaussian beams, only 0 or pi are
%                           currently possible
%   - azimuthalAngle:       scalar (float) azimuthal angle of incoming
%                           beam/wave, in radians
%   - polarization:         string (char) polarization of incoming beam/wave
%                           ('TE' or 'TM')
%   - beamWidth:            scalar (float) width of beam waist. use 0 or inf
%                      
    
function (same unit as wavelength)
%   - solver:               valid instance of celes_solver class
numerics = celes_numerics('lmax',                   3, ...
                          'polarAnglesArray',       0:pi/5e3:pi, ...
                          'azimuthalAnglesArray',   0:pi/1e2:2*pi, ...
                          'gpuFlag',                true, ...
                          'particleDistanceResolution', 1, ...
                          'solver',                 solver);

% define a grid of points where the field will be evaulated
[x,z] = meshgrid(-4000:40:4000, -3000:40:5000); y = zeros(size(x));
% initialize output class instance
%   - fieldPoints:          Nx3 array (float) points where to evaluate the
%                           electric near field
%   - fieldPointsArrayDims: 1x2 array (int) dimensions of the array, needed to
%                           recompose the computed field as a n-by-m image
output = celes_output('fieldPoints',                [x(:),y(:),z(:)], ...
                      'fieldPointsArrayDims',       size(x));

% initialize simulation class instance
%   - input:                valid instance of celes_input class
%   - numerics:             valid instance of celes_input class
%   - output:               valid instance of celes_output class
simul = celes_simulation('input',                   input, ...
                         'numerics',                numerics, ...
                         'output',                  output);

%% run simulation
simul.run;

% evaluate transmitted and reflected power
simul.evaluatePower;
fprintf('transmitted power: \t%.4f %%\n', ...
    simul.output.totalFieldForwardPower/simul.output.initialFieldPower*100)
fprintf('reflected power: \t%.4f %%\n', ...
    simul.output.totalFieldBackwardPower/simul.output.initialFieldPower*100)

% evaluate field at output.fieldPoints
simul.evaluateFields;

%% plot results
% display particles
figure('Name','Particle positions','NumberTitle','off');
plot_spheres(gca,simul.input.particles.positionArray, ...
                 simul.input.particles.radiusArray, ...
                 simul.input.particles.refractiveIndexArray)

% plot near field
comp = ["real Ex", "real Ey", "real Ez", "abs E", "real Hx", "real Hy", "real Hz", "abs H"];
fieldtype = 'total field';

figure('Name','Near-field cross-cut','NumberTitle','off');
for sp = 1:numel(comp)
    subplot(2,4,sp)
    plot_field(gca, simul, comp(sp), fieldtype);
    caxis([-2*contains(comp(sp),'real'), 2])
end