Storage rings are used in particle physics for the storing of charged particles prior to their use in other apparatus, for example in a particle collider. Instead of having to build a storage ring, it is useful to simulate one using a computer. The following program goes through the necessary steps (using a Monte Carlo method) to simulate a storage ring and includes plots and other methods to qualify the results obtained.

The first steps in producing the storage ring simulation involve giving particles random, but uniformly distributed, x and y coordinates, as well divergences, x' and y'. We then want to be able to change these values for each 'step' along the storage ring that the particle travels. These values can be calculated through the use of the following matrices:

$$\left[\begin{array}{c}x_{2}\\x'_{2}\end{array}\right] = \left[\begin{array}{cc} \cos\left(\frac{z}{r_{0}}\sqrt{(1-n)}\right) & \frac{r_{0}}{\sqrt{(1-n)}}\sin\left(\frac{z}{r_{0}}\sqrt{(1-n)}\right) \\ -\frac{\sqrt{(1-n)}}{r_{0}}\sin\left(\frac{z}{r_{0}}\sqrt{(1-n)}\right) & \cos\left(\frac{z}{r_{0}}\sqrt{(1-n)}\right) \end{array}\right] \left[\begin{array}{c}x_{1}\\x'_{1}\end{array}\right]$$$$\left[\begin{array}{c}y_{2}\\y'_{2}\end{array}\right] = \left[\begin{array}{cc} \cos\left(\frac{z}{r_{0}}\sqrt{n}\right) & \frac{r_{0}}{\sqrt{n}}\sin\left(\frac{z}{r_{0}}\sqrt{n}\right) \\ -\frac{\sqrt{n}}{r_{0}}\sin\left(\frac{z}{r_{0}}\sqrt{n}\right) & \cos\left(\frac{z}{r_{0}}\sqrt{n}\right) \end{array}\right] \left[\begin{array}{c}y_{1}\\y'_{1}\end{array}\right]$$

We must also give these four parameters minimum and maximum values if the particle is to be stable inside of the storage ring.

The following code demonstrates this, and takes us to checkpoint 2:

from numpy import sin, cos, pi, random, sqrt

#Time how long the code takes to run. This is always an important thing to do to measure efficiency.
import time
start = time.time()

# r is the radius of the storage ring, in this case 2 metres. z is the step size.
r = 2
z = 0.1*pi

# At this stage we can assign the field index to be 0.5 for simplicity
n = 0.5

# Define Minimum and Maximum values of x,y. In this case, x and y must both lie between +3cm and -3cm.
xMin = -0.03
xMax = 0.03
yMin = -0.03
yMax = 0.03
# Determine Maximum and Minimum values of divergence. In this case the divergence must not be greater than the absolute value of
# 3cm/r
dxMax = 0.03/r
dyMax = 0.03/r
dxMin = - dxMax
dyMin = - dyMax

# matrix: We now define a function to calculate the coefficients of the matrix elements. This is easier than writing all of the equations
# by hand.
def matrix(n):
    AX=cos(z*sqrt(1-n)/r)
    BX=sin(z*sqrt(1-n)/r)*r/sqrt(1-n)
    CX=-sin(z*sqrt(1-n)/r)*sqrt(1-n)/r
    
    AY=cos(z*sqrt(n)/r)
    BY=sin(z*sqrt(n)/r)*r/sqrt(n)
    CY=-sin(z*sqrt(n)/r)*sqrt(n)/r
    return AX,BX,CX,AY,BY,CY

# Step: A function to simulate the a particle going along one step. This will use values of x, y, x' and y' as well as the
# coefficients calculated in matrix.
def Step(xCoord, yCoord, xDiv, yDiv, n, AX,BX,CX,AY,BY,CY):   
    newxCoord = AX*xCoord + BX*xDiv
    newxDiv = CX*xCoord + AX*xDiv
    newyCoord = AY*yCoord + BY*yDiv
    newyDiv = CY*yCoord + AY*yDiv
    xCoord = newxCoord
    xDiv = newxDiv
    yCoord = newyCoord
    yDiv = newyDiv
    return (xCoord, yCoord, xDiv, yDiv)

# Loop: A function to see if each step is 'successful'
def Loop(RevNumber, StepNumber, xCoord, yCoord, xDiv, yDiv,n, AX,BX,CX,AY,BY,CY):
# We must initialise the number of steps that the particle has taken to 0 for each particle.
    Steps=0
    Success=False
    return (Success)
    
