%% helicalDoseCalcSetup.m

% RTF 11/4/05
% Sets up the run files needed for running a convolution
% superposition dose calculation with helical tomotherapy photon beams.

% Versions of the convolution code that run in both Windows XP, Condor, and
% directly through Matlab can all be created here.  The only real
% differences between the files are in the access conventions, specified by
% forward slashes or backslashes.

% Stephen R Bowen April 2009
% 
% The Kernels are generated from a cylindrical geometry within MCNP.
% The patient CT data is contained within the MATLAB structure Geometry
% and has the following format:
% 
% Geometry.voxel_size  (dx dy dz)
%         .start       (x y z)
%         .data        (M x N x Q single float matrix)
%         .ROIS        {1 x Nroi cell array}
%
% Each index of the Geometry.ROIS cell array should contain at least the
% following fields, even if manually created from Amira segmentation:
%
% Geometry.ROIS{1}.name         'string'
%                 .ind          (int32 sparse data array)
%                 .dim          (M N Q)
%                 .pix_dim      (dx dy dz)
%                 .start        (x y z)
%                 .start_ind    'string' 
%
% If your Geometry structure does not containg ROI information, then all
% beamlets will be included in the calculation.  Specified PTV ROIs will
% turn beamlets off that do not deposit dose within a specific margin.

%%

%%%%% Start of user supplied inputs %%%%%

% specify where kernel and geometry files are located
kernel_file = 'Kernels.mat';

geometry_file = '../plan/HN003/HN003final.mat';

% flags used to select which calculations will be done
Matlab_flag = 0;
Condor_flag = 1;
WinXP_flag = 0;

% define the overall beam field size (cm) for each beam angle in the
% tomotherapy coordinate system:
%
%   xp: MLC leaf axis
%   yp: MLC jaw axis
%
% remember yp points in the z-axis in the CT coordinate system
xpmin = -20;
xpmax = 20;
%ypmin = -0.05;  % total jaw width is 0.1 cm
%ypmax = 0.05;
%ypmin = -0.125; % total jaw width is 0.25 cm
%ypmax = 0.125;
ypmin = -0.3125;  % total jaw width is 0.625 cm
ypmax = 0.3125;
% ypmin = -0.5;    % total jaw width is 1 cm
% ypmax = 0.5;
% ypmin = -1.2250;  % total jaw width is 2.45 cm
% ypmax = 1.2250;


pitch_factor = 0.86;

% helical parameter: the number of cm the beam advances in the z-direction per gantry rotation
pitch = pitch_factor*(ypmax - ypmin);    % pitch factor times jaw width in z-direction
                                         % 0.86 determined to be optimal
                                         % (ref.Kissick)

% If a single PTV sinogram based beamlet optimization is used, specify PTV 
% ROI index number from Geometry file.  If no particular PTV is defined, 
% all PTV sinograms will be used to turn beamlets off.
indPTV = 18;

%%%%% End of user supplied inputs %%%%%

%%

load(geometry_file);

if size(Geometry.data,3) == 1
    Nrot = 1;
else
    Nrot = ceil(Geometry.voxel_size(3).*double(size(Geometry.data,3))./pitch);
end
                                  
% define the limits of the angles that will be used for the calculation
phimin = 0;  % starting angle in radians
phimax = 2*Nrot*pi;

% the number of beamlets that will fill in the above-defined grid in each direction
Mxp = 64/40*(xpmax - xpmin);
Nyp = 1;
Nphi = Nrot*51;  % number of angles used in the calculation

SAD = 85;  % source-axis distance for the x-ray source

beamlet_batch_size = 64;  % number of beamlets for condor to calculate in each batch

condor_folder = 'Condor';
winxp_folder = 'WinXP';
matlab_folder = 'Matlab';

code_folder = 'convolution';   % folder where convolution code resides, relative to this code

scripts_folder = 'scripts_m';  % dose calculation helper scripts

orig_code_folder = 'code';  % location of the convolution code that will be used

% create names for condor input and output folders
input_folder = 'input';
output_folder = 'output';
code_folder = 'source_code';

% name of the convolution/superposition executable, which will be in the
% 'code' folder of each respective run type folder
condor_exectuable_name = 'convolutionCondor';  % relative path on the cluster where code will be
winxp_executable_name = 'convolution.exe';
matlab_executable_name = 'convolution_mex';  % name of the Matlab version of the dose calculation code

