!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! MAD optics file for the CTF3 linac
! calculate initial conditions from quad scan on girder 10
! 
!
! Frank Tecker
! Mar. 2005
! last modification:
!      17.03.2005  initial version
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!OPTION, -echo;

! outout to file OUTPUT
!assign, echo=output;

skipprematch = 0;
MADXFLAG := 0;
PTCNSTEPS := 10;
ptc_setswitch, debuglevel=0, nocavity=false, 
               maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;


call, file="currents";
nex=nex*1e-6;
ney=ney*1e-6;

call, file="ctf3-2005-1.seqx";





! try an even weaker focusing
k_ls1     := -2.652174957014; ! weaker foc.
k_ls2     :=  5.007697921381; ! weaker foc.



! calculate Twiss parameters for regular cell structure
BEAM, PARTICLE=ELECTRON, ENERGY=E4/1000. , sigt=1, sige=1/1000 ;
USE, period=cell_lin;
SAVEBETA, PLACE=#S, LABEL=REGCELL;
SAVEBETA, PLACE=cell_mrk, LABEL=MIDCELL;
SAVEBETA, PLACE=cell_q1,  LABEL=TriplQ1;
SAVEBETA, PLACE=cell_q3,  LABEL=TriplQ3;
TWISS, SAVE;
value, MIDCELL->BETX;
value, MIDCELL->ALFX;
value, MIDCELL->BETY;
value, MIDCELL->ALFY;

write, table=twiss, file=stdtwiss.10-15.txt;

PLOT, noversion=true, haxis=s, vaxis=betx, bety, colour=100, spline, 
      title="Regular linac cell";



call, file="definitions.madx";

!it is possible to match an user defined variable if it is named mvar1, mvar2, mvar3, or mvar4
!mvar1 := sigx; 
!mvar2 := sigy;
!select, flag=twiss, column=name, s, energy, beta11,beta22,alfa11,alfa22, disp1,disp2, mu1,mu2, mvar1, mvar2;

!select, flag=twiss, column=name, s, energy, beta11,beta22,alfa11,alfa22, 
!                           disp1,disp2, mu1,mu2, 
!	       ptcsigx, ptcsigy, ptcsigp;


select, flag=twiss, clear;
select, flag=twiss, column=name, energy, beta11, alfa11, x, px, t, pt, disp1, disp2;

!****** initial twiss parameters ....
B0_10  :  BETA0, BETX=BET0X10, ALFX=ALF0X10, 
                 BETY=BET0Y10, ALFY=ALF0Y10;

show, B0_10;

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
!Here comes initial


!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

title, "CTF3 Last Run 2005 - Girders 10 to 15";

value, E8;
BEAM, PARTICLE=ELECTRON, ENERGY=E8/1000. , sigt=1, sige=1/1000 ;
USE, period=LINAC10; != ( SEC10Q,  SEC11, SEC12, SEC13, SEC14, SEC15 )


ptc_create_universe;
ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
ptc_setswitch, debuglevel=0, nocavity=false, 
               maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
PTC_TWISS, table=twiss, icase=6, no=2, BETA0=B0_10, betz=1.0;
ptc_end;



print, text="PLOT, noversion=true, haxis=s, vaxis=betx, bety";
PLOT, noversion=true, haxis=s, vaxis=beta11, beta22, colour=100, spline, title="Initial optics (PTC)";



print, text="###################################################################################";
print, text="###################################################################################";
print, text="###################################################################################";
print, text="Match quads at girders 10 and 11 to periodic solution at MID11, starting with B0_10";


USE, period=SEC10Q_11tomid;!: LINE (  SEC10Q, SEC11tm);

