! this statements are needed to have explicit values in the output file
SYKICK=0;
S20SY02.MKDC010.cory:= SYKICK*2;
S20SY02.MKAC020.cory:= SYKICK*3;
S20SY02.MKDC030.cory:= SYKICK*2;
S20SY02.MKAC040.cory:= SYKICK*3;
S20SY02.MKDC050.cory:= SYKICK*2;


s20bc01.mqua010.k1:=1.;
Print,text='MATCH POINT: BC2 in Linac 2';
MATCH,SEQUENCE=Swissfel,range=s20sy02$start/s20bc02$end,beta0=twisssep;
VARY,NAME=s20bc01.mqua010.k1,STEP=0.0001;
VARY,NAME=s20bc01.mqua030.k1,STEP=0.0001;
VARY,NAME=s20bc01.mqua060.k1,STEP=0.0001;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20bc02$end,betx<4;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20bc02$end,bety<90;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20bc01$end,betx<250;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20bc01$end,bety<20;
LMDIF,CALLS=1000,TOLERANCE=1.e-21;
ENDMATCH;

twiss,SEQUENCE=swissfel,range=s20sy02$start/S20BC02$END,beta0=twisssep;
plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;



muLin3=0.2;
pol = 1;
s20cb01.mqua430.k1=0.5*pol;
s20cb02.mqua430.k1=-0.5*pol;

Print,text='MATCH POINT: FODO - Linac 2';
match,SEQUENCE=swissfel,range=S20cb01$start/s20cb02$end;
Vary,name=S20CB01.MQUA430.k1,step=0.0001;
Vary,name=S20CB02.MQUA430.k1,step=0.0001;
constraint,sequence=swissfel,range=s20cb02$end,mux=muLin3;
constraint,sequence=swissfel,range=s20cb02$end,muy=muLin3;
lmdif,calls=100,tolerance=1e-21;
endmatch;

twiss,SEQUENCE=swissfel,range=s20cb01$start/s20cb02$end;

nbetx=table(twiss,s20cb02$end,betx);
nalfx=table(twiss,s20cb02$end,alfx);
nbety=table(twiss,s20cb02$end,bety);
nalfy=table(twiss,s20cb02$end,alfy);
plot,haxis=s,vaxis=betx,bety,colour=100;



s20cb03.mqua430.k1 = s20cb01.mqua430.k1;
s20cb04.mqua430.k1 = s20cb02.mqua430.k1;
s30cb01.mqua430.k1 = s20cb01.mqua430.k1;
s30cb02.mqua430.k1 = s20cb02.mqua430.k1;
s30cb03.mqua430.k1 = s20cb01.mqua430.k1;
s30cb04.mqua430.k1 = s20cb02.mqua430.k1;
s30cb05.mqua430.k1 = s20cb01.mqua430.k1;
s30cb06.mqua430.k1 = s20cb02.mqua430.k1;
s30cb07.mqua430.k1 = s20cb01.mqua430.k1;
s30cb08.mqua430.k1 = s20cb02.mqua430.k1;
s30cb09.mqua430.k1 = s20cb01.mqua430.k1;
s30cb10.mqua430.k1 = s20cb02.mqua430.k1;
s30cb11.mqua430.k1 = s20cb01.mqua430.k1;
s30cb12.mqua430.k1 = s20cb02.mqua430.k1;
s30cb13.mqua430.k1 = s20cb01.mqua430.k1;
s30cb14.mqua430.k1 = s20cb02.mqua430.k1;




Print,text='MATCH POINT: Match into Linac 2';
MATCH,SEQUENCE=Swissfel,range=s20sy02$start/s20cb04$end,beta0=twisssep;
VARY,NAME=s20ma01.mqua010.k1,STEP=0.0001;
VARY,NAME=s20ma01.mqua020.k1,STEP=0.0001;
VARY,NAME=s20ma01.mqua040.k1,STEP=0.0001;
vary,name=s20cb01.mqua430.k1,step=0.0001;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20cb03$start,betx=nbetx;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20cb03$start,bety=nbety;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20cb03$start,alfx=nalfx;
CONSTRAINT,SEQUENCE=SwissFEL,range= s20cb03$start,alfy=nalfy;
LMDIF,CALLS=1000,TOLERANCE=1.e-21;
ENDMATCH;

twiss,SEQUENCE=swissfel,range=s20sy02$start/S20cb04$end,beta0=twisssep;
plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;



