Implementing Integration Limits when using Numerical Integration

Question

I'm using numerical integration to model the trajectory of a particle passing through a non-uniform magnetic field. I'm specifically using a Markov Chain Monte Carlo algorithm (Metropolis-Hastings) that allows me to calculate model data in order to fit the actual particle's data. The issue that I have is that I want to integrate individual particles at a single time, since sometimes the fit covers the other particles trajectory:

Note: This situation illustrates two particles (an anti-particle and particle). You can see that the fit ends at just over the beginning of the origin of the other particle (travelling rightwards.)

In this situation I begin the integration at about z = 337, but I want the integration stop at around z = 550, as this is the origin of the pair creation. I have tried to introduce a break statement into the integration so the integration stops at the origin of pair creation, like so:

def evaluation(theta,phi,E,xi,yi,zi):  ### For creating model/experimental data

initial_vel = BROH(E)[0]
gamma_2 = BROH(E)[2]
relative_mass = BROH(E)[3]

first_x = np.zeros(len(actual_x))
first_y = np.zeros(len(actual_y))
first_z = np.zeros(len(actual_z))

xmodel = np.zeros(len(actual_x))   ### Store model data here
ymodel = np.zeros(len(actual_y))
zmodel = np.zeros(len(actual_z))

velocity_x = np.zeros(len(actual_x))  ### Store velocity values to calculate subsequent x,y,z model data
velocity_y = np.zeros(len(actual_y))
velocity_z = np.zeros(len(actual_z))

Bx = np.zeros(len(actual_x))
By = np.zeros(len(actual_y))
Bz = np.zeros(len(actual_z))

first_x[0] = xi         ### Initial guesses for x,y,z
first_y[0] = yi
first_z[0] = zi

velocity_x[0] = initial_vel*np.sin(theta)*np.cos(phi)  ### Initial values for velocities
velocity_y[0] = initial_vel*np.sin(theta)*np.sin(phi)
velocity_z[0] = initial_vel*np.cos(theta)

index = 0
for i in range(len(actual_x) - 1):  ### Loop over experimental/model trajectory
    
    zbzero = zradius[2][0] #for evemt 93  # for event 71 550
    zb = abs(first_z[i] - zbzero)
    if zb > 1000:
        zb = 1000
    
    global Qcharge
    Qcharge = -1.  #positive or negative charge +1 or -1 
    Bz = 1678.5 + 0.080008*zb - 0.019289*zb**2 + 1.3946e-5*zb**3 + 3.0161e-8*zb**4
    Bz = Qcharge*Bz  #for opposite/ normal charge/positive 
    
    Rr = first_x[i]**2 + first_y[i]**2
    if Rr > 1000:
        Rr = 1000
    
    Fact = np.sqrt(Rr) / 40
    Br = Fact*(6.2674e-3 + 0.67562*zb + 1.2677e-4*zb**2 - 6.8352e-6*zb**3 + 6.6604e-9*zb**4)
    Phir = np.arctan2(first_y[i],first_x[i])
    Br = Qcharge*Br #for opposite/ normal charge/positive 
    
    Bx = -2/3*Br*np.cos(Phir)
    By = -2/3*Br*np.sin(Phir)
    
    B_field = np.array([Bx,By,Bz])
    velocity = np.array([velocity_x[i],velocity_y[i],velocity_z[i]])
    cross_product = np.cross(velocity,B_field)
    
    ### Calculate subsequent velocities for model/experimental
    velocity_x[i+1] = velocity_x[i] + const*cross_product[0]*dt / relative_mass
    velocity_y[i+1] = velocity_y[i] + const*cross_product[1]*dt / relative_mass
    velocity_z[i+1] = velocity_z[i] + const*cross_product[2]*dt / relative_mass  

    first_x[i+1] = first_x[i] + velocity_x[i]*dt + 0.5*const*cross_product[0]*dt**2 / relative_mass   
    first_y[i+1] = first_y[i] + velocity_y[i]*dt + 0.5*const*cross_product[1]*dt**2 / relative_mass  
    first_z[i+1] = first_z[i] + velocity_z[i]*dt + 0.5*const*cross_product[2]*dt**2 / relative_mass
    
    if first_x[i+1] > -150 and first_x[i+1] < 150:
        if first_y[i+1] > -150 and first_y[i+1] < 150:
            if first_z[i+1] > 0 and first_z[i+1] < 1000:
                
                global index_max
                index = index + 1
                xmodel[index] = first_x[i+1] + 0.5*const*cross_product[0]*dt**2 / relative_mass 
                ymodel[index] = first_y[i+1] + 0.5*const*cross_product[1]*dt**2 / relative_mass  
                zmodel[index] = first_z[i+1] + 0.5*const*cross_product[2]*dt**2 / relative_mass
                index_max = index
                
    if zmodel[index_max] == zmax:
        break
                
return xmodel[1:index_max], ymodel[1:index_max], zmodel[1:index_max], index_max

However, this if statement never gets executed as zmodel[index_max] at no point ever equals zmax. Is there another method to set the limits when performing numerical integration that allows for each set of data to be integrated separately?

Answer 1

zmax and zmodel[index_max] are both floating point numbers and trying to execute a statement on the basis that these are equal is generally a bad idea. Further, your model data for the fit doesn't guarantee that at any point will equal the experimental data, regardless if they are integers or floats. Instead, it would be ideal to replace your if-break statement with:

if zmodel[index_max] >= zmax:
   break