!if (skipprematch == 0)
! {
   match, beta0=b0_10;
      vary, name=IQDB1005, step=0.1, lower= -12, upper=12;  !Line SEC10Q, element 2,  currnt of the quad
      vary, name=IQFC1010, step=0.1, lower= -100,upper=100; !Line SEC10Q, element 4,  currnt of the quad
      vary, name=IQDB1015, step=0.1, lower= -12, upper=12;  !Line SEC10Q, element 6,  currnt of the quad

      CONSTRAINT, range=MID10, betx = MIDCELL->BETX, alfx = MIDCELL->ALFX, 
                               bety = MIDCELL->BETY, alfy = MIDCELL->ALFY;

      jacobian, calls=100,tolerance=1e-04;
!      lmdif,    calls=100,tolerance=1e-4;
   !   migrad,  calls=100,tolerance=1e-4;
   endmatch;

   value, table(twiss,MID10,betx) - MIDCELL->BETX;
   value, table(twiss,MID10,alfx) - MIDCELL->ALFX;
   value, table(twiss,MID10,bety) - MIDCELL->BETY;
   value, table(twiss,MID10,alfy) - MIDCELL->ALFY;


   match, use_macro; !beta0=b0_10;

      vary, name=IQDB1005, step=1e-06, lower= -12, upper=12;  !Line SEC10Q, element 2,  currnt of the quad
      vary, name=IQFC1010, step=1e-06, lower= -100,upper=100; !Line SEC10Q, element 4,  currnt of the quad
      vary, name=IQDB1015, step=1e-06, lower= -12, upper=12;  !Line SEC10Q, element 6,  currnt of the quad

      vary, name=IQDB1105, step=1e-02, lower= 0, upper=12;  !Line SEC11 , element 2,  currnt of the quad
      vary, name=IQFC1110, step=1e-02, lower= 0, upper=100; !Line SEC11 , element 4,  currnt of the quad

      m11a: macro =
        {
          TWISS,  beta0=b0_10;

          value, IQDB1005; 
          value, IQFC1010; 
          value, IQDB1015; 
          value, IQDB1105; 
          value, IQFC1110;

          value, table(twiss,MID10,betx) - MIDCELL->BETX;
          value, table(twiss,MID10,alfx) - MIDCELL->ALFX;
          value, table(twiss,MID10,bety) - MIDCELL->BETY;
          value, table(twiss,MID10,alfy) - MIDCELL->ALFY;

          value, table(twiss,MID11,betx) - MIDCELL->BETX; ! Line SEC11, in element D1160
          value, table(twiss,MID11,alfx) - MIDCELL->ALFX;
          value, table(twiss,MID11,bety) - MIDCELL->BETY;
          value, table(twiss,MID11,alfy) - MIDCELL->ALFY;
        };

      CONSTRAINT, weight=1,  expr = table(twiss,MID10,betx) = MIDCELL->BETX;
      CONSTRAINT, weight=10, expr = table(twiss,MID10,alfx) = MIDCELL->ALFX;
      CONSTRAINT, weight=1,  expr = table(twiss,MID10,bety) = MIDCELL->BETY;
      CONSTRAINT, weight=10, expr = table(twiss,MID10,alfy) = MIDCELL->ALFY;

      CONSTRAINT, weight=100,  expr = table(twiss,MID11,betx) = MIDCELL->BETX;
      CONSTRAINT, weight=1000, expr = table(twiss,MID11,alfx) = MIDCELL->ALFX;
      CONSTRAINT, weight=100,  expr = table(twiss,MID11,bety) = MIDCELL->BETY;
      CONSTRAINT, weight=1000, expr = table(twiss,MID11,alfy) = MIDCELL->ALFY;

      jacobian,    calls=100,tolerance=1e-06;
   !   migrad,   calls=3000,tolerance=1e-09;
   endmatch;

   value, MIDCELL->BETX; ! Line SEC11, in element D1160
   value, MIDCELL->ALFX;
   value, MIDCELL->BETY;
   value, MIDCELL->ALFY;

   value, table(twiss,MID11,betx); ! Line SEC11, in element D1160
   value, table(twiss,MID11,alfx);
   value, table(twiss,MID11,bety);
   value, table(twiss,MID11,alfy);

   value, table(twiss,MID11,betx) - MIDCELL->BETX; ! Line SEC11, in element D1160
   value, table(twiss,MID11,alfx) - MIDCELL->ALFX;
   value, table(twiss,MID11,bety) - MIDCELL->BETY;
   value, table(twiss,MID11,alfy) - MIDCELL->ALFY;


! };


value, IQDB1005; 
value, IQFC1010; 
value, IQDB1015; 
value, IQDB1105; 
value, IQFC1110;

value, E8, E11, P11, 100 *FQC / P11;


match, use_ptcknob;

   ptc_varyknob, element=CL.QDB1005, kn=1, exactmatch=true, 
                            lower= -12 *FQB / P10, 
                            upper=  12 *FQB / P10, 
	        trustrange = 0.2, step=1e-6 ;
   ptc_varyknob, element=CL.QFC1010, kn=1, exactmatch=true, !Line SEC10Q, element 4,  currnt of the quad
                            lower= -100 *FQC / P10, 
                            upper=  100 *FQC / P10, 
	        trustrange = 0.2, step=1e-6 ;
   ptc_varyknob, element=CL.QDB1015, kn=1, exactmatch=true, !Line SEC10Q, element 6,  currnt of the quad
                            lower= -12 *FQB / P10, 
                            upper=  12 *FQB / P10, 
	        trustrange = 0.1, step=1e-6 ;

   ptc_varyknob, element=CL.QDB11, kn=1, exactmatch=false, !Line SEC11 , element 2,  currnt of the quad
                            lower = -12 *FQB / P11, 
                            upper =  0, 
	        trustrange = 0.3, step=1e-6 ;

   ptc_varyknob, element=CL.QFC1110, kn=1, exactmatch=true, !Line SEC11 , element 4,  currnt of the quad
                            lower=  0, 
                            upper=  100 *FQC / P11, 
	        trustrange = 0.3 , step=1e-6;

   ptc_create_universe;
   ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
   ptc_setswitch, debuglevel=0, nocavity=false, 
                  maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
   PTC_TWISS, table=twiss, icase=6, no=3, beta0=b0_10, betz=1.0;



   CONSTRAINT, weight=1,  expr = table(twiss,MID10,beta11) = MIDCELL->BETX;
   CONSTRAINT, weight=10, expr = table(twiss,MID10,alfa11) = MIDCELL->ALFX;
   CONSTRAINT, weight=1,  expr = table(twiss,MID10,beta22) = MIDCELL->BETY;
   CONSTRAINT, weight=10, expr = table(twiss,MID10,alfa22) = MIDCELL->ALFY;

   CONSTRAINT, weight=100,  expr = table(twiss,MID11,beta11) = MIDCELL->BETX;
   CONSTRAINT, weight=100, expr = table(twiss,MID11,alfa11) = MIDCELL->ALFX;
   CONSTRAINT, weight=100,  expr = table(twiss,MID11,beta22) = MIDCELL->BETY;
   CONSTRAINT, weight=100, expr = table(twiss,MID11,alfa22) = MIDCELL->ALFY;

