function M = calc_TransportMatrix(zob,zim, PS_Q0D, PS_Q1D, PS_Q2D, E, BEAMLINE)

% Define magnet positions
zQ0D = 1996.98244;
zQ1D = 1999.20656;
zQ2D = 2001.43099;

if ~exist('BEAMLINE')
    
    % Load the beamline file
    load 'specBeamline.mat'
end

   
% Find where everything is
Q0D =findcells(BEAMLINE,'Name','Q0D');
Q1D =findcells(BEAMLINE,'Name','Q1D');
Q2D =findcells(BEAMLINE,'Name','Q2D');

% Adjust the positions of things, and set magnet strengths
BEAMLINE{1} = BEAMLINE{Q0D(1)-2}; BEAMLINE{1}.S = zob;

BEAMLINE{2} = BEAMLINE{Q0D(1)-1};  BEAMLINE{2}.L = zQ0D-zob - 0.5;
BEAMLINE{3} = BEAMLINE{Q0D(1)};    BEAMLINE{3}.B = PS_Q0D/10/2;
BEAMLINE{4} = BEAMLINE{Q0D(2)};    BEAMLINE{4}.B = PS_Q0D/10/2;

BEAMLINE{5} = BEAMLINE{Q1D(1)-1};  BEAMLINE{5}.L = zQ1D- zQ0D - 1;
BEAMLINE{6} = BEAMLINE{Q1D(1)};    BEAMLINE{6}.B = PS_Q1D/10/2;
BEAMLINE{7} = BEAMLINE{Q1D(2)};    BEAMLINE{7}.B = PS_Q1D/10/2;

BEAMLINE{8} = BEAMLINE{Q2D(1)-1};  BEAMLINE{8}.L = zQ2D - zQ1D - 1;
BEAMLINE{9} = BEAMLINE{Q2D(1)};    BEAMLINE{9}.B = PS_Q2D/10/2;
BEAMLINE{10} = BEAMLINE{Q2D(2)};   BEAMLINE{10}.B = PS_Q2D/10/2;

BEAMLINE{11} = BEAMLINE{Q2D(2)+1}; BEAMLINE{11}.L = zim - zQ2D - 0.5;
BEAMLINE{12} = BEAMLINE{Q2D(2)+2}; 

% Fix all S positions
SetSPositions( 1, length(BEAMLINE),  BEAMLINE{1}.S);

% Set the energy of each element
for i =1:length(BEAMLINE)
    BEAMLINE{i}.P = E;
end

% Calculate the transport matrix
[~,M]=RmatAtoB(1, length(BEAMLINE));


end
