diff --git a/Scripts/Reference-SwissFEL/initTwiss.madx b/Scripts/Reference-SwissFEL/initTwiss.madx
new file mode 100644
index 0000000..cc33b02
--- /dev/null
+++ b/Scripts/Reference-SwissFEL/initTwiss.madx
@@ -0,0 +1,9 @@
+option,-echo;
+! starting condition
+Twiss0: beta0, betx = 29.971, alfx = 0.003, bety = 26., alfy = -0.288; ! location:Start
+Twisssep: beta0,betx = 48.26, alfx = 7.322, bety = 14.293, alfy = -3.513;  ! location:sinsy02$start
+TwissM1: beta0, betx = 11.2, alfx = 3.2,  bety = 0.8, alfy = 0.75; ! location:sinlh02.mqua410$end
+TwissM2: beta0, betx = 11.1557, alfx = -1.17,  bety = 50, alfy = 0.; ! location:sindi02.mqua020$start
+TwissM3: beta0, betx = 4.88, alfx = 0.5546,  bety = 16.11, alfy = -1.81; ! location:s10bc01.mqua020$start
+TwissM4: beta0, betx = 6.14, alfx = -0.83,  bety = 22.86, alfy = -1.18; ! location:s10ma01.mqua020$start
+
diff --git a/Scripts/Reference-SwissFEL/matchAramis.madx b/Scripts/Reference-SwissFEL/matchAramis.madx
new file mode 100644
index 0000000..29b1359
--- /dev/null
+++ b/Scripts/Reference-SwissFEL/matchAramis.madx
@@ -0,0 +1,354 @@
+
+
+
+
+
+y_sep=0.01;
+py_sep=0.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;
+
+
+
+use,sequence=swissfel;
+select,flag=Error,pattern="SAT.*";
+ealign,DY=0.01;
+
+MATCH,SEQUENCE=swissfel,range=s20sy02$start/s20sy02$end,beta0=Twisssep;
+VARY,NAME=s20sy02.mqua070.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=s20sy02.mqua100.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=s20sy02.mqua140.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=s20sy02.mqua180.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=SYKICK,STEP=0.0001;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,y=y_sep;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,py=py_sep;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,dy=-0.007;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,dpy=1.7e-3;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mqua180, BETY < 60;
+LMDIF,CALLS=8000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+! this statements are needed to have explicit values in the output file
+SYKICK=0;
+S20SY02.MKDC010.c1:= SYKICK*2;
+S20SY02.MKAC020.c1:= SYKICK*3;
+S20SY02.MKDC030.c1:= SYKICK*2;
+S20SY02.MKAC040.c1:= SYKICK*3;
+S20SY02.MKDC050.c1:= SYKICK*2;
+
+
+s20bc01.mqua010.k1:=1.;
+
+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.mqua050.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;
+
+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;
+
+
+
+
+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;
+
+
+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;
+
+
+!L = 9;
+!alpha = 1.5;
+!beta = alpha*L;
+
+
+!MATCH,SEQUENCE=Swissfel,range=s20sy02$start/s30cb14$end,beta0=twisssep;
+!VARY,NAME=s30cb10.mqua430.k1,STEP=0.0001;
+!VARY,NAME=s30cb11.mqua430.k1,STEP=0.0001;
+!VARY,NAME=s30cb12.mqua430.k1,STEP=0.0001;
+!VARY,NAME=s30cb13.mqua430.k1,STEP=0.0001;
+!VARY,NAME=s30cb14.mqua430.k1,STEP=0.0001;
+!CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb14.mqua430$end,betx=beta;
+!CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb14.mqua430$end,bety=beta;
+!CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb14.mqua430$end,alfx=alpha;
+!CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb14.mqua430$end,alfy=alpha;
+!CONSTRAINT,SEQUENCE=SwissFEL,range= s30cb14.mqua430$end,X<1;
+!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;
+
+
+
+Select,flag=Error,pattern="SPO.*";
+EALIGN,DY=0.01;
+
+
+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;
+
+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;
+
+
+pbetx=3.9036;
+palfx=-0.05021;
+pbety=37.6374;
+palfy=4.4433;
+
+twisspor: beta0, betx=pbetx,alfx=palfx, bety=pbety,alfy=palfy;
+
+s30cb13.mqua430.k1 = 0.456;
+s30cb14.mqua430.k1 = -1.741;
+s30cb15.mqua010.k1 = 1.1533;
+s30cb15.mqua020.k1 = -0.0961;
+
+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=pbetx;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end,bety=pbety;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end,alfx=palfx;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=s30cb15$end,alfy=palfy;
+!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;
+R56 = -000e-6;
+
+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;
+
+MATCH,SEQUENCE=swissfel,range=sarcl02$start/sarma01.mqua010,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/sarma01.mqua010,betx=bx,alfx=ax,bety=by,alfy=ay;
+
+nbetx= table(twiss,sarma01$start,betx);
+nalfx=-table(twiss,sarma01$start,alfx);
+nbety= table(twiss,sarma01$start,bety);
+nalfy=-table(twiss,sarma01$start,alfy);
+
+
+
+
+MATCH,SEQUENCE=swissfel,range=S30SY01$start/sarma01.mqua010,beta0=twisspor;
+VARY,NAME=sarcl01.mqua020.k1,STEP=0.0001;
+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;
+!CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua080,betx=beta;
+!CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua080,alfx=alpha;
+!CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua080,bety=beta;
+!CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua080,alfy=alpha;
+
+CONSTRAINT,SEQUENCE=swissfel,range=sarcl02$start,betx=nbetx;
+CONSTRAINT,SEQUENCE=swissfel,range=sarcl02$start,alfx=nalfx;
+CONSTRAINT,SEQUENCE=swissfel,range=sarcl02$start,bety=nbety;
+CONSTRAINT,SEQUENCE=swissfel,range=sarcl02$start,alfy=nalfy;
+CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua050,bety<30;
+CONSTRAINT,SEQUENCE=swissfel,range=sarcl01.mqua190,betx<30;
+LMDIF,CALLS=100,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+muAR=0.17;
+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;
+
+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;
+
+
+
+twiss,SEQUENCE=swissfel,range=#s/sarun19$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;
+
+exit;
+
+
+
diff --git a/Scripts/Reference-SwissFEL/matchAthos.madx b/Scripts/Reference-SwissFEL/matchAthos.madx
new file mode 100644
index 0000000..840b19a
--- /dev/null
+++ b/Scripts/Reference-SwissFEL/matchAthos.madx
@@ -0,0 +1,396 @@
+
+!-----------------------------------
+! switch yard
+! matching condition at the septum, the new design by Natalia on 27th July 2011
+! find the correct kicker strength to have the right dispersion and offset.
+
+y_sep=0.01;
+py_sep=0.0;
+dy_sep=-y_sep;
+dpy_sep=-py_sep;
+
+
+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;
+
+
+
+use,sequence=swissfel;
+MATCH,SEQUENCE=swissfel,range=s20sy02$start/s20sy02$end,beta0=twisssep;
+VARY,NAME=s20sy02.mqua070.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=s20sy02.mqua100.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=s20sy02.mqua140.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=s20sy02.mqua180.k1,STEP=0.0001,lower=-3,upper=3;
+VARY,NAME=SYKICK,STEP=0.0001;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,y=y_sep;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,py=py_sep;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,dy<-0.006;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mbnd200,dpy<1;
+CONSTRAINT,SEQUENCE=swissfel,range=S20SY02.mqua180, BETY < 60;
+LMDIF,CALLS=8000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+use,sequence=swissfel;
+select,flag=Error,pattern="SAT.*";
+ealign,DY=0.01;
+
+satsy01.mqua020.k1 := 0.48577;
+satsy01.mqua040.k1 := 1.0122;
+satsy01.mqua070.k1 := -0.055;
+satsy01.mqua090.k1 := -0.826;
+
+kmin=-2.1;
+kmax=2.1;
+MATCH,SEQUENCE=swissfel,range=s20SY02$start/satcl01$end,beta0=twisssep;
+VARY, NAME=satsy01.mqua020.k1, STEP=0.00001,LOWER=0,UPPER=4;
+VARY, NAME=satsy01.mqua040.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy01.mqua070.k1, STEP=0.00001,LOWER=Kmin,UPPER=KMAX;
+VARY, NAME=satsy01.mqua090.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01.mbnd200, DX < 7.05e2;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01.mbnd200, DPX =  10.e-3;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01.mbnd200, ALFX=0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01.mbnd200, ALFY=0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01.mbnd200, DPY =-1.16e-4;
+!CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01.mbnd200, DPY =0;
+
+LMDIF, CALLS=1000, TOLERANCE=1.E-21;
+ENDMATCH;
+
+satsy01.mqua210.k1=satsy01.mqua090.k1;
+satsy01.mqua230.k1=satsy01.mqua070.k1;
+satsy01.mqua260.k1=satsy01.mqua040.k1;
+satsy01.mqua280.k1=satsy01.mqua020.k1;
+satsy01.mqua300.k1 := 5.27400e-01;
+
+
+
+
+MATCH,SEQUENCE=swissfel,range=s20SY02$start/satcl01$end,beta0=twisssep;
+VARY, NAME=satsy01.mqua230.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy01.mqua210.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy01.mqua300.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy01.mqua280.k1, STEP=0.00001,LOWER=-4,UPPER=4;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01$end, DX = 0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01$end, DPX = 0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01$end, DY < -0.008;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy01$end, X<10;
+LMDIF, CALLS=1000, TOLERANCE=1.E-21;
+ENDMATCH;
+
+leakdisp = 0;
+satsy01.mqua280.k1 = satsy01.mqua280.k1*(1+leakdisp);
+
+!------------------------------------------
+! close vertical dispersion
+
+
+satsy02.mqua010.k1 := - 1.71;
+satsy02.mqua110.k1 := 1.22622;
+satsy02.mqua120.k1 := 1.64073;
+satsy02.mqua230.k1 = -1.61620;
+
+
+MATCH,SEQUENCE=swissfel,range=s20SY02$start/satcl01$end,beta0=twisssep;
+VARY, NAME=satsy02.mqua010.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy02.mqua110.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy02.mqua120.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy02$end, DY = 0.0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy02$end, DPY = 0.0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satsy02$end, X <0;
+LMDIF, CALLS=1000, TOLERANCE=1.E-21;
+ENDMATCH;
+
+
+satcl011:=-1.3;
+satcl012:= 1.1;
+satcl01.mqua120.k1:=satcl011;
+satcl01.mqua130.k1:=satcl012;
+satcl01.mqua180.k1:=satcl012;
+satcl01.mqua190.k1:=satcl011;
+
+
+
+MATCH, SEQUENCE=swissfel, range=SATCL01$START/SATCL01$END,betx=5,alfx=0,bety=30,alfy=0;
+VARY, NAME=satcl011, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satcl012, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satcl01$end, DX = 0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satcl01$end, DPX = 0;
+LMDIF, CALLS=1000, TOLERANCE=1.E-21;
+ENDMATCH;
+
+
+satsy02.mqua230.k1 = -0.12;
+satsy03.mqua010.k1 = 0.56;
+satsy03.mqua040.k1 = 0.72;
+satsy03.mqua070.k1 = 0.427;
+satsy03.mqua100.k1 = -0.01;
+satsy03.mqua130.k1 = -0.67;
+
+
+MATCH,SEQUENCE=swissfel,range=#s/#e,range=s20SY02$start/satcl01$end,beta0=twisssep;
+VARY, NAME=satsy02.mqua230.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy03.mqua010.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy03.mqua040.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy03.mqua070.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy03.mqua100.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME=satsy03.mqua130.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+CONSTRAINT, SEQUENCE=swissFEL, RANGE=satcl01.msex160,betx=10.8;
+CONSTRAINT, SEQUENCE=swissFEL, RANGE=satcl01.msex160,alfx=-0.01;
+CONSTRAINT, SEQUENCE=swissFEL, RANGE=satcl01.msex160,bety<30.8;  ! was 13.8
+CONSTRAINT, SEQUENCE=swissFEL, RANGE=satcl01.msex160,alfy=-0.01;
+constraint, sequence=swissfel, range=satsy03.mqua130$end, betx <90 ;
+constraint, sequence=swissfel, range=satsy03.mqua130$end, bety <90 ;
+LMDIF, CALLS=1000, TOLERANCE=1.E-21;
+ENDMATCH;
+
+
+MATCH, SEQUENCE=swissfel, range=s20SY02$start/satcl01$end,beta0=twisssep;
+VARY, NAME= satcl01.mqua180.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+VARY, NAME= satcl01.mqua190.k1, STEP=0.00001,LOWER=Kmin,UPPER=Kmax;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satcl01$end, DX = 0;
+CONSTRAINT, SEQUENCE=swissfel, RANGE=satcl01$end, DPX = 0;
+LMDIF, CALLS=1000, TOLERANCE=1.E-21;
+ENDMATCH;
+
+twiss,sequence=swissfel,range=s20sy02$start/satcl01$end,rmatrix,beta0=twisssep;
+
+plot,haxis=s,vaxis=betx,bety,range=s20SY02$start/#e,colour=100;
+plot,haxis=s,vaxis=dx,dy,range=s20SY02$start/#e,colour=100;
+plot,haxis=s,vaxis=re56,range=s20SY02$start/#e,colour=100;
+
+
+
+!---------------------------------------
+! match first cell of linac1RF cell of Athos
+muATL := 0.053;
+
+! Linac 1 Fodo lattice matching
+satcb01.mqua230.k1 := 1.490961973;
+satcb01.mqua430.k1 := -1.490511125;
+
+
+match,SEQUENCE=swissfel,range=satcb01$start/satcb01$end;
+Vary,name=SatCB01.MQUA230.k1,step=0.0001;
+Vary,name=SatCB01.MQUA430.k1,step=0.0001;
+constraint,sequence=satcb01,range=#e,mux=muATL;
+constraint,sequence=satcb01,range=#e,muy=muATL;
+lmdif,calls=100,tolerance=1e-21;
+endmatch;
+
+
+satcl02.mqua230.k1=satcb01.mqua230.k1;
+satcl02.mqua430.k1=satcb01.mqua430.k1;
+
+
+twiss,SEQUENCE=swissfel,range=satcb01$start/satcb01$end;
+nbetx=table(twiss,satcb01$end,betx);
+nalfx=table(twiss,satcb01$end,alfx);
+nbety=table(twiss,satcb01$end,bety);
+nalfy=table(twiss,satcb01$end,alfy);
+
+
+!---------------------------------------
+! match to reference point at satdi01-mqua250
+
+MATCH, SEQUENCE=swissfel,range=s20sy02$start/satcb01$start,beta0=twisssep;
+VARY,NAME=satdi01.mqua025.k1,STEP=0.0001,lower=-2.0,upper=2.0;
+VARY,NAME=satdi01.mqua080.k1,STEP=0.0001,lower=-2.0,upper=2.0;
+VARY,NAME=satdi01.mqua220.k1,STEP=0.0001,lower=-1.8,upper=1.8;
+VARY,NAME=satdi01.mqua230.k1,STEP=0.0001,lower=-1.8,upper=1.8;
+CONSTRAINT,SEQUENCE=swissfel,range=satdi01.mqua250$start,betx=35.76;
+CONSTRAINT,SEQUENCE=swissfel,range=satdi01.mqua250$start,alfx=-1.166;
+CONSTRAINT,SEQUENCE=swissfel,range=satdi01.mqua250$start,bety=24.087;
+CONSTRAINT,SEQUENCE=swissfel,range=satdi01.mqua250$start,alfy=0.648;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+!------------------------------------
+! match into SATCB01
+
+MATCH, SEQUENCE=swissfel,range=s20sy02$start/satcb01$end,beta0=twisssep;
+VARY,NAME=satdi01.mqua250.k1,STEP=0.0001,lower=-1.8,upper=1.8;
+VARY,NAME=satdi01.mqua260.k1,STEP=0.0001,lower=-1.8,upper=1.8;
+VARY,NAME=satdi01.mqua280.k1,STEP=0.0001,lower=-1.8,upper=1.8;
+VARY,NAME=satdi01.mqua300.k1,STEP=0.0001,lower=-1.8,upper=1.8;
+
+CONSTRAINT,SEQUENCE=swissfel,range=satcb01$start,betx=nbetx;
+CONSTRAINT,SEQUENCE=swissfel,range=satcb01$start,alfx=nalfx;
+CONSTRAINT,SEQUENCE=swissfel,range=satcb01$start,bety=nbety;
+CONSTRAINT,SEQUENCE=swissfel,range=satcb01$start,alfy=nalfy;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+twiss, SEQUENCE=swissFEL,range=s20sy02$start/satcl02$end,beta0=twisssep;
+plot,haxis=s,vaxis=betx,bety,range=satdi01$start/satcl02$end,colour=100;
+
+
+
+!!!!!!!!!!!!!!!!!!!!!! Athos undulator
+
+pol=-1;
+satun06.mqua080.k1 := -1.699673315*pol;
+satun07.mqua080.k1 := 1.746074065*pol;
+muAT=0.15;
+match,SEQUENCE=swissfel,range=satun06$start/satun07$end;
+Vary,name=satun06.mqua080.k1,step=0.0001;
+Vary,name=satun07.mqua080.k1,step=0.0001;
+constraint,sequence=swissfel,range=satun07$end,mux=muAT;
+constraint,sequence=swissfel,range=satun07$end,muy=muAT;
+lmdif,calls=100,tolerance=1e-21;
+endmatch;
+
+
+satun08.mqua080.k1 := satun06.mqua080.k1;
+satun09.mqua080.k1 := satun07.mqua080.k1;
+satun10.mqua080.k1 := satun06.mqua080.k1;
+satun11.mqua080.k1 := satun07.mqua080.k1;
+satun12.mqua080.k1 := satun06.mqua080.k1;
+satun13.mqua080.k1 := satun07.mqua080.k1;
+satun14.mqua080.k1 := satun06.mqua080.k1;
+satun15.mqua080.k1 := satun07.mqua080.k1;
+satun16.mqua420.k1 := satun06.mqua080.k1;
+satun17.mqua080.k1 := satun07.mqua080.k1;
+satun18.mqua080.k1 := satun06.mqua080.k1;
+satun19.mqua080.k1 := satun07.mqua080.k1;
+satun20.mqua080.k1 := satun06.mqua080.k1;
+satun21.mqua080.k1 := satun07.mqua080.k1;
+satun22.mqua080.k1 := satun06.mqua080.k1;
+
+
+twiss,sequence=swissfel,range=satun06$start/satun07$end;
+plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;
+ATbetx=table(twiss,satun07$end,betx);
+ATalfx=table(twiss,satun07$end,alfx);
+ATbety=table(twiss,satun07$end,bety);
+ATalfy=table(twiss,satun07$end,alfy);
+
+
+
+eebetx=8;
+eealfx=0.5;
+eebety=30;
+eealfy=-0.5;
+
+eebetx=7;
+eealfx=0.5;
+eebety=30;
+eealfy=1.5;
+
+twissee: beta0,betx=eebetx,bety=eebety,alfx=eealfx,alfy=eealfy;
+
+MATCH, SEQUENCE=swissfel,range=satma01.mqua250$end/satun07$end,beta0=twissee;
+VARY,NAME=satma01.mqua610.k1,STEP=0.0001,lower=-2.0,upper=2.0;
+VARY,NAME=satma01.mqua630.k1,STEP=0.0001,lower=-2.0,upper=2.0;
+VARY,NAME=satma01.mqua640.k1,STEP=0.0001,lower=-2.0,upper=2.0;
+VARY,NAME=satun04.mqua060.k1,STEP=0.0001,lower=-3.8,upper=3.8;
+VARY,NAME=satun05.mqua420.k1,STEP=0.0001,lower=-3.8,upper=3.8;
+CONSTRAINT,SEQUENCE=swissfel,range=satun06$start,betx=ATbetx;
+CONSTRAINT,SEQUENCE=swissfel,range=satun06$start,bety=ATbety;
+CONSTRAINT,SEQUENCE=swissfel,range=satun06$start,alfx=ATalfx;
+CONSTRAINT,SEQUENCE=swissfel,range=satun06$start,alfy=ATalfy;
+CONSTRAINT,SEQUENCE=swissfel,range=satun06$start,x<10;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+twiss,SEQUENCE=swissfel,range=satma01.mqua250$end/satun11$end,beta0=twissee;
+plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;
+
+
+
+scl=4;
+satcb02.mqua230.k1=scl*0.2;
+satcb02.mqua430.k1=-scl*0.2;
+satcb03.mqua230.k1=scl*0.2;
+satcb03.mqua430.k1=-scl*0.2;
+
+
+
+satma01.mqua120.k1=0;
+satma01.mqua140.k1=0;
+satma01.mqua160.k1=0;
+satma01.mqua180.k1=0;
+satma01.mqua230.k1=0;
+satma01.mqua250.k1=0;
+
+
+MATCH, SEQUENCE=swissfel,range=s20sy02$start/satma01$end,beta0=twisssep;
+VARY,NAME=satcb03.mqua430.k1,STEP=0.0001,lower=-3.8,upper=3.8;
+VARY,NAME=satma01.mqua120.k1,STEP=0.0001,lower=-3.8,upper=3.8;
+VARY,NAME=satma01.mqua140.k1,STEP=0.0001,lower=-3.8,upper=3.8;
+VARY,NAME=satma01.mqua160.k1,STEP=0.0001,lower=-3.8,upper=3.8;
+VARY,NAME=satma01.mqua180.k1,STEP=0.0001,lower=-3.8,upper=3.8;
+VARY,NAME=satma01.mqua230.k1,STEP=0.0001,lower=-3.0,upper=3.0;
+VARY,NAME=satma01.mqua250.k1,STEP=0.0001,lower=-3.0,upper=3.0;
+CONSTRAINT,SEQUENCE=swissfel,range=satma01.mqua250$end,betx=eebetx;
+CONSTRAINT,SEQUENCE=swissfel,range=satma01.mqua250$end,bety=eebety;
+CONSTRAINT,SEQUENCE=swissfel,range=satma01.mqua250$end,alfx=eealfx;
+CONSTRAINT,SEQUENCE=swissfel,range=satma01.mqua250$end,alfy=eealfy;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua120$start, betx <70;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua120$start, bety <70;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua140$start, betx <70;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua140$start, bety <70;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua160$start, betx <70;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua160$start, bety <70;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua180$start, betx <70;
+CONSTRAINT,SEQUENCE=swissFEL,range=satma01.mqua180$start, bety <70;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+MATCH, SEQUENCE=swissfel,range=s20sy02$start/satbd01$end,beta0=twisssep;
+VARY,NAME=satun22.mqua080.k1,STEP=0.0001,lower=-5.0,upper=5.0;
+VARY,NAME=satma02.mqua010.k1,STEP=0.0001,lower=-5.0,upper=5.0;
+VARY,NAME=satma02.mqua020.k1,STEP=0.0001,lower=-5.0,upper=5.0;
+VARY,NAME=satma02.mqua040.k1,STEP=0.0001,lower=-5.0,upper=5.0;
+VARY,NAME=satma02.mqua050.k1,STEP=0.0001,lower=-5.0,upper=5.0;
+VARY,NAME=satma02.mqua070.k1,STEP=0.0001,lower=-5.0,upper=5.0;
+
+CONSTRAINT,SEQUENCE=swissfel,range=satbd01.mqua010$start,betx=50;
+CONSTRAINT,SEQUENCE=swissfel,range=satbd01.mqua010$start,bety=50;
+CONSTRAINT,SEQUENCE=swissfel,range=satbd01.mqua010$start,alfx=0;
+CONSTRAINT,SEQUENCE=swissfel,range=satbd01.mqua010$start,alfy=0;
+CONSTRAINT,SEQUENCE=swissfel,range=satbd01.mqua010$start,x=0;
+CONSTRAINT,SEQUENCE=swissfel,range=satbd01.mqua010$start,px=0;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+SATBD01.MQUA010.K1=   5.069789237265826e-01;
+SATBD01.MQUA030.K1=  -6.184665903305938e-01 ;
+SATBD01.MQUA050.K1=   1.324854421128419e-01 ;
+SATBD01.MQUA070.K1=   4.740280173415730e-01 ;
+SATBD01.MQUA090.K1= -7.442648986155268e-01;
+SATBD02.MQUA030.K1=  -1.050806764409536e+00;
+
+twiss,SEQUENCE=swissfel,range=s20sy02$start/satbd01$end,beta0=twisssep;
+plot,haxis=s,vaxis=betx,bety,range=satdi01$start/#e,colour=100;
+
+
+twiss,SEQUENCE=swissfel,range=#s/satbd01$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;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Scripts/Reference-SwissFEL/matchInjector.madx b/Scripts/Reference-SwissFEL/matchInjector.madx
new file mode 100644
index 0000000..e887c8b
--- /dev/null
+++ b/Scripts/Reference-SwissFEL/matchInjector.madx
@@ -0,0 +1,283 @@
+! matching scripts for SwissFEL reference
+
+! into laser heater
+
+MATCH,SEQUENCE=SwissFEL,range=#s/sinbc02$start,beta0 = Twiss0;
+VARY,NAME=sinlh01.MQUA020.k1,STEP=0.0001;
+VARY,NAME=sinlh01.MQUA040.k1,STEP=0.0001;
+VARY,NAME=sinlh01.MQUA050.k1,STEP=0.0001;
+VARY,NAME=sinlh01.MQUA070.k1,STEP=0.0001;
+VARY,NAME=sinlh02.MQUA010.k1,STEP=0.0001;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinlh02.mqua410$end, alfx = TwissM1->ALFX;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinlh02.mqua410$end, alfy = TwissM1->ALFY;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinlh02.mqua410$end, betx = TwissM1->BETX;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinlh02.mqua410$end, bety = TwissM1->BETY;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb04.mqua130$end, X < 10;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+SIMPLEX,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+!---------------------------------------------------------------------------------
+! match the periodic lattice of the S-band structure
+
+
+sinsb04.mqua130.k1 :=0.7219189344;
+sinsb04.mqua230.k1 :=-0.71562666978;
+muInj := 0.2;
+MATCH,SEQUENCE=swissfel,range=sinsb04$start/sinsb04$end;
+VARY,NAME=sinsb04.MQUA130.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=sinsb04.MQUA230.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+CONSTRAINT,SEQUENCE=swissfel,range=#e,mux=muInj;
+CONSTRAINT,SEQUENCE=swissfel,range=#e,muy=muInj;
+LMDIF,CALLS=100,TOLERANCE=1.e-21;
+ENDMATCH;
+
+twiss,SEQUENCE=swissfel,range=sinsb04$start/sinsb04$end;
+nbetx=table(twiss,sinsb04$end,betx);
+nalfx=table(twiss,sinsb04$end,alfx);
+nbety=table(twiss,sinsb04$end,bety);
+nalfy=table(twiss,sinsb04$end,alfy);
+
+! Matching from Laser heater into S-band linac
+
+
+
+MATCH,SEQUENCE=SwissFEL,range=#s/sinsb04$end,beta0 = Twiss0;
+VARY,NAME=sinlh02.MQUA410.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinlh03.MQUA030.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinlh03.MQUA040.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinlh03.MQUA060.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinlh03.MQUA080.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinsb03.MQUA130.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinsb03.MQUA230.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb04$start,betx=nbetx;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb04$start,bety=nbety;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb04$start,alfx=nalfx;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb04$start,alfy=nalfy;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb03.mqua130$start,betx<50;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb03.mqua130$start,bety<50;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb03.mqua230$start,betx<50;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinsb03.mqua230$start,bety<50;
+
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+twiss,sequence=SwissFEL,range=#s/sinsb04$end,beta0 = twiss0;
+plot,haxis=s,vaxis=betx,bety,colour=100;
+
+sinsb05.mqua130.k1 := sinsb04.mqua130.k1;
+sinsb05.mqua230.k1 := sinsb04.mqua230.k1;
+
+MATCH,SEQUENCE=SwissFEL,range=#s/sindi02$end,beta0=twiss0;
+VARY,NAME=sinbc01.mqua020.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinbc01.mqua050.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinbc01.mqua070.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinbc01.mqua090.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sinbc01.mqua110.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sindi01.mqua020.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+VARY,NAME=sindi01.mqua070.k1,STEP=0.0001,lower=-3.5,upper=3.5;
+
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinbc01.mqua110$end,betx<100;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinbc02.mbnd400,bety<100;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sinbc02.mbnd400,betx<5.1;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sindi02.mqua020$start,bety=TwissM2->BETY;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sindi02.mqua020$start,alfy=TwissM2->ALFY;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sindi02.mqua020$start,betx=TwissM2->BETX;
+CONSTRAINT,SEQUENCE=SwissFEL,range=sindi02.mqua020$start,alfx=TwissM2->ALFY;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+twiss,sequence=SwissFEL,range=#s/SINDI02.MQUA020$START,beta0=twiss0;
+plot,haxis=s,vaxis=betx,bety,colour=100;
+
+
+
+
+
+
+!-------------------------------------
+! match periodic lattice of Linac 1
+
+muLin1 := 0.1883;
+
+s10cb02.mqua230.k1 := -1.490961973;
+s10cb02.mqua430.k1 := 1.490511125;
+
+MATCH,SEQUENCE=swissfel,range=s10cb02$start/s10cb02$end;
+VARY,NAME=s10cb02.MQUA130.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10cb02.MQUA230.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+CONSTRAINT,SEQUENCE=swissfel,range=#e,mux=muLin1;
+CONSTRAINT,SEQUENCE=swissfel,range=#e,muy=muLin1;
+LMDIF,CALLS=100,TOLERANCE=1.e-21;
+ENDMATCH;
+
+twiss,SEQUENCE=swissfel,range=s10cb02$start/s10cb02$end;
+nbetx=table(twiss,s10cb02$end,betx);
+nalfx=table(twiss,s10cb02$end,alfx);
+nbety=table(twiss,s10cb02$end,bety);
+nalfy=table(twiss,s10cb02$end,alfy);
+
+
+s10cb01.mqua230.k1 := s10cb02.mqua230.k1;
+s10cb01.mqua430.k1 := s10cb02.mqua430.k1;
+S10DI01.MQUA120.k1 :=S10CB02.MQUA230.k1;
+S10CB03.MQUA230.k1 :=S10CB02.MQUA430.k1;
+S10CB03.MQUA430.k1 :=S10CB02.MQUA230.k1;
+S10CB04.MQUA230.k1 :=S10CB02.MQUA430.k1;
+S10CB04.MQUA430.k1 :=S10CB02.MQUA230.k1;
+S10CB05.MQUA230.k1 :=S10CB02.MQUA430.k1;
+S10CB05.MQUA430.k1 :=S10CB02.MQUA230.k1;
+S10CB06.MQUA230.k1 :=S10CB02.MQUA430.k1;
+S10CB06.MQUA430.k1 :=S10CB02.MQUA230.k1;
+S10CB07.MQUA230.k1 :=S10CB02.MQUA430.k1;
+S10CB07.MQUA430.k1 :=S10CB02.MQUA230.k1;
+S10CB08.MQUA230.k1 :=S10CB02.MQUA430.k1;
+S10CB08.MQUA430.k1 :=S10CB02.MQUA230.k1;
+S10CB09.MQUA230.k1 :=S10CB02.MQUA430.k1;
+
+
+
+!-----------------------------------
+! match into Linac 1st
+
+
+
+MATCH,SEQUENCE=SwissFEL,range=#s/s10cb02$end,beta0=twiss0;
+VARY,NAME=sindi02.mqua020.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=sindi02.mqua030.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=sindi02.mqua050.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=sindi02.mqua060.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=sindi02.mqua090.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=s10cb01.mqua230.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=s10cb01.mqua430.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+
+CONSTRAINT,SEQUENCE=swissfel,range=s10cb02$start,betx=nbetx;
+CONSTRAINT,SEQUENCE=swissfel,range=s10cb02$start,alfx=nalfx;
+CONSTRAINT,SEQUENCE=swissfel,range=s10cb02$start,bety=nbety;
+CONSTRAINT,SEQUENCE=swissfel,range=s10cb02$start,alfy=nalfy;
+CONSTRAINT,SEQUENCE=swissfel,range=s10cb01$start,x=0;
+CONSTRAINT,SEQUENCE=swissfel,range=s10cb01$start,y=0;
+CONSTRAINT,SEQUENCE=swissfel,range=s10cb01$start,px=0;
+!CONSTRAINT,SEQUENCE=swissfel,range=s10cb01$start,py=0;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+twiss,sequence=SwissFEL,range=#s/S10BC01$start,beta0=twiss0;
+plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;
+
+
+USE,SEQUENCE=SwissFEL;
+MATCH,SEQUENCE=SwissFEL,range=#s/s20cb01$start,beta0=twiss0;
+VARY,NAME=s10cb07.mqua430.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10cb08.mqua230.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10cb08.mqua430.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10cb09.mqua230.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10bc01.mqua020$start,betx=TwissM3->BETX;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10bc01.mqua020$start,alfx=TwissM3->ALFX;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10bc01.mqua020$start,bety=TwissM3->BETY;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10bc01.mqua020$start,alfy=TwissM3->ALFY;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+!-------------------------------------------------
+! Linac Bunch Compressor
+! reference point s10ma01-mqua020
+
+
+
+USE,SEQUENCE=SwissFEL;
+MATCH,SEQUENCE=SwissFEL,range=#s/s20cb01$start,beta0 = twiss0;
+VARY,NAME=s10bc01.mqua020.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10bc01.mqua040.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10bc01.mqua060.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10bc01.mqua080.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+VARY,NAME=s10bc01.mqua100.k1,STEP=0.0001,lower=-2.5,upper=2.5;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10ma01.mqua020$start,betx=TwissM4->BETX;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10ma01.mqua020$start,alfx=TwissM4->ALFX;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10ma01.mqua020$start,bety=TwissM4->BETY;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10ma01.mqua020$start,alfy=TwissM4->ALFY;
+CONSTRAINT,SEQUENCE=swissFEL,range=s10ma01.mqua020$start,dx=0;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+
+!---------------------------------------
+! Linac 2
+
+
+s20cb01.mqua430.k1 = -0.8156528972;
+s20cb02.mqua430.k1 = 0.8156479879;
+
+muLin2=0.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=linac2,range=#e,mux=muLin2;
+constraint,sequence=linac2,range=#e,muy=muLin2;
+lmdif,calls=100,tolerance=1e-21;
+endmatch;
+
+s20cb03.mqua430.k1 = s20cb01.mqua430.k1;
+s20cb04.mqua430.k1 = s20cb01.mqua430.k1;
+
+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);
+
+
+
+!--------------------------------
+! Linac Bunch Compressor and matching into Linac 2
+
+MATCH,SEQUENCE=swissfel,beta0=twiss0;
+VARY,NAME=s10ma01.mqua020.k1,STEP=0.0001,lower=-2.4,upper=2.4;
+VARY,NAME=s10ma01.mqua050.k1,STEP=0.0001,lower=-2.4,upper=2.4;
+VARY,NAME=s10ma01.mqua070.k1,STEP=0.0001,lower=-2.4,upper=2.4;
+VARY,NAME=s10ma01.mqua110.k1,STEP=0.0001,lower=-2.4,upper=2.4;
+VARY,NAME=s10ma01.mqua130.k1,STEP=0.0001,lower=-2.4,upper=2.4;
+CONSTRAINT,SEQUENCE=Swissfel,range=s20cb01$start,betx=nbetx;
+CONSTRAINT,SEQUENCE=Swissfel,range=s20cb01$start,alfx=nalfx;
+CONSTRAINT,SEQUENCE=Swissfel,range=s20cb01$start,bety=nbety;
+CONSTRAINT,SEQUENCE=Swissfel,range=s20cb01$start,alfy=nalfy;
+CONSTRAINT,SEQUENCE=Swissfel,range=s10ma01.mqua070$start,bety<100;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+
+!-----------------
+! into septum
+
+MATCH,SEQUENCE=SwissFEL,range=#s/s20sy02$start,beta0 = Twiss0;
+VARY,NAME=s20sy01.mqua020.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=s20sy01.mqua030.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=s20sy01.mqua050.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+VARY,NAME=s20sy01.mqua080.k1,STEP=0.0001,lower=-5.5,upper=5.5;
+
+CONSTRAINT,SEQUENCE=swissfel,range=s20sy01$end,betx=Twisssep->BETX;
+CONSTRAINT,SEQUENCE=swissfel,range=s20sy01$end,alfx=Twisssep->ALFX;
+CONSTRAINT,SEQUENCE=swissfel,range=s20sy01$end,bety=Twisssep->BETY;
+CONSTRAINT,SEQUENCE=swissfel,range=s20sy01$end,alfy=Twisssep->ALFY;
+LMDIF,CALLS=1000,TOLERANCE=1.e-21;
+ENDMATCH;
+
+!
+twiss,SEQUENCE=swissfel,range=#s/s20sy01$end,beta0 = Twiss0;
+plot,haxis=s,vaxis=betx,bety,range=s10cb08$start/#e,colour=100;
+plot,haxis=s,vaxis=betx,bety,range=#s/#e,colour=100;
+
+
+
diff --git a/issues.txt b/issues.txt
index 5c3db0d..08e8c80 100644
--- a/issues.txt
+++ b/issues.txt
@@ -1,4 +1,7 @@
-Currently none
+1) Initialize all quadrupoles at the beginning!
+2) Allow for parameter inputs, e.g. phase advance
+3) make initTwiss to alist of several reference points
+4) match to beta0 and not individual elements
 
 
 
diff --git a/matchmaker.py b/matchmaker.py
index c73190e..9eda9ac 100644
--- a/matchmaker.py
+++ b/matchmaker.py
@@ -5,7 +5,7 @@ from onlinemodel.madx import CMadX
 
 class MatchMaker:
     def __init__(self):
-        self.matchlist={'SwissFEL+':'Scripts/SFPlus'}
+        self.matchlist={'Reference-SwissFEL':'Scripts/Reference-SwissFEL','SwissFEL+':'Scripts/SFPlus'}
 
     def initScripts(self,target):
         self.scriptdir = self.matchlist[target]