!   jacobian,    calls=100, tolerance=1e-04;
   lmdif,    calls=100, tolerance=1e-10;
endmatch; 

IQDB1005 = -cl.qdb1005_kn1  * P10 / FQB;
IQFC1010 =  cl.qfc1010_kn1  * P10 / FQC;
IQDB1015 = -cl.qdb1015_kn1  * P10 / FQB;
IQDB1105 = -cl.qdb11_kn1    * P11 / FQB;
IQFC1110 =  cl.qfc1110_kn1  * P11 / FQC;


value, IQDB1005;
value, IQFC1010;
value, IQDB1015;
value, IQDB1105;
value, IQFC1110;

write, table=twiss, file=SECS10-11.twiss;


atmid11: beta0,  betx=table(twiss,MID11,beta11),
                 alfx=table(twiss,MID11,alfa11),
                 bety=table(twiss,MID11,beta22),
                 alfy=table(twiss,MID11,alfa22),
                 mux =table(twiss,MID11,mu1),
                 muy =table(twiss,MID11,mu2),
                 x   =table(twiss,MID11,x),
                 px  =table(twiss,MID11,px),
                 y   =table(twiss,MID11,y),
                 py  =table(twiss,MID11,py),
                 t   = 0,
                 pt  = 0,
!                 t   =table(twiss,MID11,t),
!                 pt  =table(twiss,MID11,pt),
                 dx  =table(twiss,MID11,disp1),
                 dpx =table(twiss,MID11,disp2),
                 dy  =table(twiss,MID11,disp3),
                 dpy =table(twiss,MID11,disp4);

betzatmid11=table(twiss,MID11,beta33);
alfzatmid11=table(twiss,MID11,alfa33);

eneatmid11 = table(twiss,MID11,energy);

show, atmid11;
value, eneatmid11;




print, text="###################################################################################";
print, text="###################################################################################";
print, text="###################################################################################";
print, text="Track Twiss parameters from MID11 till the end of girder 11, get init12 this way";


BEAM, PARTICLE=ELECTRON, ENERGY=eneatmid11, sigt=1, sige=1/1000 ;

use, period=SEC11fm;!= ( D1160_2, CL.ACS1160, D1190, CL.BPM1190 );

ptc_create_universe;
ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
ptc_setswitch, debuglevel=0, nocavity=false, 
               maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
PTC_TWISS, table=twiss, icase=6, no=2, BETA0=atmid11, betz=betzatmid11, alfz=alfzatmid11;
ptc_end;


write, table=twiss, file=SEC11fm.twiss;

init12: beta0,  betx=table(twiss,SEC11fm$END,beta11),
                alfx=table(twiss,SEC11fm$END,alfa11),
                bety=table(twiss,SEC11fm$END,beta22),
                alfy=table(twiss,SEC11fm$END,alfa22),
                mux =table(twiss,SEC11fm$END,mu1),
                muy =table(twiss,SEC11fm$END,mu2),
                x   =table(twiss,SEC11fm$END,x),
                px  =table(twiss,SEC11fm$END,px),
                y   =table(twiss,SEC11fm$END,y),
                py  =table(twiss,SEC11fm$END,py),
                t   = 0,
                pt  = 0,
!                t   =table(twiss,SEC11fm$END,t),
!                pt  =table(twiss,SEC11fm$END,pt),
                dx  =table(twiss,SEC11fm$END,disp1),
                dpx =table(twiss,SEC11fm$END,disp2),
                dy  =table(twiss,SEC11fm$END,disp3),
                dpy =table(twiss,SEC11fm$END,disp4);


eneinit12 = table(twiss,SEC11fm$END,energy);

show, init12;
value, eneinit12;

print, text="###################################################################################";
print, text="###################################################################################";
print, text="###################################################################################";
print, text="Match quads at girders 12 to periodic solution at MID12, starting with init12";


BEAM, PARTICLE=ELECTRON, ENERGY=eneinit12 , sigt=1, sige=1/1000 ;
use, period=SEC12tm;