# This checks to see if the particle's motion along the step is valid. The Step function is called for the particle, the new
# values of x, y, x' and y' are checked to see if they are still valid. If at any point they are not valid the particle stops
# being simulated (break). Otherwise the particle is asked to complete another step. This loop ends if 2000 steps have been
# taken, as this is the point at which the particle is considered stable.
    TN=RevNumber*StepNumber
    for TotalNumber in range(0,TN):
        xCoord,yCoord,xDiv, yDiv = Step(xCoord,yCoord,xDiv, yDiv,n, AX,BX,CX,AY,BY,CY)
        if xCoord<xMin or xCoord>xMax or xDiv<dxMin or xDiv>dxMax or yCoord<yMin or yCoord>yMax or yDiv<dyMin or yDiv>dyMax :
            break 
        else:
            Steps += 1
    if Steps == 2000:
        Success = True
    return (Success)


# Now that our functions are defined, we can start a simulation!          
AX,BX,CX,AY,BY,CY = matrix(n)
# How many particles do we want to simulate?
ParticleNumber = 10000

#For each particle we must give it a random x, y, x' and y'. This is done through random.uniform
for j in range(0,ParticleNumber):
    xCoord = random.uniform(-0.03,0.03)
    yCoord = random.uniform(-0.03,0.03)
    xDiv = random.uniform(dxMin, dxMax)
    yDiv = random.uniform(dyMin, dyMax)
    Result = Loop(50,40, xCoord, yCoord, xDiv, yDiv,n,AX,BX,CX,AY,BY,CY)
    

end = time.time()
print('The time taken to run this program was:',end - start)

The time taken to run this program was: 0.042661428451538086

Now that particles are travelling through the storage ring, it would be useful to know how many particles are successfully completing the requirement of 2000 steps that we have introduced for us to consider them stable.

The way to do this is to measure acceptance, the ratio of stable orbit particles to the total number of particles introduced to the ring.

It would also be useful to make a note of how many particles are stable relative to x, x' and how many are stable to y, y'. This is done through making graphs of acceptance in x and y.

It will be useful to also make a note of at what point (step) it is into the simulation the particles that become unstable do so.

The following code will provide graphs of these and take us to checkpoint 4:

from numpy import sin, cos, pi, random, sqrt

# We must also now import a graphing method.
from matplotlib import pyplot as plt

r = 2
z = pi

import time
start = time.time()

xMin = -0.03
xMax = 0.03
yMin = -0.03
yMax = 0.03
dxMax = 0.03/r
dyMax = 0.03/r
dxMin = - dxMax
dyMin = - dyMax

def matrix(n):
    AX=cos(z*sqrt(1-n)/r)
    BX=sin(z*sqrt(1-n)/r)*r/sqrt(1-n)
    CX=-sin(z*sqrt(1-n)/r)*sqrt(1-n)/r
    
    AY=cos(z*sqrt(n)/r)
    BY=sin(z*sqrt(n)/r)*r/sqrt(n)
    CY=-sin(z*sqrt(n)/r)*sqrt(n)/r
    return AX,BX,CX,AY,BY,CY

def Step(xCoord, yCoord, xDiv, yDiv, n, AX,BX,CX,AY,BY,CY):   
    newxCoord = AX*xCoord + BX*xDiv
    newxDiv = CX*xCoord + AX*xDiv
    newyCoord = AY*yCoord + BY*yDiv
    newyDiv = CY*yCoord + AY*yDiv
    xCoord = newxCoord
    xDiv = newxDiv
    yCoord = newyCoord
    yDiv = newyDiv
    return (xCoord, yCoord, xDiv, yDiv)

def Loop(RevNumber, StepNumber, xCoord, yCoord, xDiv, yDiv,n, AX,BX,CX,AY,BY,CY):
    Steps=0
    Success=False
# To check whether the particle is stable in x, y we must introduce more booleans.
    AccY = False
    AccX = False
    
    TN=RevNumber*StepNumber
    for TotalNumber in range(0,TN):
        xCoord,yCoord,xDiv, yDiv = Step(xCoord,yCoord,xDiv, yDiv,n, AX,BX,CX,AY,BY,CY)
        if xCoord<xMin or xCoord>xMax or xDiv<dxMin or xDiv>dxMax or yCoord<yMin or yCoord>yMax or yDiv<dyMin or yDiv>dyMax:
            
