import sympy as sp
from math import radians
import sympy.physics.mechanics as me
from pydy.system import System
import numpy as np

f = sp.Symbol('f')
t = sp.Symbol('t')
x, y = me.dynamicsymbols('x y')          # cartesian position
v, gamma = me.dynamicsymbols('v gamma')  # polar velocity
xd, yd       = me.dynamicsymbols('x y', 1)
kd_eqs       = [ xd - v*sp.cos(gamma), 
                 yd - v*sp.sin(gamma) ]

BaseFrame = me.ReferenceFrame('BaseFrame')
origin  = me.Point('origin')
bodycentre  = origin.locatenew('bodycentre', x*BaseFrame.x + y*BaseFrame.y)

# either one will do
bodycentre.set_vel(BaseFrame, bodycentre.pos_from(origin).diff(t, BaseFrame))
#bodycentre.set_vel(BaseFrame, v*sp.cos(gamma)*BaseFrame.x + v*sp.sin(gamma)*BaseFrame.y)

Izz = me.outer(BaseFrame.z, BaseFrame.z)*10
body = me.RigidBody(name="body", masscenter=bodycentre, 
                       frame=BaseFrame, mass=10, inertia=(Izz, bodycentre))

forces = [ (bodycentre, f*BaseFrame.x) ]
bodies = [ body ]

KM = me.KanesMethod(BaseFrame, q_ind=[x,y], u_ind=[v,gamma], kd_eqs=kd_eqs)
fr, frstar = KM.kanes_equations(bodies, forces)

times = np.linspace(0.0, 1, 101)   # 0.01 second increments
sys = System(KM,
             constants={f:0},  # hardcoding f=0 causes issue #427
             initial_conditions={x:0, y:0, v:5, gamma:radians(90) },
             times=times)

y1 = sys.integrate()
print(times[:2])
y1[:2]