if (skipprematch == 0)
 {

   match, use_macro; 
      vary,name=IQDC1205,step=0.01, lower=  0   , upper=100;!Line SEC12, element 2,  currnt of the quad
      vary,name=IQFD1210,step=0.01, lower=  0   , upper=200;!Line SEC12, element 4,  currnt of the quad

       m12: macro =
        {
          TWISS, beta0=init12;

          value, IQDC1205;
          value, IQFD1210;

          value, table(twiss,MID12,betx) - MIDCELL->BETX;
          value, table(twiss,MID12,alfx) - MIDCELL->ALFX;
          value, table(twiss,MID12,bety) - MIDCELL->BETY;
          value, table(twiss,MID12,alfy) - MIDCELL->ALFY;
        };

      CONSTRAINT, expr = table(twiss,MID12,betx) = MIDCELL->BETX;
      CONSTRAINT, expr = table(twiss,MID12,alfx) = MIDCELL->ALFX;
      CONSTRAINT, expr = table(twiss,MID12,bety) = MIDCELL->BETY;
      CONSTRAINT, expr = table(twiss,MID12,alfy) = MIDCELL->ALFY;

      jacobian,calls=100, tolerance=1e-04;
   !   lmdif,calls=1000, tolerance=1e-12;
   !   migrad,calls=1000, tolerance=1e-12;
   endmatch;
 }

match, use_ptcknob;
!12
   ptc_varyknob, element=CL.QDC12, kn=1, exactmatch=false, !Line SEC12 , element 4,  currnt of the quad
                            lower= -100 *FQC / P12, 
                            upper=  0,
	        trustrange = 0.2, step=1e-6 ;
   ptc_varyknob, element=CL.QFD1210, kn=1, exactmatch=true, !Line SEC12 , 
                            lower=  0, 
                            upper=  200 *FQD / P12, 
	        trustrange = 0.2, step=1e-6 ;

   ptc_create_universe;
   ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
   ptc_setswitch, debuglevel=0, nocavity=false, 
                  maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
   PTC_TWISS, table=twiss, icase=6, no=3, 
              beta0=init12, betz=1.0;


   CONSTRAINT, weight=10,  expr = table(twiss,MID12,beta11) = MIDCELL->BETX;
   CONSTRAINT, weight=100, expr = table(twiss,MID12,alfa11) = MIDCELL->ALFX;
   CONSTRAINT, weight=10,  expr = table(twiss,MID12,beta22) = MIDCELL->BETY;
   CONSTRAINT, weight=100, expr = table(twiss,MID12,alfa22) = MIDCELL->ALFY;

!   jacobian,    calls=100, tolerance=1e-06;
   lmdif,    calls=1000, tolerance=1e-09;
endmatch;


IQDC1205 = -cl.qdc12_kn1    * P12 / FQC;
IQFD1210 =  cl.qfd1210_kn1  * P12 / FQD;

value, IQDC1205;
value, IQFD1210;


write, table=twiss, file=SEC12.twiss;


!remember beta0 at MID12
atMID12: beta0,  betx=table(twiss,MID12,beta11),
                 alfx=table(twiss,MID12,alfa11),
                 bety=table(twiss,MID12,beta22),
                 alfy=table(twiss,MID12,alfa22),
                 mux =table(twiss,MID12,mu1),
                 muy =table(twiss,MID12,mu2),
                 x   =table(twiss,MID12,x),
                 px  =table(twiss,MID12,px),
                 y   =table(twiss,MID12,y),
                 py  =table(twiss,MID12,py),
                 t   = 0,
                 pt  = 0,
!                 t   =table(twiss,MID12,t),
!                 pt  =table(twiss,MID12,pt),
                 dx  =table(twiss,MID12,disp1),
                 dpx =table(twiss,MID12,disp2),
                 dy  =table(twiss,MID12,disp3),
                 dpy =table(twiss,MID12,disp4);

eneatmid12 = table(twiss,MID12,energy);

show, atMID12;
value, eneatmid12;

print, text="###################################################################################";
print, text="###################################################################################";
print, text="###################################################################################";
print, text="Track Twiss parameters from MID11 till the end of girder 11, get init12 this way";


!get twiss parameters at the end of girder 12, starting from MID12


BEAM, PARTICLE=ELECTRON, ENERGY=eneatmid12, sigt=1, sige=1/1000 ;
use, period=SEC12fm;!= ( D1160_2, CL.ACS1160, D1190, CL.BPM1190 );

ptc_create_universe;
ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
ptc_setswitch, debuglevel=0, nocavity=false, 
               maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
PTC_TWISS, table=twiss, icase=6, no=2, BETA0=atMID12, betz=1.0;
ptc_end;

write, table=twiss, file=SEC12fm.twiss;

!parameters at the end od SEC12 are initial for sec13
init13: beta0,  betx=table(twiss,SEC12fm$END,beta11),
                alfx=table(twiss,SEC12fm$END,alfa11),
                bety=table(twiss,SEC12fm$END,beta22),
                alfy=table(twiss,SEC12fm$END,alfa22),
                mux =table(twiss,SEC12fm$END,mu1),
                muy =table(twiss,SEC12fm$END,mu2),
                x   =table(twiss,SEC12fm$END,x),
                px  =table(twiss,SEC12fm$END,px),
                y   =table(twiss,SEC12fm$END,y),
                py  =table(twiss,SEC12fm$END,py),
                t   = 0,
                pt  = 0,