# At this point the particle has become unstable, as it has satisfied one of the above criteria. To determine Acceptance in x,y
# we must check if the failure was to do with x, if so, y is successful, if not x is successful.

            if xCoord<xMin or xCoord>xMax or xDiv<dxMin or xDiv>dxMax:
                AccY = True
            else:
                AccX = True
            break 
        else:
            Steps += 1
    if Steps == TN:
        Success = True
        AccY = True
        AccX = True
    return (Success, AccX, AccY, Steps)

# Create lists for field index and acceptance for use in graphing. Field index, Acceptance, Acceptance in x, Acceptance in y.
nList = []
AList = []
XList = []
YList = []

# To figure out which step the particles that fail do so, we create a list that we can append the values of 'steps' that fail.
StepsList = []
RevNumber = 50
StepNumber = 40
TN = RevNumber*StepNumber


# As we are considering the where the particles fail for both n = 0.1 and n = 0.5, we should use a list (Lost) which contains
# two lists of particles that fail: 1) when considering n = 0.1, 2) when considering n = 0.5.
Lost=[ [], [] ]

# Here we assign the StepsList list to contain the values 0 to 2000. This will provide us with an x-axis.
# We also initialise both of the Lost lists to contain just as many values, but all of the values are set to 0. These will
# serve as our y-axes.
for s in range(0,TN):
    StepsList.append(s)
    Lost[0].append(0)
    Lost[1].append(0)   

# We want to consider the simulation for particles with varying field index, from 0 to 1, with increments of 0.05. As n = 0
# will not work nicely with our functions for particle motion we use n = 0.01 here instead.
for i in range(0, 20):
    if i == 0:
        n = 0.01
    else:
        n = i/20
    
# We fill nList with the values of n as we loop through them. We also initialise the acceptances to 0.    
    nList.append(n)
    AX,BX,CX,AY,BY,CY = matrix(n)
    SuccessNumber = 0
    AcceptXNo = 0
    AcceptYNo = 0
    ParticleNumber = 10000

    
    for j in range(0,ParticleNumber):
        xCoord = random.uniform(-0.03,0.03)
        yCoord = random.uniform(-0.03,0.03)
        xDiv = random.uniform(dxMin, dxMax)
        yDiv = random.uniform(dyMin, dyMax)
        Result = Loop(RevNumber, StepNumber, xCoord, yCoord, xDiv, yDiv,n,AX,BX,CX,AY,BY,CY)

# If, for a particle, our Loop function returns that the particle is successful, we add one to the number of successes.
# Similarly for the successes in x and successes in y.
        if Result[0] == True:
            SuccessNumber +=1
        if Result[1] == True:
            AcceptXNo +=1
        if Result[2] == True:
            AcceptYNo +=1
        
# When a particle loops, we return the step to which it progressed. We now assign this value to a variable, Steps, and add one
# to the value stored at the 'Steps-1' index of Lost (for the appropriate list in Lost.) We use Steps-1 to account for the 
# discrepancies in our counting to 2000 from 1, and the computer starting from 0.
        Steps = Result[3]
        if n == 0.1:
            Lost[0][Steps-1] += 1
        if n == 0.5:
            Lost[1][Steps-1] += 1
        
# To stop the graph from showing a spike at Steps = 2000 (where all of the successful particles lie), we assign this value to 0.
    Lost[0][TN-1] = 0
    Lost[1][TN-1] = 0
    
# We next append our lists to contain the correct values, so that we can use them to plot graphs.
    Acceptance = SuccessNumber/ParticleNumber
    AcceptanceX = AcceptXNo/ParticleNumber
    AcceptanceY = AcceptYNo/ParticleNumber
    AList.append(Acceptance)
    XList.append(AcceptanceX)
    YList.append(AcceptanceY)


# Define a function to replace instances of '0' with the string 'nan'. This removes the tail of 0's from our graphs of losses.    
def zero_to_nan(Lostzero):
    """Replace every 0 with 'nan' and return a copy."""
    return [float('nan') if x==0 else x for x in Lostzero] 


# We now plot the graph for Overall Acceptance, along with an appropriate title and axes labels.
plt.plot(nList, AList,'r')
plt.suptitle("Measured values of Acceptance against Field Index (n)")
plt.xlabel("Field Index (n)")
plt.ylabel("Acceptance")