% set the beam parameters, assuming a helical beam trajectory
% folders that will be inside the 'input' folder
beamspec_folder = 'beamspecfiles';    % directory where beam files will be stored
beamspec_batches_folder = 'beamspecbatches';
beamspec_batch_base_name = 'beamspecbatch';  % base name for a beamlet batch file
geometry_folder = 'geometryfiles';
kernel_folder = 'kernelfiles';  % folder where kernel information will be saved
kernel_filenames_condor = 'kernelFilenamesCondor.txt';
geometry_filenames_condor = 'geometryFilenamesCondor.txt';
kernel_filenames_winxp = 'kernelFilenamesWinXP.txt';
geometry_filenames_winxp = 'geometryFilenamesWinXP.txt';

% output folders
beamlet_batches_folder = 'beamletbatches';  % folder where resulting beamlet batches will be stored
beamlet_batch_base_name = 'beamletbatch';  % base name for a dose batch file
batch_output_folder = 'batchoutput';  % folder to which stdout will be printed

geometry_header_filename = 'geometryHeader.txt';
geometry_density_filename = 'density.bin';  % save the density, not the Hounsfield units!

% name for the condor submit file that will be used
condor_submit_file = 'convolutionSubmit.txt';
winxp_submit_file = 'convolutionSubmit.bat';  % WinXP run script format not known yet
matlab_submit_file = 'matlabSubmit.m';  % Matlab run script


% check the validity of the user-defined variables
if xpmin >= xpmax
    error('xpmin must be less than xpmax.');
end

if ypmin >= ypmax
    error('ypmin must be less than ypmax.');
end

if phimin > phimax
    error('phimin must be less than or equal to phimax.');
end

if Mxp <= 0 | Nyp <= 0 | Nphi <= 0
    error('Nxp, Nyp, and Nphi must be greater than zero.');
end

if SAD < 50
    error('It is recommended that the SAD be greater than 50 cm.');
end

% add the Matlab script and dos calculation folders to the path
path(path,[pwd '/' scripts_folder]);

% the xy plane is perpendicular to the isocenter axis of the linac gantry

% size of each beam aperture, making them vectors so extension to
% non-uniform aperture sizes becomes obvious
del_xp = (xpmax - xpmin)/Mxp;
del_yp = (ypmax - ypmin)/Nyp;

% Calculate the xp and yp offsets, which lie at the centers of the
% apertures.
xp = [xpmin:del_xp:xpmax-del_xp] + del_xp/2;
yp = [ypmin:del_yp:ypmax-del_yp] + del_yp/2;

% Gantry angles used.  These start at phimin and go to phimax, which is
% different from the way xp and yp are calculatd.
if Nphi == 1  % use the average of phimin and phimax of only one gantry angle desired
    phi = (phimin + phimax)/2;
else
    dphi = (phimax - phimin)/Nphi;
    phi = phimin:dphi:phimax;
end

% load data required for the dose calculator
load(kernel_file);

x_iso = Geometry.start(1) + (Geometry.voxel_size(1).*double(size(Geometry.data,1)))./2;
y_iso = Geometry.start(2) + (Geometry.voxel_size(2).*double(size(Geometry.data,2)))./2;
z_iso = Geometry.start(3);

iso = [x_iso y_iso z_iso];

if isfield(Geometry.ROIS) == 0
    sinogram = ones(Mxp,Nphi);   % all beamlets are turned on
elseif isfield(indPTV) == 1
% PTV sinogram specifies which beams are turned on
    RtomoAll = tomoBeamOnOffSinogram(Geometry.ROIS{indPTV},xpmin,xpmax,Mxp,Nphi,Nrot);
    sinogram = RtomoAll;
else
    sinogram = zeros(Mxp,Nphi);
    for i=1:length(Geometry.ROIS)
       if strncmp(Geometry.ROIS{i}.name,'PTV',3) == 1
           sinogram = sinogram + tomoBeamOnOffSinogram(Geometry.ROIS{i},xpmin,xpmax,Mxp,Nphi,Nrot);
       end
    end
    sinogram = double(~~sinogram);  
