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));

C=QF*D*EB*B*EB*D*QD*D*EB*B*EB*D*QF;

Beta_a_x = sqrt(-C.matrix(1,2)*C.matrix(2,2)/(C.matrix(2,1)*C.matrix(1,1)));
Beta_b_x = -C.matrix(1,2)/(Beta_a_x*C.matrix(2,1));

Beta_a_z = sqrt(-C.matrix(3,4)*C.matrix(4,4)/(C.matrix(4,3)*C.matrix(3,3)));
Beta_b_z = -C.matrix(3,4)/(Beta_a_z*C.matrix(4,3));

Disp_a_x = -C.matrix(2,5)/C.matrix(2,1);
Disp_b_x = -C.matrix(1,1)*C.matrix(2,5)/C.matrix(2,1)+C.matrix(1,5);

B1= OptFunc(0., Beta_a_x, 0., Beta_a_z, 0., Disp_a_x, 0.);
B2=Beta(C, B1);

Accelerator = list(QF,D,EB,B,EB,D,QD,D,EB,B,EB,D,QF);

M = OptElem(0., eye(5,5));

s = 0;
beta_x = Beta_a_x;
beta_z = Beta_a_z;
d = Disp_a_x;
for K=Accelerator
  M = K*M;
  B2 = Beta(M, B1);
  s = [s, B2.s];
  beta_x = [beta_x, B2.Beta(1,1)];
  beta_z = [beta_z, B2.Beta(3,3)];
  d = [d, B2.Disp(1)];
end

xbasc();

plot2d(s', [beta_x' beta_z' d'], leg="beta_x@beta_z@disp", style=[5 2 3]);
xtitle("Accelerator functions", "s [m]", "beta,disp [m]");
xgrid();
