pycollimate.py 文件源码

def check_mcs_step_septum(septum, proton, step_in):
    """
    From K2 scattering routine. It is basically treatment of edge effect inside the septum jaws
    @param septum: septum object
    @param proton: Proton object
    @param step_in: Random step chosen for this iteration
    """

    theta_0 = 13.6e-3 / (proton.mom_p) / np.sqrt(septum.rad_length) * np.sqrt(step_in) * (1. + 0.038 * np.log(step_in / septum.rad_length))



    x_temp = proton.x - septum.aper_u if proton.x <= (septum.aper_u + septum.thickness / 2.) else proton.x - septum.aper_u - septum.thickness

    a = 9 * theta_0 ** 2 / 3.
    b = -1 * proton.x_p ** 2
    c = -2 * proton.x_p * x_temp
    d = -x_temp ** 2
    # print 'a', a, 'b', b, 'c', c, 'd', d
    # test = np.linspace(-10, 10, 100)
    # test_y = a * test ** 3 + b * test ** 2 + test * c + d
    # plt.plot(test, test_y, 'lime')
    sol = np.roots([a, b, c, d])
    # print sol
    sol_re = [100000.]
    for i in sol:
        if np.imag(i) == 0:
            sol_re.append(float(np.real(i)))
    min_sol = min(sol_re)
    if min_sol < step_in:
        if min_sol < septum.rad_length * 0.1e-2:
            return septum.rad_length * 0.1e-2
        else:
            return min(sol_re)
    else:
        return step_in