end


% ensure that all data are in single precision
f = fieldnames(Geometry);
for k=1:length(f)
    if isnumeric(getfield(Geometry,f{k}))
        Geometry = setfield(Geometry,f{k},single(getfield(Geometry,f{k})));    
    end
end

f = fieldnames(Kernels);
for k=1:length(f)
    if isnumeric(getfield(Kernels,f{k}))
        Kernels = setfield(Kernels,f{k},single(getfield(Kernels,f{k})));    
    end
end

% size of the CT data grid
[M,N,Q] = size(Geometry.data);

% Shift the isocenter
Geometry.start(1) = Geometry.start(1) - iso(1);
Geometry.start(2) = Geometry.start(2) - iso(2);

% Make double precision copies of Geometry attributes so Matlab
% mathematical operations can be performed:

% voxel size
dx = double(Geometry.voxel_size(1));
dy = double(Geometry.voxel_size(2));
dz = double(Geometry.voxel_size(3));

% start location
x0 = double(Geometry.start(1));
y0 = double(Geometry.start(2));
z0 = double(Geometry.start(3));

% calculate the coordinate system for the CT in case it's needed later
x = x0:dx:x0+(M-1)*dx;
y = y0:dy:y0+(N-1)*dy;
z = z0:dz:z0+(Q-1)*dz;

% find the total number of beams
Nbeam = Nphi*Mxp*Nyp;

batch_num = 0;  % start the count for the number of total batches

% check the geometry file to ensure that it's not in Hounsfield units
if length(find(Geometry.data > 20)) | length(find(Geometry.data < 0))
    error('Double check the Geometry structure, it may still be in Hounsfield units!');
end

% fill up a cell array of beam structures, grouped by batch
clear batches;
batch_num = 0;

batches = cell(1);  % start the batches cell array (cell array of beam batches)
batches{1} = cell(1);  % cell array of beams

overall_beam_num = -1;      % the first overall_beam_num will be zero

% calculate beams for all source directions and apertures
for k=1:Nphi  % loop through all gantry angles
    % calculate the source location for a helical trajectory
    beam.SAD = single(SAD);
    beam.y_vec = single([SAD*sin(phi(k)) -SAD*cos(phi(k)) pitch/(2*pi)*abs(phi(k)) + iso(3)]);
    % the kp vector ks the beam direction, ip and jp span the beam's eye view
    beam.ip = single([-cos(phi(k)) -sin(phi(k)) 0]);
    beam.jp = single([0 0 1]);
    beam.kp = single([-sin(phi(k)) cos(phi(k)) 0]);
    
    for n=1:Nyp  % loop through all apertures in yp-direction
        for m=1:Mxp  % loop through all apertures in the xp-direction
            % calculate the beam if the tomotherapy fluence value is non-zero
            if sinogram(m,k) > 0
                num = m + (n-1)*Mxp + (k-1)*Mxp*Nyp - 1;         % beamlet number (overall)
                overall_beam_num = overall_beam_num + 1;            
                batch_num = floor(overall_beam_num/beamlet_batch_size);       % batch number
                beam_num = mod(overall_beam_num,beamlet_batch_size);          % beam number in this batch
                
                % set the beam aperture parameters
                beam.del_xp = single(del_xp);
                beam.del_yp = single(del_yp);
                beam.xp = single(-xp(m));   % this parameter is flipped in order to get the tomotherapy delivery geometry correct
                beam.yp = single(yp(n));
                beam.num = single(num);  % record the beam number to avoid any later ambiguity
                
                batches{batch_num+1}{beam_num+1} = beam;
            end
        end
    end
end

% Everything else in this file is related to saving the batches in a
% useable form.
if Matlab_flag == 1
    % delete the old submission file
    err = rmdir(matlab_folder,'s');
    
    % create new folders for the submission
    mkdir([matlab_folder '/' input_folder '/' beamspec_batches_folder]);
    mkdir([matlab_folder '/' input_folder '/' kernel_folder]);
    mkdir([matlab_folder '/' input_folder '/' geometry_folder]);
    mkdir([matlab_folder '/' output_folder '/' batch_output_folder]);
    mkdir([matlab_folder '/' output_folder '/' beamlet_batches_folder]);
    
    % copy the convolution/superposition code to the submission folder
    copyfile(orig_code_folder,[matlab_folder '/' code_folder]);
    
    % save the kernels and geometry structures to the submission folder
    save([matlab_folder '/' input_folder '/' kernel_folder '/Kernels.mat'],'Kernels');
    fprintf(['Successfully saved Matlab kernels to ' input_folder '/' kernel_folder '\n']);
    
    save([matlab_folder '/' input_folder '/' geometry_folder '/Geometry.mat'],'Geometry');
    fprintf(['Successfully saved Matlab geometry to ' input_folder '/' geometry_folder '\n']);
