import bikewheelcalc as bc
import matplotlib.pyplot as plt


# Create an example wheel
wheel = bc.BicycleWheel()
wheel.hub = bc.Hub(diam1=0.04, width1=0.03)
wheel.rim = bc.Rim.general(radius=0.3,
                           area=82.0e-6,
                           I11=5620e-12,
                           I22=1187e-12,
                           I33=1124e-12,
                           Iw=0.0,
                           young_mod=69.0e9,
                           shear_mod=26.0e9)

wheel.lace_cross(n_spokes=36, n_cross=3, diameter=2.0e-3,
                 young_mod=210e9, offset=0.0)

# Create finite-element model from wheel
fem = bc.BicycleWheelFEM(wheel, verbose=True)

# Create a rigid body for the hub nodes and constrain it
R1 = bc.RigidBody('hub', [0, 0, 0], fem.get_hub_nodes())
fem.add_rigid_body(R1)
fem.add_constraint(R1.node_id, range(6))

# Add an upward force, distributed over bottom 3 nodes
fem.add_force(0, 1, 500)
fem.add_force(1, 1, 250)
fem.add_force(35, 1, 250)

soln = fem.solve(pretension=1000)

# Draw deformed wheel
f1 = plt.figure(1)
f_def = soln.plot_deformed_wheel(rel_scale=0.1)
plt.axis('off')

f2 = plt.figure(2)
soln.plot_spoke_tension(fig=f2)
f2.gca().set_yticklabels([])

plt.show()

from bikewheelcalc import BicycleWheel, Rim, Hub, ModeMatrix
import matplotlib.pyplot as plt
import numpy as np


# Create an example wheel and rim
wheel = BicycleWheel()
wheel.hub = Hub(width=0.05, diameter=0.05)
wheel.rim = Rim(radius=0.3, area=100e-6,
                I_lat=200./69e9, I_rad=100./69e9, J_tor=25./26e9, I_warp=0.0,
                young_mod=69e9, shear_mod=26e9)
wheel.lace_cross(n_spokes=36, n_cross=3, diameter=2.0e-3, young_mod=210e9)


# Create a ModeMatrix model with 24 modes
mm = ModeMatrix(wheel, N=24)

# Create a 500 Newton pointing radially inwards at theta=0
F_ext = mm.F_ext(0., np.array([0., 500., 0., 0.]))

# Calculate stiffness matrix
K = mm.K_rim(tension=False) + mm.K_spk(smeared_spokes=True, tension=False)

# Solve for the mode coefficients
dm = np.linalg.solve(K, F_ext)

# Calculate the rotational stiffness
# Get radial deflection
theta = np.linspace(-np.pi, np.pi, 100)
rad_def = mm.rim_def_rad(theta, dm)


# Draw deformed wheel
f1 = plt.figure(1)
plt.plot(theta, 1000.*rad_def)

plt.xlabel('theta')
plt.ylabel('Radial deflection [mm]')

# f2 = plt.figure(2)
# soln.plot_spoke_tension(fig=f2)
# f2.gca().set_yticklabels([])

plt.tight_layout()
plt.show()

Update radial force example to use ModeMatrix