!                t   =table(twiss,SEC12fm$END,t),
!                pt  =table(twiss,SEC12fm$END,pt),
                dx  =table(twiss,SEC12fm$END,disp1),
                dpx =table(twiss,SEC12fm$END,disp2),
                dy  =table(twiss,SEC12fm$END,disp3),
                dpy =table(twiss,SEC12fm$END,disp4);


eneinit13 = table(twiss,SEC12fm$END,energy);

show, init13;
value, eneinit13;

print, text="###################################################################################";
print, text="###################################################################################";
print, text="###################################################################################";
print, text="Match quads at girders 13 to periodic solution at MID13, starting with init13";


BEAM, PARTICLE=ELECTRON, ENERGY=eneinit13, sigt=1, sige=1/1000 ;
use, period=SEC13tm;

if (skipprematch == 0)
 {

    match, use_macro; 
       vary,name=IQDD1305,step=0.01, lower=  0   , upper=100;
       vary,name=IQFD1310,step=0.01, lower=  0   , upper=200;

        m13a: macro =
         {
           TWISS, beta0=init13;

           value, IQDD1305, IQFD1310;
           value, table(twiss,MID13,betx) - MIDCELL->BETX,
                  table(twiss,MID13,alfx) - MIDCELL->ALFX,
                  table(twiss,MID13,bety) - MIDCELL->BETY,
                  table(twiss,MID13,alfy) - MIDCELL->ALFY;

           value, table(twiss,MID13,betx) - MIDCELL->BETX,
                  table(twiss,MID13,alfx) - MIDCELL->ALFX,
                  table(twiss,MID13,bety) - MIDCELL->BETY,
                  table(twiss,MID13,alfy) - MIDCELL->ALFY;
         };

       CONSTRAINT, expr = table(twiss,MID13,betx) = MIDCELL->BETX;
       CONSTRAINT, weight=10, expr = table(twiss,MID13,alfx) = MIDCELL->ALFX;
       CONSTRAINT, expr = table(twiss,MID13,bety) = MIDCELL->BETY;
       CONSTRAINT, weight=10, expr = table(twiss,MID13,alfy) = MIDCELL->ALFY;

       lmdif, calls=100, tolerance=1e-06;
    !   migrad,   calls=300, tolerance=1e-04;
    endmatch;

 }



match, use_ptcknob;

!13
   ptc_varyknob, element=CL.QDD13, kn=1, exactmatch=false, 
                            lower= -100 *FQD / P13, 
                            upper=  0, 
	        trustrange = 0.2, step=1e-3 ;
   ptc_varyknob, element=CL.QFD1310, kn=1, exactmatch=true, 
                            lower=  0, 
                            upper=  200 *FQD / P13, 
	        trustrange = 0.2, step=1e-3 ;

   ptc_create_universe;
   ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
   ptc_setswitch, debuglevel=0, nocavity=false, 
                  maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
   PTC_TWISS, table=twiss, icase=6, no=3, beta0=init13, betz=1.0;


   CONSTRAINT, weight=10,  expr = table(twiss,MID13,beta11) = MIDCELL->BETX;
   CONSTRAINT, weight=100, expr = table(twiss,MID13,alfa11) = MIDCELL->ALFX;
   CONSTRAINT, weight=10,  expr = table(twiss,MID13,beta22) = MIDCELL->BETY;
   CONSTRAINT, weight=100, expr = table(twiss,MID13,alfa22) = MIDCELL->ALFY;

   jacobian,    calls=100, tolerance=1e-08;
!   migrad,    calls=100, tolerance=1e-04;
endmatch;


IQDD1305 = -cl.qdd13_kn1    * P13 / FQD;
IQFD1310 =  cl.qfd1310_kn1  * P13 / FQD;


value, IQDD1305;
value, IQFD1310;

write, table=twiss, file=SEC13.twiss;


!remember beta0 at MID12
atMID13: beta0,  betx=table(twiss,MID13,beta11),
                 alfx=table(twiss,MID13,alfa11),
                 bety=table(twiss,MID13,beta22),
                 alfy=table(twiss,MID13,alfa22),
                 mux =table(twiss,MID13,mu1),
                 muy =table(twiss,MID13,mu2),
                 x   =table(twiss,MID13,x),
                 px  =table(twiss,MID13,px),
                 y   =table(twiss,MID13,y),
                 py  =table(twiss,MID13,py),
                 t   = 0,
                 pt  = 0,
 !                t   =table(twiss,MID13,t),
!                 pt  =table(twiss,MID13,pt),
                 dx  =table(twiss,MID13,disp1),
                 dpx =table(twiss,MID13,disp2),
                 dy  =table(twiss,MID13,disp3),
                 dpy =table(twiss,MID13,disp4);

eneatmid13 = table(twiss,MID13,energy);

show, atMID13;
value, eneatmid13;

print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"
print, text="Track Twiss parameters from MID13 till the end of girder 13, get init14 this way";

!get twiss parameters at the end of girder 12, starting from MID12