# The next figure will have two subplots, so that we are easily able to compare the acceptance in x and in y.
f, (xlist, ylist) = plt.subplots(2, 1)
xlist.plot(nList, XList, 'r')
ylist.plot(nList, YList, 'r')
xlist.set_title("Comparison of Acceptance in X and Acceptance in Y")
plt.xlabel("Field Index (n)")
plt.ylabel("                                  Acceptance in Y            Acceptance in X")

# This figure will again have two subplots, used to compare the losses at n = 0.1 and n = 0.5. They are plot on an appropriate
# logarithmic axis.
f, (pointone, pointfive) = plt.subplots(2, 1)
pointone.semilogx(StepsList, zero_to_nan(Lost[0]))
pointfive.semilogx(StepsList, zero_to_nan(Lost[1]))
pointone.set_title("Comparison of the Step Number which a Particle becomes Unstable when n = 0.1 (top) and n = 0.5 (bottom)")
plt.xlabel("Step Number (logarithmic)")
plt.ylabel("                                  Number of Occurances of Unstabliity")



end = time.time()
print('The time taken to run this program was:',end - start)

The time taken to run this program was: 218.66129279136658

Now that we have graphs for the overall acceptance, acceptance in x and acceptance in y against field index, it would be useful to compare these with graphs from theory. Theory predicts that the overall acceptance is given by:

$$ \frac{\pi^{2}}{16}\sqrt{n(1-n)}$$

This can be factorised into the contribution to acceptance from x (acceptance in x) and the contribution to acceptance from y (acceptance in y).

The acceptance in x is theorised to be: $$ \frac{\pi}{4}\sqrt{(1-n)}$$

The acceptance in y is theorised to be: $$ \frac{\pi}{4}\sqrt{n}$$

We can superimpose these results onto our graphs we have already obtained and compare them.

It will also be useful to quantify how the graph of acceptance compares with the theoretical prediction. This will be done through the chi-squared method. It is useful to see that the chi-squared coefficient is minimised as we simulate more particles, showing that the results obtained are more accurate if we use more particles.

Another way of thinking of the acceptance is through the concept of acceptance sets. These sets contain the allowed values of x, x' or y, y' for a given field index. These are shown to obey the form of an ellipse and will also be plotted. We will complete checkpoint 6 here.

from numpy import sin, cos, pi, random, sqrt
from matplotlib import pyplot as plt

r = 2
z = 0.1*pi

import time
start = time.time()

xMin = -0.03
xMax = 0.03
yMin = -0.03
yMax = 0.03
dxMax = 0.03/r
dyMax = 0.03/r
dxMin = - dxMax
dyMin = - dyMax

def matrix(n):
    AX=cos(z*sqrt(1-n)/r)
    BX=sin(z*sqrt(1-n)/r)*r/sqrt(1-n)
    CX=-sin(z*sqrt(1-n)/r)*sqrt(1-n)/r
    
    AY=cos(z*sqrt(n)/r)
    BY=sin(z*sqrt(n)/r)*r/sqrt(n)
    CY=-sin(z*sqrt(n)/r)*sqrt(n)/r
    return AX,BX,CX,AY,BY,CY

def Step(xCoord, yCoord, xDiv, yDiv, n, AX,BX,CX,AY,BY,CY):   
    newxCoord = AX*xCoord + BX*xDiv
    newxDiv = CX*xCoord + AX*xDiv
    newyCoord = AY*yCoord + BY*yDiv
    newyDiv = CY*yCoord + AY*yDiv
    xCoord = newxCoord
    xDiv = newxDiv
    yCoord = newyCoord
    yDiv = newyDiv
    return (xCoord, yCoord, xDiv, yDiv)

def Loop(RevNumber, StepNumber, xCoord, yCoord, xDiv, yDiv,n, AX,BX,CX,AY,BY,CY):
    Steps=0
    Success=False
    AccY = False
    AccX = False
    AX,BX,CX,AY,BY,CY = matrix(n)
    TN=RevNumber*StepNumber
    for TotalNumber in range(0,TN):
        xCoord,yCoord,xDiv, yDiv = Step(xCoord,yCoord,xDiv, yDiv,n, AX,BX,CX,AY,BY,CY)
        if xCoord<xMin or xCoord>xMax or xDiv<dxMin or xDiv>dxMax or yCoord<yMin or yCoord>yMax or yDiv<dyMin or yDiv>dyMax :
            if xCoord<xMin or xCoord>xMax or xDiv<dxMin or xDiv>dxMax:
                AccY = True
            else:
                AccX = True
            break 
        else:
            Steps += 1
    if Steps == 2000:
        Success = True
        AccY = True
        AccX = True
    return (Success, AccX, AccY, Steps, TN)