Print,text='MATCH POINT: Corrected Matching into Linac 3';
MATCH,SEQUENCE=Swissfel,range=s20sy02$start/s30cb07$end,beta0=twisssep;
VARY,NAME=s30cb01.mqua430.k1,STEP=0.0001;
VARY,NAME=s30cb02.mqua430.k1,STEP=0.0001;
VARY,NAME=s30cb03.mqua430.k1,STEP=0.0001;
VARY,NAME=s30cb04.mqua430.k1,STEP=0.0001;
CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb05$start,betx=nbetx;
CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb05$start,bety=nbety;
CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb05$start,alfx=nalfx;
CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb05$start,alfy=nalfy;
LMDIF,CALLS=1000,TOLERANCE=1.e-21;
ENDMATCH;




twiss,SEQUENCE=swissfel,range=s20sy02$start/S30cb14$END,beta0=twisssep;
plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;





PSYKICK:=0.1;
S30cb15.MKDC010.cory := PSYKICK*2;
S30cb15.MKAC020.cory := PSYKICK*3;
S30cb15.MKDC030.cory := PSYKICK*2;
S30cb15.MKAC040.cory := PSYKICK*3;
S30CB15.MKDC050.cory := PSYKICK*2;
pqf:=1;

S30cb15.mqua100.k1 := -0.2*pqf;
S30cb15.mqua200.k1 := -pqf;

Print,text='MATCH POINT: Kicker Strength for Porthos Septum';
MATCH,SEQUENCE=Swissfel,range=S30CB15$START/s30cb15$END,BETA0 = TWISSsep;
VARY, NAME=PSYKICK, STEP=0.00001;
VARY, NAME=PQF, STEP=0.00001;
CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end, Y = 0.01;
CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end, PY = 0.;
LMDIF, CALLS=1000, TOLERANCE=1.E-21;
ENDMATCH;


s30cb13.mqua430.k1 = 0.456;
s30cb14.mqua430.k1 = -1.741;
s30cb15.mqua010.k1 = 1.1533;
s30cb15.mqua020.k1 = -0.0961;



Print,text='MATCH POINT: Match into Porthos Septum';
MATCH,SEQUENCE=Swissfel,range=s20sy02$start/s30sy01$end,beta0=twisssep;
VARY,NAME=s30cb13.mqua430.k1,STEP=0.0001;
VARY,NAME=s30cb14.mqua430.k1,STEP=0.0001;
VARY,NAME=s30cb15.mqua010.k1,STEP=0.0001;
VARY,NAME=s30cb15.mqua020.k1,STEP=0.0001;
CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end,betx=twisspor->betx;
CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end,bety=twisspor->bety;
CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end,alfx=twisspor->alfx;
CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end,alfy=twisspor->alfy;
!LMDIF, CALLS=1000, TOLERANCE=1.E-21;
Simplex, CALLS=1000, TOLERANCE=1.E-21;
ENDMATCH;

twiss,SEQUENCE=swissfel,range=s20sy02$start/s30cb15$end,beta0=twisssep;
plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;
plot,haxis=s,vaxis=y,dy,range=s30cb15$start/#e,colour=100;



S30cb15.MKDC010.cory := 0;
S30cb15.MKAC020.cory := 0;
S30cb15.MKDC030.cory := 0;
S30cb15.MKAC040.cory := 0;
S30CB15.MKDC050.cory := 0;




