使用数值积分时实施积分限制

Question

我正在使用数值积分来模拟粒子通过非均匀磁场的轨迹。我特别使用了马尔可夫链蒙特卡罗算法 (Metropolis-Hastings)，它允许我计算模型数据以拟合实际粒子的数据。我遇到的问题是我想一次整合单个粒子，因为有时拟合会覆盖其他粒子轨迹：

注：这种情况说明了两个粒子（一个反粒子和一个粒子）。您可以看到拟合在另一个粒子的起点（向右移动）刚好结束。

在这种情况下，我在 z = 337 左右开始积分，但我希望积分在 z = 550 左右停止，因为这是对创建的起源。我试图在集成中引入一个 break 语句，以便集成在对创建的起点处停止，如下所示：

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

但是，这个 if 语句永远不会被执行，因为 zmodel[index_max] 在任何时候都不会等于 zmax。是否有另一种方法可以在执行数值积分时设置限制，允许单独积分每组数据？