(* This program computes the multiturn injection in an acceptance(aksept)
 in pi mmmrad.The beam is injected at sent(mm),the bump starts at
 sbump(mm) and decreases by ds(mm) every turn. The injection septum has 
 a thickness of eps(mm) and is placed at septpos mm(outside position). The 
 horizontal dimension of the incoming beam is x(mm)(given by the emitance emitt[pimmmrad]
  considered at 2 sigmas,
  it's why the distribution considered gaussian is with a 2 instead of 0.5). 
  The tune is q. The number of turn we try to inject is jmax given by the beam
   dimension and the acceptance. The mismatch is defined by ex. The number of turn 
   observed is jt and given by the decreases of the bump observed(btot). 
   beta is the horizontal Twiss parameter
 at the injection point. emitt(pimmmrad) is the emittance( 2 sigma) of the incoming 
beam, emmax(pimmmrad) is the acceptance of the machine The drawings are 
shown in the normalised phase space.

	M Chanel editor last modif. 21/08/2001.*)


sent=45;sbump=40;ds:=-1.3;eps=1;
septpos=40;
btot=40;beta=1.6;emitt=20;emmax=180;
ex=Sqrt[2];aksept=Sqrt[emmax*beta];x=Sqrt[emitt*beta];nx=10;dx=x/nx;
q=2*Pi*2.20//N;septum=septpos-eps;
jt=Round[Abs[btot/ds]]+3;sent=sent-x*(ex-1)//N;
er=sent-sbump;jmax=Floor[(sent-aksept-sbump)/ds]+1;

u=2Pi//N;nu=30;du=u/(nu-1);
m=N[{{Cos[q],Sin[q]},{-Sin[q],Cos[q]}}];

Print["jmax=",jmax," jt=",jt]
g0=Graphics[Line[{{septum,-40},{septpos,-40},{septpos,40},
		{septum,40},{septum,-40}}]];
gc[p1_]:=Graphics[Disk[p1,1]];
ga[p1_]:=Graphics[Circle[p1,aksept]]

dist[x_,y_]:=Sqrt[(y[[1]]-x[[1]])^2+(y[[2]]-x[[2]])^2];
coord[w_,b_,bu_]:=(1-w)*bu+w*b;

(*bump à chaque tour*)
bu=Table[{Max[sbump+(j-1)*ds,sbump-btot],0},{j,1,jt}];
(*une couleur*)
k1[a_,b_,c_,r_]:=Module[{z},z=1-Exp[-2*(r/x)^2];
	RGBColor[a+(1-a)*z,b+(1-b)*z,c+(1-c)*z]]
(*une ellipse*)
l[a_,b_,or_]:=Module[{r1,r2},r1=Cos[or];r2=Sin[or];
	Table[{a*r1*Sin[t]+b*r2*Cos[t],
	a*r2*Sin[t]+b*r1*Cos[t]}//N,{t,0,u,du}]];

Print["le faisceau injecté"]
zi=Table[l[y/ex,y*ex,0],{y,x,dx,-dx}];ztot=zi;
Do[If[zi[[i,k,1]]<septpos-sent,zi[[i,k,1]]=septpos-sent],{i,1,nx},{k,1,nu}]
bi=Table[{0,0},{p,1,jt},{j,1,jt},{i,1,nx},{k,1,nu}];
Do[bi[[p,p,i,k]]=zi[[i,k]]+{sent,0},{p,1,jmax},
				{i,1,nx},{k,1,nu}];
Print["test aksept"]
Do[Module[{z,w},z=dist[bi[[p,p,i,k]],bu[[p]]];w=aksept/z;
	If[w<1,bi[[p,p,i,k]]=coord[w,bi[[p,p,i,k]],bu[[p]]] ]],
	{p,1,jmax},{i,1,nx},{k,1,nu}]

Print["Calcul les tours"]
Do[bi[[p,j,i,k]]=m.(bi[[p,j-1,i,k]]-bu[[j-1]])+bu[[j]];
	If[bi[[p,j,i,k,1]]>septum,bi[[p,j,i,k,1]]=septum],
		{p,1,jmax},{j,p+1,jt},{i,1,nx},{k,1,nu}];

Print["calcul graphics"]
bg=Table[Graphics[{k1[1-p/jmax,p/jmax,0,(nx+1-i)*dx],
	Polygon[bi[[p,j,i]]]}],{p,1,jmax},
			{j,1,jt},{i,1,nx}];

Print["draw it"]
graphdim={{-20,60},{-20,20}};
graphratio=(graphdim[[2,2]]-graphdim[[2,1]])/(graphdim[[1,2]]-graphdim[[1,1]]);
(*Do[Show[{g0,bg[[1,j]],gc[bu[[j]]],ga[bu[[j]]]},
	Frame->True,PlotRange->{{-20,60},{-30,30}},
	AspectRatio->graphratio],{j,1,jt}]
Do[Show[{g0,bg[[2,j]],gc[bu[[j]]],ga[bu[[j]]]},
	Frame->True,PlotRange->graphdim,
	AspectRatio->graphratio],{j,2,jt}]
Do[Show[{g0,bg[[3,j]],gc[bu[[j]]],ga[bu[[j]]]},
	Frame->True,PlotRange->graphdim,
	AspectRatio->graphratio],{j,3,jt}]
Do[Show[{g0,bg[[4,j]],gc[bu[[j]]],ga[bu[[j]]]},
	Frame->True,PlotRange->graphdim,
	AspectRatio->graphratio],{j,4,jt}]*)

Do[Module[{z},z=Min[j,jmax];
	Show[{g0,Table[bg[[p,j]],{p,1,z}],
	gc[bu[[j]]],ga[bu[[j]]]},
	Frame->True,PlotRange->graphdim,
	AspectRatio->graphratio]],{j,1,jt}]

ClearAll[a1,a2,onesurf,n,ss,surftot,npart,ellips]
(*bi[[p,j,i]],p est le numéro du tour injecté jusqu'à jmax,
	j est le tour dans la machine de p jusqu'à nt, 
	i est chacune des nx ellipses du faisceau
	Cette partie calcule les efficacités d'injection*)
onesurf[a1_,a2_,g_]:=Abs[a2[[1]]*(a1[[2]]-g[[2]])+a1[[1]]*(g[[2]]-a2[[2]])+g[[1]]*(a2[[2]]-a1[[2]])];
ss[pol_]:=Module[{n,surf,g},n=Length[pol];g={Sum[pol[[i,1]],{i,1,n}]/n,Sum[pol[[i,2]],{i,1,n}]/n};
	surf=onesurf[pol[[1]],pol[[n]],g];
	Do[surf=surf+onesurf[pol[[i]],pol[[i+1]],g],{i,1,n-1}];surf/2]
surftot[ellips_]:=Module[{ssd},ssd=Table[ss[ellips[[i]]],{i,1,nx}];Do[ssd[[i]]=ssd[[i]]-ssd[[i+1]],{i,1,nx-1}];ssd]

npart[ellips_]:=Module[{surface},surface=surftot[ellips];Sum[surface[[i]]*(Exp[-2*((nx+1-i)*dx/x)^2]),{i,1,nx}]]
npartlinac=npart[ztot]


reste=Table[{j,npart[bi[[p,j]]]/npartlinac},{p,1,jmax},{j,1,jt}];

(* cette partie dessine les le nombre de particules injectées sachant que un faisceau linac sur un tour représente 1 unité...c.a.d. si 12 tours sont injectés et on trouve 8 unités cela veut dir qu'il y a 8 tours efficaces injectés ou 8/12=66% d'efficacité*)
trev=1;
restep=Table[Graphics[Line[reste[[p]]]],{p,1,jmax}]
accumul=Table[{0,0},{k,1,2*jt}]
Do[accumul[[2*j-1]]={(j-1)*trev,Sum[reste[[p,j,2]],{p,1,jmax}]};
	accumul[[2j]]={j*trev,Sum[reste[[p,j,2]],{p,1,jmax}]},{j,1,jt}]
Show[restep,Frame->True]
ListPlot[accumul,PlotJoined->True,PlotStyle->{RGBColor[0,0,1],Thickness[0.01]},Frame->True,GridLines->Automatic]