!--------------------------------------------------
! Match ecol

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
! energy collimator
if (ECOLasBC ==0.) {
    sarcl02.mqua130.k1=1.8;
    R56 = -000e-6;

    Print,text='MATCH POINT: Isochronous Setup of Energy Collimator';
    MATCH,RMATRIX,SEQUENCE=swissfel,range=SARCL02$START/SARCL02$END,betx=10,alfx=0,bety=10,alfy=0;
    VARY,NAME=sarcl02.mqua130.k1,STEP=0.0001;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl02.mqua210,re56=R56*0.5;
    LMDIF,CALLS=100,TOLERANCE=1.e-21;
    ENDMATCH;

    sarcl02.mqua210.k1=-sarcl02.mqua130.k1*0.9;
    sarcl02.mqua250.k1= sarcl02.mqua130.k1;

    Print,text='MATCH POINT: Dispersion Symmetry Point of Energy Collimator';
    MATCH,SEQUENCE=swissfel,range=SARCL02$START/SARCL02$END,betx=10,alfx=0,bety=10,alfy=0;
    VARY,NAME=sarcl02.mqua250.k1,STEP=0.0001;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl02.dbpm260$end,dpx=0;
    LMDIF,CALLS=100,TOLERANCE=1.e-21;
    ENDMATCH;

    sarcl02.mqua310.k1=sarcl02.mqua250.k1;
    sarcl02.mqua350.k1=sarcl02.mqua210.k1;
    sarcl02.mqua460.k1=sarcl02.mqua130.k1;

    dL=0.625+0.25;
    bx0=25;
    bx=bx0+dL*dL/bx0;
    ax=-dL/bx0;
    by0=0.8;
    by=by0+dL*dL/by0;
    ay=-dL/by0;

    twiss,sequence=swissfel,range=sarcl02.mqsk300/sarcl02$end,betx=bx,alfx=ax,bety=by,alfy=ay;
    plot,haxis=s,vaxis=betx,bety,colour=100;

    nbetx= table(twiss,sarcl02$end,betx);
    nalfx=-table(twiss,sarcl02$end,alfx);
    nbety= table(twiss,sarcl02$end,bety);
    nalfy=-table(twiss,sarcl02$end,alfy);


    sarcl02.msex255.k2:=110;
    sarcl02.msex305.k2:=110;

    Print,text='MATCH POINT: Sextupoles in Energy Collimator';
    MATCH,chrom,SEQUENCE=sarcl02,range=#s/#e,betx=nbetx,alfx=nalfx,bety=nbety,alfy=nalfy;
    VARY,NAME=ecol_k2,STEP=0.0001;
    CONSTRAINT,SEQUENCE=Sarcl02,range=#e,ddx=0;
    LMDIF,CALLS=100,TOLERANCE=1.e-21;
    ENDMATCH;

    sarcl02.msex255.k2=ecol_k2;
    sarcl02.msex305.k2=ecol_k2;

    twiss,chrom,rmatrix,range=SARCL02$START/SARCL02$END,betx=nbetx,alfx=nalfx,bety=nbety,alfy=nalfy;
    plot,haxis=s,vaxis=betx,colour=100;
    plot,haxis=s,vaxis=bety,colour=100;
    plot,haxis=s,vaxis=dx,colour=100;
    plot,haxis=s,vaxis=ddx,colour=100;
    plot,haxis=s,vaxis=re56,colour=100;


} else {
    sarcl02.mqua130.k1=0;
    sarcl02.mqua210.k1=0;
    sarcl02.mqua250.k1=0;
    sarcl02.mqua310.k1=0;
    sarcl02.mqua350.k1=0;
    sarcl02.mqua460.k1=0;
    sarcl02.msex255.k2=0;
    sarcl02.msex305.k2=0;
}







!---------------------------------
! match into energy collimator
! reference point is sarma01-mqua010


sarcl01.mqua020.k1 := 0.4550820417;
sarcl01.mqua050.k1 := -0.02350776516;
sarcl01.mqua080.k1 := -0.3705426893;
sarcl01.mqua100.k1 := 0.3937186328;
sarcl01.mqua140.k1 := 1.214406656;
sarcl01.mqua190.k1 := -1.475278354;


use,sequence=swissfel;
Print,text='MATCH POINT: Match into Energy Collimator';
MATCH,SEQUENCE=swissfel,range=s20SY02$start/sarma02$start,beta0=twisssep;
VARY,NAME=sarcl01.mqua020.k1,STEP=0.0001,lower=-1., upper =1.;
VARY,NAME=sarcl01.mqua050.k1,STEP=0.0001;
VARY,NAME=sarcl01.mqua080.k1,STEP=0.0001;
VARY,NAME=sarcl01.mqua100.k1,STEP=0.0001;
VARY,NAME=sarcl01.mqua140.k1,STEP=0.0001;
VARY,NAME=sarcl01.mqua190.k1,STEP=0.0001;
if (ECOLasBC ==0.){
    CONSTRAINT,SEQUENCE=swissfel,range=sarma01.mqua010$start,betx=twissARECOL->BETX;
    CONSTRAINT,SEQUENCE=swissfel,range=sarma01.mqua010$start,alfx=twissARECOL->ALFX;
    CONSTRAINT,SEQUENCE=swissfel,range=sarma01.mqua010$start,bety=twissARECOL->BETY;
    CONSTRAINT,SEQUENCE=swissfel,range=sarma01.mqua010$start,alfy=twissARECOL->ALFY;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua050,bety<30;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua190,betx<30;
 } else {
     CONSTRAINT,SEQUENCE=swissfel,range=sarcl02.dscr280$start,betx  = 10;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl02.dscr280$start,alfx=0.14;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl02.dscr280$start,bety  = 10;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl02.dscr280$start,alfy=0.14;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua050,bety<70;
    CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua190,betx<70;
 }
LMDIF,CALLS=100,TOLERANCE=1.e-21;
ENDMATCH;