end

if WinXP_flag == 1
    % delete the old submission file
    err = rmdir(winxp_folder,'s');
    
    % create folders where batch information will be sent
    mkdir([winxp_folder '/' input_folder '/' beamspec_batches_folder]);
    mkdir([winxp_folder '/' output_folder '/' beamlet_batches_folder]);
    mkdir([winxp_folder '/' output_folder '/' batch_output_folder]);
    
    % copy the convolution/superposition code to the submission folder
    copyfile(orig_code_folder,[winxp_folder '/' code_folder]);
    % copy the winxp executable to the submission folder
    copyfile([winxp_folder '/' code_folder '/Debug/' winxp_executable_name],[winxp_folder '/' winxp_executable_name]);
    
    % save the kernels and geometry files
    save_kernels(Kernels,[winxp_folder '/' input_folder '/' kernel_folder]);
    fprintf(['Successfully saved WinXP kernels to ' input_folder '/' kernel_folder '\n']);

    save_geometry(Geometry,[winxp_folder '/' input_folder '/' geometry_folder],geometry_header_filename,geometry_density_filename);
    fprintf(['Successfully saved WinXP geometry to ' input_folder '/' geometry_folder '\n']);
end

if Condor_flag == 1
    % delete the old submission file
    err = rmdir(condor_folder,'s');
    
    % create folders where batch information will be sent
    mkdir([condor_folder '/' input_folder '/' beamspec_batches_folder]);
    mkdir([condor_folder '/' output_folder '/' beamlet_batches_folder]);
    mkdir([condor_folder '/' output_folder '/' batch_output_folder]);
    
    % copy the current code folder to the new code folder
    copyfile(orig_code_folder,[condor_folder '/' code_folder]);
    
    % save the kernels and geometry, as this only needs to be done once
    save_kernels(Kernels,[condor_folder '/' input_folder '/' kernel_folder]);
    fprintf(['Successfully saved Condor kernels to ' input_folder '/' kernel_folder '\n']);
    
    save_geometry(Geometry,[condor_folder '/' input_folder '/' geometry_folder],geometry_header_filename,geometry_density_filename);
    fprintf(['Successfully saved Condor geometry to ' input_folder '/' geometry_folder '\n']);
end

% write the batches to files
for n=1:batch_num+1
    batch = batches{n};  % current batch
    
    if Matlab_flag == 1
        save([matlab_folder '/' input_folder '/' beamspec_batches_folder '/' beamspec_batch_base_name num2str(n) '.mat'],'batch');
    end
    
    if WinXP_flag == 1
        save_beamspec_batch(batch,[winxp_folder '/' input_folder '/' beamspec_batches_folder],[beamspec_batch_base_name num2str(n-1) '.txt']);
    end
    
    if Condor_flag == 1
        save_beamspec_batch(batch,[condor_folder '/' input_folder '/' beamspec_batches_folder],[beamspec_batch_base_name num2str(n-1) '.txt']);
    end
end

