function [y]=Drift(s)   // Free drift
  y=[1, s;0, 1];
endfunction

function [y]=Quadrupole(k,s)    // quadrupole function
  if k<0 then // focussing
    Omega=sqrt(-k)*s
    y=[cos(Omega),sin(Omega)/sqrt(-k); -sqrt(-k)*sin(Omega) cos(Omega)];
  else  //defocussing
    Omega=sqrt(k)*s
    y=[cosh(Omega),sinh(Omega)/sqrt(k); sqrt(k)*sinh(Omega) cosh(Omega)];
  end
endfunction

// typed list definition of an optical element
function [y]=OptElem( s, mat )
  offset = [0.; 0.];
  y=tlist(['ele';'s';'matrix';'offset'],s, mat, offset)
endfunction

// overloaded multiplication of two optical elements
function [z]=%ele_m_ele(a,b)  
  z=OptElem( a.s+b.s, a.matrix*b.matrix )
endfunction

// overloaded transport of a vector through an optical element
function [z]=%ele_m_s(a,x)  
  z = a.matrix*(x-a.offset)+a.offset
endfunction

sigma_survey = 0.0003;  // position accuracy

//  define generic components
ds=20;qs=1.; k=0.05;
O=OptElem( ds, Drift(ds) );
F=OptElem( qs, Quadrupole(-k, qs ) );
D=OptElem( qs, Quadrupole(k, qs) );
d=OptElem( 0., [1 0; 0 1] );

F.offset = [rand( sigma_survey, 'normal' ) * sigma_survey; 0.];
D.offset = [rand( sigma_survey, 'normal' ) * sigma_survey; 0.];
FODO = list( F, O, D, O );
accelerator = FODO;
N = 200 - 1;
for k=1:N
  F.offset = [rand( sigma_survey, 'normal' ) * sigma_survey; 0.];
  D.offset = [rand( sigma_survey, 'normal' ) * sigma_survey; 0.];
  FODO = list( F, O, D, O );
  accelerator= lstcat( accelerator, FODO );
end

x=[ 0.; 0. ];
s=0.;
sigma = x;
for k=accelerator
  d = k*d;
  s = [s,d.s];
  x = k*x;
  sigma = [sigma, x ];
end

xbasc()
xtitle( 'Exercise 3', 's [m]', 'sigma' );
plot( s, sigma );

