QF = OptElem(0.2, Quadrupole(-1.2, 0.2));
QD = OptElem(0.4, Quadrupole(1.2, 0.4));
B  = OptElem(1.5, Dipole(3.8197, 1.5));
EB = OptElem(0.0, Edge(3.8197, 0.1964));
D  = OptElem(0.55, Drift(0.55));

accelerator=list(QF,D,EB,B,EB,D,QD,D,EB,B,EB,D,QF,QF,D,EB,B,EB,D,QD,D,EB,B,EB,D,QF,QF,D,EB,B,EB,D,QD,D,EB,B,EB,D,QF,QF,D,EB,B,EB,D,QD,D,EB,B,EB,D,QF);
m=OptElem(0.,eye(5,5));

x=[0.;0.0001;0.;0.;0.];
s=0.;
for k=accelerator
  m=k*m;
  s=[s,m.s];
  x=[x,m.matrix*x(:,1)];
end

xbasc()
plot(s,x(1,:))
xtitle("Beam evolution", "s [m]", "x [m]")
grid(2)