BEAM, PARTICLE=ELECTRON, ENERGY=eneatmid13 , sigt=1, sige=1/1000 ;
use, period=SEC13fm;!= ( D1160_2, CL.ACS1160, D1190, CL.BPM1190 );

ptc_create_universe;
ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
ptc_setswitch, debuglevel=0, nocavity=false, 
               maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
PTC_TWISS, table=twiss, icase=6, no=2, BETA0=atMID13, betz=1.0;
ptc_end;


write, table=twiss, file=SEC13fm.twiss;

!parameters at the end od SEC13 are initial for sec14
init14: beta0,  betx=table(twiss,SEC13fm$END,beta11),
                alfx=table(twiss,SEC13fm$END,alfa11),
                bety=table(twiss,SEC13fm$END,beta22),
                alfy=table(twiss,SEC13fm$END,alfa22),
                mux =table(twiss,SEC13fm$END,mu1),
                muy =table(twiss,SEC13fm$END,mu2),
                x   =table(twiss,SEC13fm$END,x),
                px  =table(twiss,SEC13fm$END,px),
                y   =table(twiss,SEC13fm$END,y),
                py  =table(twiss,SEC13fm$END,py),
                t   = 0,
                pt  = 0,
!                t   =table(twiss,SEC13fm$END,t),
!                pt  =table(twiss,SEC13fm$END,pt),
                dx  =table(twiss,SEC13fm$END,disp1),
                dpx =table(twiss,SEC13fm$END,disp2),
                dy  =table(twiss,SEC13fm$END,disp3),
                dpy =table(twiss,SEC13fm$END,disp4);


eneinit14 = table(twiss,SEC13fm$END,energy);

show, init14;
value, eneinit14;

print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"
print, text="Match quads at girders 14 to periodic solution at MID14, starting with init14";


BEAM, PARTICLE=ELECTRON, ENERGY=eneinit14, sigt=1, sige=1/1000 ;
use, period=SEC14tm;

if (skipprematch == 0)
 {
   match, use_macro; 
      vary,name=IQDD1405,step=0.01, lower=  0   , upper=100;
      vary,name=IQFD1410,step=0.01, lower=  0   , upper=200;

       m14a: macro =
        {
          TWISS, beta0=init14;

          value, IQDD1405;
          value, IQFD1410;

          value, table(twiss,MID14,betx) - MIDCELL->BETX;
          value, table(twiss,MID14,alfx) - MIDCELL->ALFX;
          value, table(twiss,MID14,bety) - MIDCELL->BETY;
          value, table(twiss,MID14,alfy) - MIDCELL->ALFY;
        };

      CONSTRAINT, expr = table(twiss,MID14,betx) = MIDCELL->BETX;
      CONSTRAINT, expr = table(twiss,MID14,alfx) = MIDCELL->ALFX;
      CONSTRAINT, expr = table(twiss,MID14,bety) = MIDCELL->BETY;
      CONSTRAINT, expr = table(twiss,MID14,alfy) = MIDCELL->ALFY;

      lmdif, calls=100, tolerance=1e-03;
   !   migrad,   calls=1000,tolerance=1e-04;
   endmatch;

 }

match, use_ptcknob;

   ptc_varyknob, element=CL.QDD14, kn=1, exactmatch=false, 
                            lower= -100 *FQD / P14, 
                            upper=  0, 
	        trustrange = 0.5 , step=1e-3;
   ptc_varyknob, element=CL.QFD1410, kn=1, exactmatch=true, 
                            lower=  0, 
                            upper=  200 *FQD / P14, 
	        trustrange = 0.5, step=1e-3 ;


   ptc_create_universe;
   ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
   ptc_setswitch, debuglevel=0, nocavity=false, 
                  maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
   PTC_TWISS, table=twiss, icase=6, no=3, beta0=init14, betz=1.0;


   CONSTRAINT, weight=10,  expr = table(twiss,MID14,beta11) = MIDCELL->BETX;
   CONSTRAINT, weight=100, expr = table(twiss,MID14,alfa11) = MIDCELL->ALFX;
   CONSTRAINT, weight=10,  expr = table(twiss,MID14,beta22) = MIDCELL->BETY;
   CONSTRAINT, weight=100, expr = table(twiss,MID14,alfa22) = MIDCELL->ALFY;

!   jacobian,    calls=100, tolerance=1e-04;
   lmdif,    calls=100, tolerance=1e-04;
endmatch;

IQDD1405 = -cl.qdd14_kn1    * P14 / FQD;
IQFD1410 =  cl.qfd1410_kn1  * P14 / FQD;

value, IQDD1405;
value, IQFD1410;

write, table=twiss, file=SEC14.twiss;


!remember beta0 at MID12
atMID14: beta0,  betx=table(twiss,MID14,beta11),
                 alfx=table(twiss,MID14,alfa11),
                 bety=table(twiss,MID14,beta22),
                 alfy=table(twiss,MID14,alfa22),
                 mux =table(twiss,MID14,mu1),
                 muy =table(twiss,MID14,mu2),
                 x   =table(twiss,MID14,x),
                 px  =table(twiss,MID14,px),
                 y   =table(twiss,MID14,y),
                 py  =table(twiss,MID14,py),
                 t   = 0,
                 pt  = 0,