# Here we create more lists for the theoretical values of acceptance we will obtain. We also create three lists for our 
# plots of the x-x' phase space and the y-y' phase space.
nList = []
AList = []
XList = []
YList = []
TList = []
StepsList = []
RevNumber = 50
StepNumber = 40
TN = RevNumber*StepNumber
Lost=[ [], [] ]
TX = []
TY = []
xAccSet = [ [], [], [] ]
xPAccSet = [ [], [], [] ]
yAccSet = [ [], [], [] ]
yPAccSet = [ [], [], [] ]

for s in range(0,TN):
    StepsList.append(s)
    Lost[0].append(0)
    Lost[1].append(0)

# We need to initialise our Chi-squared parameter to 0.
Chi = 0    
    
for i in range(0, 20):
    if i == 0:
        n = 0.01
    else:
        n = i/20
    
    nList.append(n)
    AX,BX,CX,AY,BY,CY = matrix(n)
    SuccessNumber = 0
    AcceptXNo = 0
    AcceptYNo = 0
    ParticleNumber = 10000

    
    for j in range(0,ParticleNumber):
        xCoord = random.uniform(-0.03,0.03)
        yCoord = random.uniform(-0.03,0.03)
        xDiv = random.uniform(dxMin, dxMax)
        yDiv = random.uniform(dyMin, dyMax)
        Result = Loop(RevNumber, StepNumber, xCoord, yCoord, xDiv, yDiv,n,AX,BX,CX,AY,BY,CY)
        if Result[0] == True:
            SuccessNumber +=1
        if Result[1] == True:
            AcceptXNo +=1
        if Result[2] == True:
            AcceptYNo +=1
        
            
        Steps = Result[3]
        if n == 0.1:
            Lost[0][Steps-1] += 1
        if n == 0.5:
            Lost[1][Steps-1] += 1
        
# In order to plot the maximum values of x, x', y and y' for successful particles we must store this information. This is done
# by appending it to the lists we made.
        if Result[0] == True:
            if n == 0.1:
                xAccSet[0].append(xCoord)
                xPAccSet[0].append(xDiv)
                yAccSet[0].append(yCoord)
                yPAccSet[0].append(yDiv)
            if n == 0.5:
                xAccSet[1].append(xCoord)
                xPAccSet[1].append(xDiv)
                yAccSet[1].append(yCoord)
                yPAccSet[1].append(yDiv)
            if n == 0.9:
                xAccSet[2].append(xCoord)
                xPAccSet[2].append(xDiv)
                yAccSet[2].append(yCoord)
                yPAccSet[2].append(yDiv)

    Lost[0][TN-1] = 0
    Lost[1][TN-1] = 0
    Acceptance = SuccessNumber/ParticleNumber
    AcceptanceX = AcceptXNo/ParticleNumber
    AcceptanceY = AcceptYNo/ParticleNumber
    AList.append(Acceptance)
    XList.append(AcceptanceX)
    YList.append(AcceptanceY)
    
# For each value of field index we need to store the theoretical value of acceptance (overall, in x, and in y).
    theory= pi*pi*sqrt(n*(1-n))/16
    theoryx = pi*sqrt(1-n)/4
    theoryy = pi*sqrt(n)/4
    TList.append(theory)
    TX.append(theoryx)
    TY.append(theoryy)
    
# The Chi Squared distribution is given: Sum((Observed-Expected)**2 / Expected) This is calculated here.
    Contribution = (Acceptance - theory)**2 / theory
    Chi += Contribution
    
def zero_to_nan(values):
    """Replace every 0 with 'nan' and return a copy."""
    return [float('nan') if x==0 else x for x in values] 


# Here we have added the theoretical values of acceptance to the graph obtained earlier. We give it a blue dashed line to 
# distinguish it from the solid red line of obtained results.
plt.plot(nList, AList,'r')
plt.plot(nList, TList,'b--')
plt.suptitle("Measured (red) and Theoretical (Blue) values of Acceptance against Field Index (n)")
plt.xlabel("Field Index (n)")
plt.ylabel("Acceptance")

