<s> from sympy import symbols , Matrix , solve , Poly <EOL> from sympy . physics . mechanics import * <EOL> t , r , m , g , I , J = symbols ( '<STR_LIT>' ) <EOL> Fx , Fy , Fz = symbols ( '<STR_LIT>' ) <EOL> q = dynamicsymbols ( '<STR_LIT>' ) <EOL> qd = [ qi . diff ( t ) for qi in q ] <EOL> u = dynamicsymbols ( '<STR_LIT>' ) <EOL> ud = [ ui . diff ( t ) for ui in u ] <EOL> ud_zero = { udi : <NUM_LIT:0> for udi in ud } <EOL> ua = dynamicsymbols ( '<STR_LIT>' ) <EOL> ua_zero = { uai : <NUM_LIT:0> for uai in ua } <EOL> N = ReferenceFrame ( '<STR_LIT:N>' ) <EOL> A = N . orientnew ( '<STR_LIT:A>' , '<STR_LIT>' , [ q [ <NUM_LIT:0> ] , N . z ] ) <EOL> B = A . orientnew ( '<STR_LIT:B>' , '<STR_LIT>' , [ q [ <NUM_LIT:1> ] , A . x ] ) <EOL> C = B . orientnew ( '<STR_LIT:C>' , '<STR_LIT>' , [ q [ <NUM_LIT:2> ] , B . y ] ) <EOL> C . set_ang_vel ( N , u [ <NUM_LIT:0> ] * B . x + u [ <NUM_LIT:1> ] * B . y + u [ <NUM_LIT:2> ] * B . z ) <EOL> C . set_ang_acc ( N , C . ang_vel_in ( N ) . diff ( t , B ) <EOL> + cross ( B . ang_vel_in ( N ) , C . ang_vel_in ( N ) ) ) <EOL> P = Point ( '<STR_LIT:P>' ) <EOL> O = P . locatenew ( '<STR_LIT:O>' , - r * B . z ) <EOL> P . set_vel ( N , ua [ <NUM_LIT:0> ] * A . x + ua [ <NUM_LIT:1> ] * A . y + ua [ <NUM_LIT:2> ] * A . z ) <EOL> O . v2pt_theory ( P , N , C ) <EOL> O . set_acc ( N , O . vel ( N ) . subs ( ua_zero ) . diff ( t , B ) <EOL> + cross ( B . ang_vel_in ( N ) , O . vel ( N ) . subs ( ua_zero ) ) ) <EOL> w_c_n_qd = qd [ <NUM_LIT:0> ] * A . z + qd [ <NUM_LIT:1> ] * B . x + qd [ <NUM_LIT:2> ] * B . y <EOL> kindiffs = Matrix ( [ dot ( w_c_n_qd - C . ang_vel_in ( N ) , uv ) for uv in B ] ) <EOL> qd_kd = solve ( kindiffs , qd ) <EOL> mprint ( kindiffs ) <EOL> steady_conditions = solve ( kindiffs . subs ( { qd [ <NUM_LIT:1> ] : <NUM_LIT:0> } ) , u ) <EOL> steady_conditions . update ( { qd [ <NUM_LIT:1> ] : <NUM_LIT:0> } ) <EOL> print ( steady_conditions ) <EOL> partial_w_C = [ C . ang_vel_in ( N ) . diff ( ui , N ) for ui in u + ua ] <EOL> partial_v_O = [ O . vel ( N ) . diff ( ui , N ) for ui in u + ua ] <EOL> partial_v_P = [ P . vel ( N ) . diff ( ui , N ) for ui in u + ua ] <EOL> print ( partial_w_C ) <EOL> print ( partial_v_O ) <EOL> print ( partial_v_P ) <EOL> F_O = m * g * A . z <EOL> F_P = Fx * A . x + Fy * A . y + Fz * A . z <EOL> Fr = [ dot ( F_O , pv_o ) + dot ( F_P , pv_p ) for pv_o , pv_p in <EOL> zip ( partial_v_O , partial_v_P ) ] <EOL> R_star_O = - m * O . acc ( N ) <EOL> I_C_O = inertia ( B , I , J , I ) <EOL> T_star_C = - ( dot ( I_C_O , C . ang_acc_in ( N ) ) + cross ( C . ang_vel_in ( N ) , dot ( I_C_O , C . ang_vel_in ( N ) ) ) ) <EOL> Fr_star = [ dot ( R_star_O , pv ) + dot ( T_star_C , pav ) for pv , pav in <EOL> zip ( partial_v_O , partial_w_C ) ] <EOL> Fr_star_steady = [ Fr_star_i . subs ( ud_zero ) . subs ( steady_conditions ) . expand ( ) <EOL> for Fr_star_i in Fr_star ] <EOL> mprint ( Fr ) <EOL> mprint ( Fr_star_steady ) <EOL> steady_turning_dynamic_equation = Fr [ <NUM_LIT:0> ] + Fr_star_steady [ <NUM_LIT:0> ] <EOL> p = Poly ( steady_turning_dynamic_equation , qd [ <NUM_LIT:0> ] ) <EOL> a , b , c = p . coeffs ( ) <EOL> discriminant = b * b - <NUM_LIT:4> * a * c </s>