!                 t   =table(twiss,MID14,t),
!                 pt  =table(twiss,MID14,pt),
                 dx  =table(twiss,MID14,disp1),
                 dpx =table(twiss,MID14,disp2),
                 dy  =table(twiss,MID14,disp3),
                 dpy =table(twiss,MID14,disp4);

eneatmid14 = table(twiss,MID14,energy);

show, atMID14;
value, eneatmid14;

print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"
print, text="Track Twiss parameters from MID14 till the end of girder 14, get init14 this way";

!get twiss parameters at the end of girder 12, starting from MID12

BEAM, PARTICLE=ELECTRON, ENERGY=eneatmid14 , sigt=1, sige=1/1000 ;
use, period=SEC14fm;!= ( D1160_2, CL.ACS1160, D1190, CL.BPM1190 );

ptc_create_universe;
ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
ptc_setswitch, debuglevel=0, nocavity=false, 
               maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
PTC_TWISS, table=twiss, icase=4, no=2, BETA0=atMID14, betz=1.0;
ptc_end;


!parameters at the end od SEC14 are initial for sec14
init15: beta0,  betx=table(twiss,SEC14fm$END,beta11),
                alfx=table(twiss,SEC14fm$END,alfa11),
                bety=table(twiss,SEC14fm$END,beta22),
                alfy=table(twiss,SEC14fm$END,alfa22),
                mux =table(twiss,SEC14fm$END,mu1),
                muy =table(twiss,SEC14fm$END,mu2),
                x   =table(twiss,SEC14fm$END,x),
                px  =table(twiss,SEC14fm$END,px),
                y   =table(twiss,SEC14fm$END,y),
                py  =table(twiss,SEC14fm$END,py),
                t   = 0,
                pt  = 0,
!                t   =table(twiss,SEC14fm$END,t),
!                pt  =table(twiss,SEC14fm$END,pt),
                dx  =table(twiss,SEC14fm$END,disp1),
                dpx =table(twiss,SEC14fm$END,disp2),
                dy  =table(twiss,SEC14fm$END,disp3),
                dpy =table(twiss,SEC14fm$END,disp4);


eneinit15 = table(twiss,SEC14fm$END,energy);

show, init15;
value, eneinit15;


print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"
print, text="Match quads at girders 15 to periodic solution at MID14, starting with init14";


BEAM, PARTICLE=ELECTRON, ENERGY=eneinit15 , sigt=1, sige=1/1000 ;
use, period=SEC15tm;

match, use_macro; 
   vary,name=IQDD1505,step=0.01, lower=  0   , upper=100;
   vary,name=IQFD1510,step=0.01, lower=  0   , upper=200;

    m15a: macro =
     {
       TWISS, beta0=init15;

       value, IQDD1505;
       value, IQFD1510;
       
       value, table(twiss,MID15,betx) - 4;
       value, table(twiss,MID15,alfx) - 0;
       value, table(twiss,MID15,bety) - 4;
       value, table(twiss,MID15,alfy) - 0;
     };

   CONSTRAINT, expr = table(twiss,MID15,betx) = 4.0;
   CONSTRAINT, expr = table(twiss,MID15,alfx) = 0;
   CONSTRAINT, expr = table(twiss,MID15,bety) = 4.0;
   CONSTRAINT, expr = table(twiss,MID15,alfy) = 0;

   lmdif, calls=100, tolerance=1e-03;
!   migrad,   calls=1000,tolerance=1e-04;
endmatch;


match, use_ptcknob;

   ptc_varyknob, element=CL.QDD15, kn=1, exactmatch=false, 
                            lower= -100 *FQD / P15, 
                            upper=  0, 
	        trustrange = 0.5, step=1e-3 ;


   ptc_varyknob, element=CL.QFD1510, kn=1, exactmatch=true, 
                            lower=  0, 
                            upper=  200 *FQD / P15, 
	        trustrange = 0.5, step=1e-3 ;


   ptc_create_universe;
   ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
   ptc_setswitch, debuglevel=0, nocavity=false, 
                  maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
   PTC_TWISS, table=twiss, icase=6, no=3,beta0=init15, betz=1.0;

   CONSTRAINT, weight=10,  expr = table(twiss,MID15,beta11) = 4.;
   CONSTRAINT, weight=100, expr = table(twiss,MID15,alfa11) = 0.;
   CONSTRAINT, weight=10,  expr = table(twiss,MID15,beta22) = 4.;
   CONSTRAINT, weight=100, expr = table(twiss,MID15,alfa22) = 0.;

   lmdif, calls=500,tolerance=1e-06;
!   migrad,   calls=200,tolerance=1e-04;
endmatch;

IQDD1505 = -cl.qdd15_kn1    * P15 / FQD;
IQFD1510 =  cl.qfd1510_kn1  * P15 / FQD;

value, IQDD1505;
value, IQFD1510;

write, table=twiss, file=SEC15.twiss;

