rm(list=ls())  # clear all variable values

# set up matrices and define a complete cell:
Mf = matrix(c( 1,  0, 0, -1/25, 1, 0    , 0, 0, 1), byrow=TRUE, ncol=3)
Md = matrix(c( 1,  0, 0,  1/25, 1, 0    , 0, 0, 1), byrow=TRUE, ncol=3)
Mo = matrix(c( 1, 15, 0,   0 ,  1, 0    , 0, 0, 1), byrow=TRUE, ncol=3)
Mb = matrix(c( 1,  0, 0,   0 ,  1, 0.032, 0, 0, 1), byrow=TRUE, ncol=3)
#     this last element is a "bend" by an angle of 32 mrad

cell = list( Mf, Mo, Mb, Mo, Md, Mo, Mb, Mo)
ds   =    c(  0, 15,  0, 15,  0, 15,  0, 15)
s    =   cumsum(ds)

# find the matrix of a complete cell:
Mcell = diag(1,3)
i = 0
while(i<length(cell)){
   i=i+1
  Mcell =cell[[i]] %*% Mcell 
}
Mcell
Mcell2 = Mcell[1:2,1:2]

# compute the periodic functions at the end point of the cell:
mu   = acos(sum(diag(Mcell2))/2)
bet0 = Mcell2[1,2]/sin(mu)
alf0 = (Mcell2[1,1]-Mcell2[2,2])/2/sin(mu)
D0   = solve( diag(1,2)-Mcell2, Mcell[1:2,3] )
# Note:  the function solve()
#     solve(M,V)  solves M X = V   for X
#     solve(M)    finds the inverse of M

# results:
bet0; alf0; D0[1]

# Track the functions
# first, beta...
K = matrix( c(bet0, -alf0, -alf0, (1+alf0^2)/bet0), byrow=TRUE,ncol=2)
bet = K[1,1]
D = c(D0[1],D0[2],1)
Dsp  = D[1]
i=0
while(i<length(cell)){
   i=i+1
   K=cell[[i]][1:2,1:2] %*% K %*% t(cell[[i]][1:2,1:2])
   bet[i] = K[1,1]
   D = cell[[i]] %*% D
   Dsp[i] = D[1]
}

# plot the results:
plot(s,bet,typ="l",ylim=c(0,120),main="beta [black] and dispersion*10 [blue]")
points(s,Dsp*10,type="l",col="blue")