if Condor_flag == 1
    % create geometry filenames files
    fid = fopen([condor_folder '/' input_folder '/' geometry_filenames_condor],'w');
    fprintf(fid,'geometry_header\n');
    fprintf(fid,['./' input_folder '/' geometry_folder '/' geometry_header_filename '\n']);

    fprintf(fid,'geometry_density\n');
    fprintf(fid,['./' input_folder '/' geometry_folder '/' geometry_density_filename '\n']);
    fclose(fid);

    % create kernel filenames files
    fid = fopen([condor_folder '/' input_folder '/' kernel_filenames_condor],'w');
    fprintf(fid,'kernel_header\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/kernel_header.txt\n']);
    fprintf(fid,'kernel_radii\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/radii.bin\n']);
    fprintf(fid,'kernel_angles\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/angles.bin\n']);
    fprintf(fid,'kernel_energies\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/energies.bin\n']);
    fprintf(fid,'kernel_primary\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/primary.bin\n']);
    fprintf(fid,'kernel_first_scatter\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/first_scatter.bin\n']);
    fprintf(fid,'kernel_second_scatter\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/second_scatter.bin\n']);
    fprintf(fid,'kernel_multiple_scatter\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/multiple_scatter.bin\n']);
    fprintf(fid,'kernel_brem_annih\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/brem_annih.bin\n']);
    fprintf(fid,'kernel_total\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/total.bin\n']);
    fprintf(fid,'kernel_fluence\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/fluence.bin\n']);
    fprintf(fid,'kernel_mu\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/mu.bin\n']);
    fprintf(fid,'kernel_mu_en\n');
    fprintf(fid,['./' input_folder '/' kernel_folder '/mu_en.bin\n']);
    fclose(fid);

    % write the condor submit file
    beamspec_batch_filename = ['./' input_folder '/' beamspec_batches_folder '/' beamspec_batch_base_name '$(Process).txt'];
    beamlet_batch_filename = ['./' output_folder '/' beamlet_batches_folder '/' beamlet_batch_base_name '$(Process).bin'];
    fid = fopen([condor_folder '/' condor_submit_file],'w');
    fprintf(fid,'###############################################################\n');
    fprintf(fid,'# Condor submission script for convolution/superposition code\n');
    fprintf(fid,'###############################################################\n\n');
    fprintf(fid,'copy_to_spool = False\n');
    fprintf(fid,['Executable   = ' code_folder '/' condor_exectuable_name '\n']);
    fprintf(fid,['arguments    = ' input_folder '/' kernel_filenames_condor ' ' input_folder '/' geometry_filenames_condor ' ' beamspec_batch_filename ' ' beamlet_batch_filename '\n']);
    fprintf(fid,['Output       = ./' output_folder '/' batch_output_folder '/batchout$(Process).txt\n']);
    fprintf(fid,['Log          = ./' output_folder '/log.txt\n']);
    fprintf(fid,['Queue  ' num2str(batch_num + 1)]);
    fclose(fid);
    
    fprintf('Tomotherapy beamlet dose calculation setup complete.\n');
end