!remember beta0 at MID12
atMID15: beta0,  betx=table(twiss,MID15,beta11),
                 alfx=table(twiss,MID15,alfa11),
                 bety=table(twiss,MID15,beta22),
                 alfy=table(twiss,MID15,alfa22),
                 mux =table(twiss,MID15,mu1),
                 muy =table(twiss,MID15,mu2),
                 x   =table(twiss,MID15,x),
                 px  =table(twiss,MID15,px),
                 y   =table(twiss,MID15,y),
                 py  =table(twiss,MID15,py),
                 t   = 0,
                 pt  = 0,
!                 t   =table(twiss,MID15,t),
!                 pt  =table(twiss,MID15,pt),
                 dx  =table(twiss,MID15,disp1),
                 dpx =table(twiss,MID15,disp2),
                 dy  =table(twiss,MID15,disp3),
                 dpy =table(twiss,MID15,disp4);

eneatmid15 = table(twiss,MID15,energy);

show, atMID15;
value, eneatmid15;




print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"
print, text="###################################################################################"

! recalculate and plot optics

value, E8;
BEAM, PARTICLE=ELECTRON, ENERGY=E8/1000. , sigt=1, sige=1/1000 ;
USE, period=LINAC10; != ( SEC10Q,  SEC11, SEC12, SEC13, SEC14, SEC15 )

!ptc_setswitch, debuglevel=2;

ptc_create_universe;
ptc_create_layout, model=1, method=6, nst=PTCNSTEPS, exact=true, closed_layout=false;
ptc_setswitch, debuglevel=0, nocavity=false, 
               maxacceleration=true, fringe=true, exact_mis=true, time=true, totalpath=true;
PTC_TWISS, table=twiss, icase=6, no=2, BETA0=B0_10, betz=1.0;
ptc_end;



print, text="PLOT, noversion=true, haxis=s, vaxis=betx, bety";
PLOT, noversion=true, haxis=s, vaxis=beta11, beta22, colour=100, spline, title="Rematched optics (PTC)";

PLOT, noversion=true, haxis=s, vaxis=beta33, colour=100, spline, title="Rematched optics (PTC)";

print, text="PLOT, noversion=true, haxis=s, vaxis=sigx, sigy";
PLOT, noversion=true, haxis=s, vaxis=ptcsigx, ptcsigy, colour=100, spline, title="Rematched optics (PTC)";

print, text="PLOT, noversion=true, haxis=s, vaxis=sigp";
PLOT, noversion=true, haxis=s, vaxis=ptcsigp, colour=100, spline, title="Rematched optics (PTC)";

print, text="PLOT, noversion=true, haxis=s, vaxis=energy";
PLOT, noversion=true, haxis=s, vaxis=energy, colour=100, title="Rematched optics (PTC)";

write, table=twiss, file="LINAC10.rematched.twiss";

select,flag=save, clear;
select,flag=save, pattern="iq.*";
save, file="linac-10-15.madout";
save, file="TUNES";



value, MIDCELL->BETX, MIDCELL->ALFX, MIDCELL->BETY, MIDCELL->ALFY;
value, ENERGY=E8/1000., eneatmid11, eneinit12, eneatmid12, eneinit13, eneatmid13, 
             eneinit14, eneatmid14, eneinit15, eneatmid15;

value, atmid11->betx, atmid12->betx,  atmid13->betx,  atmid14->betx, atmid15->betx;
value, table(twiss,MID11,beta11),
       table(twiss,MID12,beta11),
       table(twiss,MID13,beta11),
       table(twiss,MID14,beta11),
       table(twiss,MID15,beta11);

value, atmid11->bety, atmid12->bety,  atmid13->bety,  atmid14->bety, atmid15->bety;
value, table(twiss,MID11,beta22),
       table(twiss,MID12,beta22),
       table(twiss,MID13,beta22),
       table(twiss,MID14,beta22),
       table(twiss,MID15,beta22);

value, atmid11->alfx, atmid12->alfx,  atmid13->alfx,  atmid14->alfx, atmid15->alfx;
value, table(twiss,MID11,alfa11),
       table(twiss,MID12,alfa11),
       table(twiss,MID13,alfa11),
       table(twiss,MID14,alfa11),
       table(twiss,MID15,alfa11);

value, atmid11->alfy, atmid12->alfy,  atmid13->alfy,  atmid14->alfy, atmid15->alfy;
value, table(twiss,MID11,alfa22),
       table(twiss,MID12,alfa22),
       table(twiss,MID13,alfa22),
       table(twiss,MID14,alfa22),
       table(twiss,MID15,alfa22);

value, init12->betx, init13->betx, init14->betx, init15->betx;
value, init12->bety, init13->bety, init14->bety, init15->bety;
value, init12->alfx, init13->alfx, init14->alfx, init15->alfx;
value, init12->alfy, init13->alfy, init14->alfy, init15->alfy;



system,"cp madx.ps linac-10-15_gbg_par.ps";


print, text="Here we STOP";
stop;
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  SSSSS  TTTTT   OOO   PPP
!  S        T    O   O  P  P
!  SSSSS    T    O   O  PPP
!      S    T    O   O  P
!  SSSSS    T     OOO   P
!