twiss,sequence=swissfel,range=s20SY02$start/sarcl01$start,beta0=twisssep;
plot,haxis=s,vaxis=betx,bety,range=S30CB15$START/SARCL02$END,colour=100;





!-----------------
muAR=0.16;
Print,text='MATCH POINT: FODO - Aramis Undulator';
match,SEQUENCE=swissfel,range=sarun02$start/sarun03$end;
Vary,name=sarun02.mqua080.k1,step=0.0001;
Vary,name=sarun03.mqua080.k1,step=0.0001;
constraint,sequence=swissfel,range=sarun03$end,mux=muAR;
constraint,sequence=swissfel,range=sarun03$end,muy=muAR;
lmdif,calls=100,tolerance=1e-21;
endmatch;

twiss,sequence=swissfel,range=sarun02$start/sarun03$end;
mbetx=table(twiss,sarun02$end,betx);
malfx=table(twiss,sarun02$end,alfx);
mbety=table(twiss,sarun02$end,bety);
malfy=table(twiss,sarun02$end,alfy);

plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;
sarun04.mqua080.k1=sarun02.mqua080.k1;
sarun05.mqua080.k1=sarun03.mqua080.k1;
sarun06.mqua080.k1=sarun02.mqua080.k1;
sarun07.mqua080.k1=sarun03.mqua080.k1;
sarun08.mqua080.k1=sarun02.mqua080.k1;
sarun09.mqua080.k1=sarun03.mqua080.k1;
sarun10.mqua420.k1=sarun02.mqua080.k1;
sarun11.mqua080.k1=sarun03.mqua080.k1;
sarun12.mqua080.k1=sarun02.mqua080.k1;
sarun13.mqua080.k1=sarun03.mqua080.k1;
sarun14.mqua080.k1=sarun02.mqua080.k1;
sarun15.mqua080.k1=sarun03.mqua080.k1;
sarun16.mqua080.k1=sarun02.mqua080.k1;
sarun17.mqua080.k1=sarun03.mqua080.k1;
sarun18.mqua080.k1=sarun02.mqua080.k1;
sarun19.mqua080.k1=sarun03.mqua080.k1;

SARUN01.mqua080.k1=sarun03.mqua080.k1;
sarma02.mqua050.k1=0;
sarma02.mqua120.k1=0.7*sarun02.mqua080.k1;


sarma01.mqua010.k1 := -1.435783623;
sarma01.mqua060.k1 := 1.29993551;
sarma01.mqua080.k1 := -1.1356922498;
sarma01.mqua120.k1 := 0.4142170392;
sarma01.mqua140.k1 := -0.4428419456;

Print,text='MATCH POINT: Match into Aramis Undulator';
MATCH,SEQUENCE=swissfel,range=S30SY01$start/sarun02$start,beta0=twisspor;
VARY,NAME=sarma01.mqua010.k1,STEP=0.0001;
VARY,NAME=sarma01.mqua060.k1,STEP=0.0001;
VARY,NAME=sarma01.mqua080.k1,STEP=0.0001;
VARY,NAME=sarma01.mqua120.k1,STEP=0.0001;
VARY,NAME=sarma01.mqua140.k1,STEP=0.0001;
CONSTRAINT,SEQUENCE=swissfel,range=sarun01$start,betx=mbetx;
CONSTRAINT,SEQUENCE=swissfel,range=sarun01$start,bety=mbety;
CONSTRAINT,SEQUENCE=swissfel,range=sarun01$start,alfx=malfx;
CONSTRAINT,SEQUENCE=swissfel,range=sarun01$start,alfy=malfy;
CONSTRAINT,SEQUENCE=swissfel,range=sarun02$start,x<5;
LMDIF,CALLS=100,TOLERANCE=1.e-21;
ENDMATCH;

! from eduard
SARUN18.MQUA080.K1=  -6.305232722687153e-01;
SARUN19.MQUA080.K1=   1.167256186348939e+00;
SARUN20.mqua220.K1=   8.363412365675689e-01;
SARBD01.MQUA010.K1=  -2.099999992065428e+00;
SARBD01.MQUA020.K1=  -2.099998228837625e+00;
SARBD02.MQUA030.K1=  -2.099998942347312e+00;

!The matching routine imposes bx=50 and ax=0 at the TDS,
!sin(mux)~1 between the TDS and the dump screen,
!and minimization of by/dispy^2 at the dump screen.
!I find a solution fulfilling the conditions with by=1.136313m  and dispy= 2.534262e-01 m at the screen.

twiss,SEQUENCE=swissfel,range=#s/sarbd02.dscr050$end,beta0=twiss0;
plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;
plot,haxis=s,vaxis=dx,dy,range=#s/#e,colour=100;