# Similarly, we plot the theoretical values of acceptance in x, y on the subplots obtained earlier.
f, (xlist, ylist) = plt.subplots(2, 1)
xlist.plot(nList, XList, 'r')
xlist.plot(nList, TX, 'b--')
ylist.plot(nList, YList, 'r')
ylist.plot(nList, TY, 'b--')
xlist.set_title("Comparison of Acceptance in X and Acceptance in Y")
plt.xlabel("Field Index (n)")
plt.ylabel("                                  Acceptance in Y            Acceptance in X")

f, (pointone, pointfive) = plt.subplots(2, 1)
pointone.semilogx(StepsList, zero_to_nan(Lost[0]))
pointfive.semilogx(StepsList, zero_to_nan(Lost[1]))
pointone.set_title("Comparison of the Step Number which a Particle becomes Unstable when n = 0.1 (top) and n = 0.5 (bottom)")
plt.xlabel("Step Number (logarithmic)")
plt.ylabel("                                  Number of Occurances of Unstabliity")


# Figure 4 is the phase space plot for x and x'. We give each value of field index a different colour to differentiate between
# them.
plt.figure(4)
plt.plot(xAccSet[0], xPAccSet[0], 'ko')
plt.plot(xAccSet[1], xPAccSet[1], 'ro')
plt.plot(xAccSet[2], xPAccSet[2], 'bo')
plt.suptitle("Acceptance Sets in x for n = 0.1 (black), n = 0.5 (yellow), n = 0.9 (blue)")
plt.xlabel("x")
plt.ylabel("x'")

# Figure 5 is similar to figure 4, but shows the phase space plot for y, y'.
plt.figure(5)
plt.plot(yAccSet[2], yPAccSet[2], 'bo')
plt.plot(yAccSet[1], yPAccSet[1], 'ro')
plt.plot(yAccSet[0], yPAccSet[0], 'ko')
plt.suptitle("Acceptance Sets in y for n = 0.1 (black), n = 0.5 (yellow), n = 0.9 (blue)")
plt.xlabel("y")
plt.ylabel("y'")


plt.show()


print('The Chi Squared "Goodness of Fit" was found to be', Chi)

end = time.time()
print('The time taken to run this program was:',end - start)

The Chi Squared "Goodness of Fit" was found to be 0.00108929487446
The time taken to run this program was: 214.6746861934662

In reality, Storage Rings need a gap in them for the injection and extraction of particles. We will now implement this gap, which will have a size of 10% of the ring and will start about halfway around the ring.

Whilst the particles travel along the steps associated with this gap the field index will be changed. This will effect the function used to change the values of x, x', y and y' when the particle travels along a step. The calculation of new matrix elements will be implemented, using the old matrices identities, but with n = 0. This is given:

$$\left[\begin{array}{c}x_{2}\\x'_{2}\end{array}\right] = \left[\begin{array}{cc} \cos\left(\frac{z}{r_{0}}\right) & r_{0}\sin\left(\frac{z}{r_{0}}\right) \\ -\frac{1}{r_{0}}\sin\left(\frac{z}{r_{0}}\right) & \cos\left(\frac{z}{r_{0}}\right) \end{array}\right] \left[\begin{array}{c}x_{1}\\x'_{1}\end{array}\right]$$$$\left[\begin{array}{c}y_{2}\\y'_{2}\end{array}\right] = \left[\begin{array}{cc} 1 & z \\ 0 & 1 \end{array}\right] \left[\begin{array}{c}y_{1}\\y'_{1}\end{array}\right]$$

The Sin series expansion was used for the second term, to avoid division by 0.

from numpy import sin, cos, pi, random, sqrt
from matplotlib import pyplot as plt

r = 2
z = 0.1*pi

import time
start = time.time()

xMin = -0.03
xMax = 0.03
yMin = -0.03
yMax = 0.03
dxMax = 0.03/r
dyMax = 0.03/r
dxMin = - dxMax
dyMin = - dyMax

