! 2011/11/22 R. De Maria Generate optimized triplet layouts

Option, -warn,-info,-echo;

lq1=0.50;
lq2=0.50;
endseq=165;
ffocus: sequence,refer = CENTRE, L := endseq;
q11 :quadrupole,l:=l1/2,k1:= k1,at:=lstar+l1/4;
q12 :quadrupole,l:=l1/2,k1:= k1,at:=lstar+l1*3/4+lq1;
q2a1:quadrupole,l:=l2/2,k1:=-k2,at:=lstar+l1+lq1+d12+l2/4;
q2a2:quadrupole,l:=l2/2,k1:=-k2,at:=lstar+l1+lq1+lq2+d12+l2*3/4;
q2b1:quadrupole,l:=l2/2,k1:=-k2,at:=lstar+l1+lq1+lq2+d12+d22+l2*5/4;
q2b2:quadrupole,l:=l2/2,k1:=-k2,at:=lstar+l1+lq1+2*lq2+d12+d22+l2*7/4;
q31 :quadrupole,l:=l3/2,k1:= k3,at:=lstar+l1+lq1+2*lq2+d12+2*l2+d22+d23+l3/4;
q32 :quadrupole,l:=l3/2,k1:= k3,at:=lstar+l1+2*lq1+2*lq2+d12+2*l2+d22+d23+l3*3/4;
endsequence;
beam;

r_tol_IT=0.0006;
h_tol_IT=0.0010;
v_tol_IT=0.0010;

calc_triplet(bxs,bys): macro= {
   use,sequence=ffocus;
   twiss,betx=bxs,bety=bys;
   bq2a1=table(twiss,q2a1,betymax);
   bq2a2=table(twiss,q2a2,betymax);
   bq2b1=table(twiss,q2b1,betymax);
   bq2b2=table(twiss,q2b2,betymax);
   bq31=table(twiss,q31,betxmax);
   bq32=table(twiss,q32,betxmax);
   bq2max=bq2a1;
   if (bq2a2>bq2max){ bq2max=bq2a2;};
   if (bq2b1>bq2max){ bq2max=bq2b1;};
   if (bq2b2>bq2max){ bq2max=bq2b2;};
   bq3max=bq32;
   alfxend=table(twiss,FFOCUS$END,alfx);
   alfyend=table(twiss,FFOCUS$END,alfy);
   betxend=table(twiss,FFOCUS$END,betx);
   betyend=table(twiss,FFOCUS$END,bety);
};

mk_triplet(bxs,bys,grad,posq4,betxq4,betyq4): macro={
imbmax=1;
endseq=posq4;

tbetx:=table(twiss,betx); tbety:=table(twiss,bety);
talfx:=table(twiss,alfx); talfy:=table(twiss,alfy);
tk1l :=table(twiss,k1l); tl:=table(twiss,l);
betxmax:=tbetx+talfx^2/tbetx/abs(tk1l/tl);
betymax:=tbety+talfy^2/tbety/abs(tk1l/tl);
select,flag=twiss,clear;
select,flag=twiss,column=name,s,l,angle,k1l,hkick,vkick,kick,betx,bety,alfx,alfy,dx,dpx,dy,dpy,mux,muy,x,y,px,py,betxmax,betymax;

brho=7e12/clight; k=grad/brho;

lstar=23;

l1=7.145; l2=6.735;
if (grad<120 )  { l1=9.145; l2=7.735; };
d12:=d23; l3:=l1;

k1:=k; k2:=k; k3:=k;
match,use_macro;
vary,name=l1;
vary,name=l2;
vary,name=k3;
use_macro,name=calc_triplet(bxs,bys);
constraint, expr= bq3max = imbmax*bq2max;
constraint, expr= betxend  =betxq4;
constraint, expr= betyend  =betyq4;
constraint, expr= 1e5*alfxend    >2e5;
constraint, expr= 1e5*alfyend    >2e5;
jacobian,calls=60,tolerance=1e-14,bisec=9;
endmatch;

if (k3>k) {
k1:=k; k2=k; k3=k;
match,use_macro;
vary,name=l1;
vary,name=l2;
vary,name=k1;
use_macro,name=calc_triplet(bxs,bys);
constraint, expr= bq3max = imbmax*bq2max;
constraint, expr= betxend  =betxq4;
constraint, expr= betyend  =betyq4;
constraint, expr= 1e5*alfxend    >2e5;
constraint, expr= 1e5*alfyend    >2e5;
jacobian,calls=60,tolerance=1e-14,bisec=9;
endmatch;
};

if (tar>1e-10){stop;};

bq3orig=bq3max;

exec,round(l1,1000);
exec,round(l2,1000);
exec,round(l3,1000);
exec,round(d12,1000);
exec,round(d22,1000);
exec,round(d23,1000);

match,use_macro;
vary,name=k1;
vary,name=k2;
vary,name=k3;
use_macro,name=calc_triplet(bxs,bys);
constraint, expr= bq3max = imbmax*bq2max;
constraint, expr= betxend  =betxq4;
constraint, expr= betyend  =betyq4;
constraint, expr= 1e5*alfxend    >2e5;
constraint, expr= 1e5*alfyend    >2e5;
jacobian,calls=60,tolerance=1e-14,bisec=9;
endmatch;

l.mqxl:=l1;l.mqx:=l2;
dq1q2a:=d12;dq2aq2b:=d22;dq2bq3:=d23;dq1aq1b=d11;
kqx1.l1 =-k1;  kqx1.r1=k1;  kqx1.l5=-k1;  kqx1.r5=k1;
kqx2a.l1=-k2; kqx2a.r1=k2; kqx2a.l5=-k2; kqx2a.r5=k2;
kqx2b.l1=-k2; kqx2b.r1=k2; kqx2b.l5=-k2; kqx2b.r5=k2;
kqx3.l1 =-k3;  kqx3.r1=k3;  kqx3.l5=-k3;  kqx3.r5=k3;
exec,calc_triplet(bxs,bys);
g1=kqx1.r1*brho;g2=kqx2a.r1*brho; g3=kqx3.r5*brho;
value,g1,g2,g3;
totlen:=lstar+l1+d12+2*l2+d22+d23+l3+4*lqq;
show,l.mqxl,l.mqx,totlen;
value,l.mqxl,l.mqx,dq1q2a,dq2aq2b,dq2bq3,lqq,totlen;
value,bq3orig,bq3max;
value,bq2a1,bq2a2,bq2b1,bq2b2,bq31,bq32;
value,alfxend;
value,alfyend;
value,bq3max*bxs;
};

__round(xv,xn): macro={ if ( xv-xn > 0.5) { xv=xn+1; } else { xv=xn; }; };
round(xv,digits): macro={
  xtmp=xv*digits;
  exec,__round(xtmp,$xtmp);
  xv=xtmp/digits;
};

exec,mk_triplet(0.60,0.60,150,120,1913,1234);