TrajectoryFactory.py 文件源码

def trajectory_from_magnetic_field_method_ODE(self, source,t=None):
        gamma = source.Lorentz_factor()

        B=source.magnetic_field
        #   trajectory =
        #   [t........]
        #   [ X/c......]
        #   [ Y/c ......]
        #   [ Z/c ......]
        #   [ Vx/c .....]
        #   [ Vy/c .....]
        #   [ Vz/c .....]
        #   [ Ax/c .....]
        #   [ Ay/c .....]
        #   [ Az/c .....]

        if t is None :
            time_calc = source.construct_times_vector(initial_contition=self.initial_condition,Nb_pts=self.Nb_pts)
        else :
            self.Nb_pts=len(t)
            time_calc =t

        trajectory = np.zeros((10, self.Nb_pts))
        trajectory[0]=time_calc

        cst = -codata.e / (codata.m_e * gamma)
        initial_condition_for_ODE=self.initial_condition

        #TODO rtol et a tol modifiable
        #atol=np.array([1e-10,1e-10,1e-10,1e-10,1e-10,1e-10])
        rtol=source.rtol_for_ODE_method()
        atol=source.atol_for_ODE_method()
        res = odeint(func=fct_ODE_magnetic_field,y0=initial_condition_for_ODE, t=trajectory[0],
                     args=(cst,B.Bx,B.By,B.Bz),rtol=rtol,atol=atol,mxstep=5000,full_output=True)
        traj = res[0]
        info = res[1]
        print("1 : nonstiff problems, Adams . 2: stiff problem, BDF")
        print(info.get('mused'))
        traj = np.transpose(traj)
        trajectory[4] = traj[0]
        trajectory[5] = traj[1]
        trajectory[6] = traj[2]
        trajectory[1] = traj[3]
        trajectory[2] = traj[4]
        trajectory[3] = traj[5]
        trajectory[7] = - cst * B.By(trajectory[3], trajectory[2],trajectory[1]) * trajectory[6]
        trajectory[9] = cst * B.By(trajectory[3], trajectory[2],trajectory[1]) * trajectory[4]
        T=self.create_from_array(trajectory)
        T.multiply_by((1.0/codata.c))
        return T