def matrix(n):
# Here we implement the new values of the matrix coefficients for if the particle is undergoing a step within the gap.
    if n == 0:
        AX = cos(z/r)
        BX = r*sin(z/r)
        CX = -sin(z/r)/r
    
        AY = 1
        BY = z
        CY = 0
    else:
        AX=cos(z*sqrt(1-n)/r)
        BX=sin(z*sqrt(1-n)/r)*r/sqrt(1-n)
        CX=-sin(z*sqrt(1-n)/r)*sqrt(1-n)/r
        
        AY=cos(z*sqrt(n)/r)
        BY=sin(z*sqrt(n)/r)*r/sqrt(n)
        CY=-sin(z*sqrt(n)/r)*sqrt(n)/r
    return AX, BX, CX, AY, BY, CY

def Step(xCoord, yCoord, xDiv, yDiv, n, AX,BX,CX,AY,BY,CY):
    #AX,BX,CX,AY,BY,CY = matrix(n)
    newxCoord = AX*xCoord + BX*xDiv
    newxDiv = CX*xCoord + AX*xDiv
    newyCoord = AY*yCoord + BY*yDiv
    newyDiv = CY*yCoord + AY*yDiv
    xCoord = newxCoord
    xDiv = newxDiv
    yCoord = newyCoord
    yDiv = newyDiv
    return (xCoord, yCoord, xDiv, yDiv)

def Loop(RevNumber, StepNumber, xCoord, yCoord, xDiv, yDiv,n, AX,BX,CX,AY,BY,CY):
    Steps=0
    Success=False
    AccY = False
    AccX = False
    RevCounter = 1
    TN=RevNumber*StepNumber
    
    for TotalNumber in range(0,TN):
# Here we want different conditions if the particle is along the gaps for injection/extraction.
        if RevCounter*(StepNumber/2)-1 < Steps < RevCounter*(StepNumber/2)+2  :  
            n = 0
        else:
            n = n
        
        xCoord,yCoord,xDiv, yDiv = Step(xCoord,yCoord,xDiv, yDiv,n, AX,BX,CX,AY,BY,CY)
        if xCoord<xMin or xCoord>xMax or xDiv<dxMin or xDiv>dxMax or yCoord<yMin or yCoord>yMax or yDiv<dyMin or yDiv>dyMax :
            if xCoord<xMin or xCoord>xMax or xDiv<dxMin or xDiv>dxMax:
                AccY = True
            else:
                AccX = True
            break 
        else:
            Steps += 1
            if Steps % StepNumber == 0:
                if RevCounter != 1:
                    RevCounter += 1
    if Steps == TN:
        Success = True
        AccY = True
        AccX = True
    return (Success, AccX, AccY, Steps, TN)


nList = []
AList = []
XList = []
YList = []
TList = []
StepsList = []
RevNumber = 50
StepNumber = 40
TN = RevNumber*StepNumber
Lost=[ [], [] ]
TX = []
TY = []
xAccSet = [ [], [], [] ]
xPAccSet = [ [], [], [] ]
yAccSet = [ [], [], [] ]
yPAccSet = [ [], [], [] ]

for s in range(0,TN):
    StepsList.append(s)
    Lost[0].append(0)
    Lost[1].append(0)


Chi = 0    
    
for i in range(0, 20):
    if i == 0:
        n = 0.01
    else:
        n = i/20
    
    nList.append(n)
    AX,BX,CX,AY,BY,CY = matrix(n)
    SuccessNumber = 0
    AcceptXNo = 0
    AcceptYNo = 0
    ParticleNumber = 10000

    
    for j in range(0,ParticleNumber):
        xCoord = random.uniform(-0.03,0.03)
        yCoord = random.uniform(-0.03,0.03)
        xDiv = random.uniform(dxMin, dxMax)
        yDiv = random.uniform(dyMin, dyMax)
        Result = Loop(RevNumber, StepNumber, xCoord, yCoord, xDiv, yDiv,n,AX,BX,CX,AY,BY,CY)
        if Result[0] == True:
            SuccessNumber +=1
        if Result[1] == True:
            AcceptXNo +=1
        if Result[2] == True:
            AcceptYNo +=1
        
            
        Steps = Result[3]
        if n == 0.1:
            Lost[0][Steps-1] += 1
        if n == 0.5:
            Lost[1][Steps-1] += 1
        

        if Result[0] == True:
            if n == 0.1:
                xAccSet[0].append(xCoord)
                xPAccSet[0].append(xDiv)
                yAccSet[0].append(yCoord)
                yPAccSet[0].append(yDiv)
            if n == 0.5:
                xAccSet[1].append(xCoord)
                xPAccSet[1].append(xDiv)
                yAccSet[1].append(yCoord)
                yPAccSet[1].append(yDiv)
            if n == 0.9:
                xAccSet[2].append(xCoord)
                xPAccSet[2].append(xDiv)
                yAccSet[2].append(yCoord)
                yPAccSet[2].append(yDiv)

    Lost[0][TN-1] = 0
    Lost[1][TN-1] = 0
    Acceptance = SuccessNumber/ParticleNumber
    AcceptanceX = AcceptXNo/ParticleNumber
    AcceptanceY = AcceptYNo/ParticleNumber
    AList.append(Acceptance)
    XList.append(AcceptanceX)
    YList.append(AcceptanceY)
    
