! Simple tracking test of tracking with a matrix.
! tstTrack_matrix.madx

option,-echo,-info,warn;

! **********  Elements (FODO cell for simple tracking test) 
corh:  hkicker,    kick= 0.000;
qf:    multipole,  knl={0,0.01001};
qd:    multipole,  knl={0,-0.01};
dip:   multipole,  knl={twopi/50.0}, lrad= 1.1e3;
halfd: drift,      l= 5.0;
mk:    marker;
endmk: marker;


! ********** Define sequences
smach: sequence, l= 40.0;
corh,  at= 0;
qf,    at= 0;
dip,   at= 5.0;
qd,    at= 10.0;
dip,   at= 15.0;
mk,    at= 20.0;
qf,    at= 20.0;
dip,   at= 25.0;
qd,    at= 30.0;
dip,   at= 35.0;
mk,    at= 40.0;
endsequence;


! ********** Define beam and use the sequence with it 
beam,  sequence= smach, particle= proton, energy= 450;
use,   sequence= smach;
!show,  beam%smach;


! ********** TWISS
SAVEBETA, label= BETA1, place= mk[1]; ! using the current sequence
twiss; ! using the current sequence

! ********** Plot optics
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= BETX,BETY, VAXIS2= DX,DY;
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= ALFX,ALFY, VAXIS2= DPX,DPY;
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= MUX,MUY;
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= X,Y,VAXIS2= PX,PY;


! ********** Define matrix
BETAXX1= BETA1->BETX; ALFAXX1= BETA1->ALFX; DXX1= BETA1->DX;  DXXP1= BETA1->DPX;
BETAYY1= BETA1->BETY; ALFAYY1= BETA1->ALFY; MUXX= BETA1->MUX; MUYY=  BETA1->MUY;
MUXX1:= MUXX*TWOPI/360.;
MUYY1:= MUYY*TWOPI/360.;
R111:= cos(MUXX1)+ALFAXX1*sin(MUXX1);
R112:= BETAXX1*sin(MUXX1);
R116:= (1.-R111)*DXX1-R112*DXXP1;
R121:= -sin(MUXX1)/BETAXX1*(1+ALFAXX1^2);
R122:= cos(MUXX1)-ALFAXX1*sin(MUXX1);
R126:= (1.-R122)*DXXP1-R121*DXX1;
R151:= R121*R116-R111*R126;
R152:= R122*R116-R112*R126;
R133:= cos(MUYY1)+ALFAYY1*sin(MUYY1);
R134:= BETAYY1*sin(MUYY1);
R143:= -sin(MUYY1)/BETAYY1*(1+ALFAYY1^2);
R144:= cos(MUYY1)-ALFAYY1*sin(MUYY1);

value, DXX1,DXXP1,R111,R112,R116,R121,R122,R126;

MATSUP1:  MATRIX,RM11=R111,RM12=R112,RM21=R121,RM22=R122,
                 RM16=R116,RM26=R126,RM51=R151,RM52=R152,
                 RM33=R133,RM34=R134,RM43=R143,RM44=R144,
                 RM55=1.,RM66=1.,l=10.0;

! ********** Define new sequence including the matrix
smach2: sequence,l= 50.0;
corh,   at= 0;
qf,     at= 0;
dip,    at= 5.0;
qd,     at= 10.0;
dip,    at= 15.0;
mk,     at= 20.0;
qf,     at= 20.0;
dip,    at= 25.0;
qd,     at= 30.0;
dip,    at= 35.0;
mk,     at= 40.0;
MATSUP1,at= 45.0;
endmk,  at= 50.0;
endsequence;


! ********** Define beam and use the sequence with it 
beam,  sequence= smach2, particle= proton, energy= 450;
use,   sequence= smach2;
!show, beam%smach2;

! ********** TWISS
! This will give twice the tune as before
twiss; ! using the current sequence

! ********** Plot optics
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= BETX,BETY, VAXIS2= DX,DY;
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= ALFX,ALFY, VAXIS2= DPX,DPY;
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= MUX,MUY;
!PLOT, TABLE= TWISS, noversion, COLOUR= 100, HAXIS= S, VAXIS1= X,Y,VAXIS2= PX,PY;


! ********** Create user defined table
create, table= mytab, column= method, particle, x, px, y, py, t, pt;

! ********** Track 2 particles
track, onepass, dump; ! onepass = no closed orbit computed before tracking!
  start, x= 0.001,px= 4.492437849e-06;
  start, x= 0.1,  px= 4.492437849e-04;
  observe, place= endmk;
  run,   turns= 10;
endtrack;

method = "TRACK";
particle = 1;
x  = table(track.obs0002.p0001,x,1);
px = table(track.obs0002.p0001,px,1);
y  = table(track.obs0002.p0001,y,1);
py = table(track.obs0002.p0001,py,1);
t  = table(track.obs0002.p0001,t,1);
pt = table(track.obs0002.p0001,pt,1);
fill, table= mytab;

method = "TRACK";
particle = 2;
x  = table(track.obs0002.p0002,x,1);
px = table(track.obs0002.p0002,px,1);
y  = table(track.obs0002.p0002,y,1);
py = table(track.obs0002.p0002,py,1);
t  = table(track.obs0002.p0002,t,1);
pt = table(track.obs0002.p0002,pt,1);
fill, table= mytab;

!write, table= mytab, file= "out_tstTrack_matrix.tfs";
!system, "cat out_tstTrack_matrix.tfs";
!system, "rm track.obs*";


! ********** Check trajectory computation with twiss, Particle 1
select, flag= twiss, clear;
!select, flag= twiss, class= MATSUP1;
select, flag= twiss, class= endmk, column= name,s,x,px,y,py,t,pt,l;
TWISS,  betx= 110, bety= 110,
  x= 0.001, px= 4.492437849e-06, file= "twiss1.tfs";

method = "TWISS";
particle = 1;
x  = table(twiss,endmk,x);
px = table(twiss,endmk,px);
y  = table(twiss,endmk,y);
py = table(twiss,endmk,py);
t  = table(twiss,endmk,t);
pt = table(twiss,endmk,pt);
fill, table= mytab;
  
  
! ********** Check trajectory computation with twiss, Particle 2
select, flag= twiss, clear;
!select, flag= twiss, class= MATSUP1;
select, flag= twiss, class= endmk, column= name,s,x,px,y,py,t,pt,l;
TWISS,  betx= 110, bety= 110,
  x= 0.1, px= 4.492437849e-04, file= "twiss2.tfs";

method = "TWISS";
particle = 2;
x  = table(twiss,endmk,x);
px = table(twiss,endmk,px);
y  = table(twiss,endmk,y);
py = table(twiss,endmk,py);
t  = table(twiss,endmk,t);
pt = table(twiss,endmk,pt);
fill, table= mytab;


!*******************************************************************
!***  test traj computation with twiss particle(TRACK) ******

write, table= mytab, file= "mytab.tfs";

return;