if WinXP_flag == 1
    % write the WinXP filename files
    fid = fopen([winxp_folder '/' input_folder '/' geometry_filenames_winxp],'w');
    fprintf(fid,'geometry_header\n');
    fprintf(fid,[input_folder '\\' geometry_folder '\\' geometry_header_filename '\n']);

    fprintf(fid,'geometry_density\n');
    fprintf(fid,[input_folder '\\' geometry_folder '\\' geometry_density_filename '\n']);
    fclose(fid);

    fid = fopen([winxp_folder '/' input_folder '/' kernel_filenames_winxp],'w');
    fprintf(fid,'kernel_header\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\kernel_header.txt\n']);
    fprintf(fid,'kernel_radii\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\radii.bin\n']);
    fprintf(fid,'kernel_angles\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\angles.bin\n']);
    fprintf(fid,'kernel_energies\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\energies.bin\n']);
    fprintf(fid,'kernel_primary\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\primary.bin\n']);
    fprintf(fid,'kernel_first_scatter\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\first_scatter.bin\n']);
    fprintf(fid,'kernel_second_scatter\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\second_scatter.bin\n']);
    fprintf(fid,'kernel_multiple_scatter\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\multiple_scatter.bin\n']);
    fprintf(fid,'kernel_brem_annih\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\brem_annih.bin\n']);
    fprintf(fid,'kernel_total\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\total.bin\n']);
    fprintf(fid,'kernel_fluence\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\fluence.bin\n']);
    fprintf(fid,'kernel_mu\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\mu.bin\n']);
    fprintf(fid,'kernel_mu_en\n');
    fprintf(fid,[input_folder '\\' kernel_folder '\\mu_en.bin\n']);
    fclose(fid);

    % line that will be used for running the code for the zeroth WinXP batch
    beamspec_batch_filename = [input_folder '\' beamspec_batches_folder '\' beamspec_batch_base_name num2str(0) '.txt'];
    beamlet_batch_filename = [output_folder '\' beamlet_batches_folder '\' beamlet_batch_base_name num2str(0) '.bin'];
    runline = [winxp_executable_name ' ' pwd '\' winxp_folder '\' input_folder '\' kernel_filenames_winxp ' ' pwd '\' winxp_folder '\' input_folder '\' geometry_filenames_winxp ' ' beamspec_batch_filename ' ' beamlet_batch_filename];
end

if Matlab_flag == 1
    % create a script to run the Matlab batches
    fid_matlab_submit = fopen([matlab_folder '\' matlab_submit_file],'w');

    % put all of the code into a cell array format for later printing
    sq = '''';  % single quote

    A = {...
        ['%% code to run matlab files'];
        ['Nbatches = ' num2str(prod(size(batches))) ';'];
        ['beamspec_batches_folder = ' sq beamspec_batches_folder sq ';'];
        ['beamlet_batches_folder = ' sq beamlet_batches_folder sq ';'];
        ['beamspec_batch_base_name = ' sq  beamspec_batch_base_name sq '; %% base name for a beamspec batch file'];
        ['beamlet_batch_base_name = ' sq  beamlet_batch_base_name sq '; %% base name for a beamlet batch file'];
        ['batch_output_folder = ' sq batch_output_folder sq ';  %% errors sent here in future version'];
        ['input_folder = ' sq input_folder sq ';'];
        ['output_folder = ' sq output_folder sq ';'];
        ['kernel_folder = ' sq kernel_folder sq ';'];
        ['geometry_folder = ' sq geometry_folder sq ';'];
        ['code_folder = ' sq code_folder sq ';'];
        ['matlab_executable_name = ' sq matlab_executable_name sq ';'];
        ['%% build the Matlab code and copy the executable file to current directory'];
        ['curr_dir = pwd;  %% current directory'];
        ['cd(code_folder);'];
        ['build;   %% build command for convolution code'];
        ['copyfile([matlab_executable_name ' sq '.dll' sq '],curr_dir);'];
        ['cd(curr_dir);  %% go back to the original directory'];
        [' '];
        ['%% load Geometry and Kernels'];
        ['load([input_folder ' sq '/' sq ' geometry_folder ' sq '/Geometry.mat' sq ']);'];
        ['load([input_folder ' sq '/' sq ' kernel_folder ' sq '/Kernels.mat' sq ']);'];
        [' '];
        ['for batch_num=1:Nbatches'];
        ['    fid_batch_error = fopen([output_folder ' sq '/' sq ' batch_output_folder ' sq '/batcherr' sq ' num2str(batch_num) ' sq '.txt' sq ' ' '],' sq 'w' sq ');  %% open error file for this batch'];
        ['    fprintf([' sq 'Calculating beamlets for batch ' sq ' num2str(batch_num) ' sq ' of ' sq ' num2str(Nbatches) ' sq '.\\n' sq ']);'];
        ['    %% open the next batch file'];
        ['    load([input_folder ' sq '/' sq ' beamspec_batches_folder ' sq '/' sq ' beamspec_batch_base_name num2str(batch_num) ' sq '.mat' sq ']);'];
        ['    Geometry.beam = batch;'];
        ['    try'];
        ['        eval([' sq 'beamlets = ' sq ' matlab_executable_name ' sq '(Kernels,Geometry);' sq ']);'];
        ['        save([output_folder ' sq '/' sq ' beamlet_batches_folder ' sq '/' sq ' beamlet_batch_base_name num2str(batch_num) ' sq '.mat' sq '], ' sq 'beamlets' sq ');'];
        ['    catch'];
        ['        fprintf([' sq 'Error in the calculation for batch ' sq ' num2str(batch_num) ' sq '\\n' sq ']); %% print to stdout'];
        ['        fprintf(fid_batch_error,[' sq 'Error in the calculation for batch ' sq ' num2str(batch_num) ' sq '\\n' sq ']); %% print to error file'];
        ['    end'];
        ['end']};
    
    % print the cell array to the matlab submission file
    for k=1:length(A)
        fprintf(fid_matlab_submit,[A{k} '\n']);    
    end
    fclose(fid_matlab_submit);
end