Blatt 1 - experienced

Stella Weruth, Jan Scheunemann, Sven Bollweg

Simulating a Cyclotron

basic functions and formulas

cfix = 3.*(10^8);

The time dependent voltage of the radio frequency

UHF[t_, U0_, ωRF_]=U0*Cos[ωRF t]

get β written with γ

β=β/.Solve[γ== (1/(Sqrt[1-(β^2)])), β][[2]]

the Lorentz force

FL=q β c B

Centripetal force attending relativistic particles

FZ=((γ m (β^2) (c^2))/r)

get the formula to calculate the radius in the cyclotron for a given velocity

rZtemp= r/.Flatten[Solve[FL==FZ, r]]

make a function from this formula

rZ[γ2_, B2_, q2_, m2_]=rZtemp/.{γ->γ2, B->B2, q->q2, m->m2}

rZ[γ,B,q,m]

Duration for a half revolution

tTemp[γ_, B_, q_, m_]=((π rZ[γ,B, q, m])/(β c))

A function to calculate the duration of a half revolution

tFunction[γ_, B_, q_, m_]=((m π γ)/(B q))

A function to calculate γ (the velocity) for a given kinetic energy

γFunction[Ekin_, m_]=(Ekin/(m (c^2)))+1

the cyclotron module

A module to simulate the cyclotron.

Cyclotron[nmax_, q_, ωRF_, U0_, m_, B_]:=Block[{Ekin=0, n=0, t=0, tn, γn, rZn, ωn, Un, backwards=1, EkinList={}, rList={}, UList={},c=3.*(10^8)},
(*repeat all calculations for nmax half revolutions*)
For[n=0, n<nmax, n++,
(*calculate the accelerating voltage*)
Un=UHF[t,U0, ωRF] *(((-1))^n);
(* increase the energy*)
Ekin= Ekin+q*Un*backwards;
(*we still want to have a positive energy, if the particle beginns to travel backwards*)
If[Ekin<0,
backwards=-backwards;
Ekin=-Ekin;];
(*calculate the γ for the new energy*)
γn=γFunction[Ekin, m];
(*and the new time for a half revolution*)
tn= tFunction[γn, B, q, m];
(*update the time for the next acceleration*)
t= t+tn;
(*calculate the radius*)
rZn=rZ[γn, B, q, m] ;
(*calculate the cyclotron frequency for this revolution*)
ωn=(π/tn);
(*append the calculated variables to lists*)
EkinList=Append[EkinList, Ekin];
rList=Append[rList, rZn];
UList=Append[UList, Un]];
Return[{EkinList, rList, UList}]]

some numbers and results

mp=9.38272*(10^8)/(c^2)

qP=1

U0P=100000

Bp=1

the frequency for a particle at rest

ωRFP=(π/((mp π)/(Bp*qP)))

results=Cyclotron[400, qP, ωRFP, U0P,mp, Bp];

?ListPlot

ListPlot[results[[1]], AxesLabel->{"n", "Ekin/eV"}]

ListPlot[results[[2]], AxesLabel->{"n", "r/m"}]

ListPlot[results[[3]], AxesLabel->{"n", "U/V"}]

the maximal energy in eV

Max[results[[1]]]

the maximal radius

Max[results[[2]]]

find the best ωHF

A module to find the maximal kinetic energy in eV for particles in a given cyclotron

maxEkin[nmax_, q_, ωRF_, U0_, m_, B_]:=Module[{res},
res=Cyclotron[nmax, q, ωRF, U0, m, B];
Max[res[[1]]]*6.242*(10^18)]

maxEkin[400, qP, ωRFP, U0P,mp, Bp]

A list with possible values for ωHF in units of ωRFP

ωList=Range[0.99, 1.01, 0.0005]

calculate the maximum energy for all possible ωHF values

maxEkinList=maxEkin[400, qP, #, U0P,mp, Bp]&/@ (ωList*ωRFP)

ListPlot[Thread[{ωList, maxEkinList}], AxesLabel->{"((ωHF)/(ωRFP))", "Ekin/eV"}]

maxmaxEkin=Max[maxEkinList]

maxPosition=Flatten[Position[maxEkinList, maxmaxEkin]][[1]]

ωfactor=ωList[[maxPosition]]

I had the idea to use the following, but FindMaximum does not work or is very slow

Plot[maxEkin[400, qP, x*ωRFP, U0P,mp, 1], {x, 0.99, 1.01}, AxesLabel->{"((ωHF)/(ωRFP))", "Ekin/eV"}]

(*FindMaximum[maxEkin[400, qP, x*ωRFP, U0P,mp, 1], {x, 0.99}]*)

results for the best ωHF

results2=Cyclotron[400, qP, ωfactor ωRFP, U0P,mp, Bp];

ListPlot[results2[[1]], AxesLabel->{"n", "Ekin/eV"}]

ListPlot[results2[[2]], AxesLabel->{"n", "r/m"}]

ListPlot[results2[[3]], AxesLabel->{"n", "U/V"}]

Max[results2[[1]]]

Max[results2[[2]]]

Comparison of both

ListPlot[{results[[1]], results2[[1]]}, AxesLabel->{"n", "Ekin/eV"}]

ListPlot[{results[[2]], results2[[2]]}, AxesLabel->{"n", "r/m"}]

ListPlot[{results[[3]], results2[[3]]}, AxesLabel->{"n", "U/V"}]