# For each value of field index we need to store the theoretical value of acceptance (overall, in x, and in y).
    theory= pi*pi*sqrt(n*(1-n))/16
    theoryx = pi*sqrt(1-n)/4
    theoryy = pi*sqrt(n)/4
    TList.append(theory)
    TX.append(theoryx)
    TY.append(theoryy)
    
# The Chi Squared distribution is given: Sum((Observed-Expected)**2 / Expected) This is calculated here.
    Contribution = (Acceptance - theory)**2 / theory
    Chi += Contribution
    
def zero_to_nan(values):
    """Replace every 0 with 'nan' and return a copy."""
    return [float('nan') if x==0 else x for x in values] 


# Here we have added the theoretical values of acceptance to the graph obtained earlier. We give it a blue dashed line to 
# distinguish it from the solid red line of obtained results.
plt.plot(nList, AList,'r')
plt.plot(nList, TList,'b--')
plt.suptitle("Measured (red) and Theoretical (Blue) values of Acceptance against Field Index (n)/n The Measured values have a gap of 10% inserted halfway around the ring")
plt.xlabel("Field Index (n)")
plt.ylabel("Acceptance")

# Similarly, we plot the theoretical values of acceptance in x, y on the subplots obtained earlier.
f, (xlist, ylist) = plt.subplots(2, 1)
xlist.plot(nList, XList, 'r')
xlist.plot(nList, TX, 'b--')
ylist.plot(nList, YList, 'r')
ylist.plot(nList, TY, 'b--')
xlist.set_title("Comparison of Acceptance in X and Acceptance in Y")
plt.xlabel("Field Index (n)")
plt.ylabel("                                  Acceptance in Y            Acceptance in X")

f, (pointone, pointfive) = plt.subplots(2, 1)
pointone.semilogx(StepsList, zero_to_nan(Lost[0]))
pointfive.semilogx(StepsList, zero_to_nan(Lost[1]))
pointone.set_title("Comparison of the Step Number which a Particle becomes Unstable when n = 0.1 (top) and n = 0.5 (bottom)")
plt.xlabel("Step Number (logarithmic)")
plt.ylabel("                                  Number of Occurances of Unstabliity")


# Figure 4 is the phase space plot for x and x'. We give each value of field index a different colour to differentiate between
# them.
plt.figure(4)
plt.plot(xAccSet[0], xPAccSet[0], 'ko')
plt.plot(xAccSet[1], xPAccSet[1], 'ro')
plt.plot(xAccSet[2], xPAccSet[2], 'bo')
plt.suptitle("Acceptance Sets in x for n = 0.1 (black), n = 0.5 (red), n = 0.9 (blue)")
plt.xlabel("x")
plt.ylabel("x'")

# Figure 5 is similar to figure 4, but shows the phase space plot for y, y'.
plt.figure(5)
plt.plot(yAccSet[2], yPAccSet[2], 'bo')
plt.plot(yAccSet[1], yPAccSet[1], 'ro')
plt.plot(yAccSet[0], yPAccSet[0], 'ko')
plt.suptitle("Acceptance Sets in y for n = 0.1 (black), n = 0.5 (red), n = 0.9 (blue)")
plt.xlabel("y")
plt.ylabel("y'")


plt.show()


print('The Chi Squared "Goodness of Fit" for the overall Acceptance was found to be', Chi)

end = time.time()
print('The time taken to run this program was:',end - start)

The Chi Squared "Goodness of Fit" for the overall Acceptance was found to be 0.00144363988245
The time taken to run this program was: 302.8053352832794