so im looking at this code, intended to look at a given force at the contact patch of a tire, and derive the forces within the a-arms in a double wishbone suspension. i would like you to look THOROUGHLY over the entire code, as well as add functionality for viewing the a-arm geometry in 3D, and upper ball joint forces. hope the code makes sense.

import numpy as np 
import math
from matplotlib import pyplot as plt
from mpl_toolkits import mplot3d
from matplotlib import collections  as mc
from mpl_toolkits.mplot3d import axes3d
from numpy.linalg import inv
g = 9.81

# car values

# mass
mass = 310

# wheelbase
wb = 1550

# track width
tw = 1200

# center of mass height
cm = 280

# static weight distribution (front load)
wd = 0.5

# input values front (ONLY RUN THIS BLOCK OR NEXT BLOCK FOR F/R FORCES)

# list of coordinates

# f = fore, a = aft, t = top, b = bottom, d = damper, tr = tie rod, cP = contact patch
# coordinates ([upright side], [chassis side])

# A
fb = np.array([[3.7628,-37.0299,108.4733],[208.68,-386.8966,124.4613]])

# B
ab = np.array([[3.7628,-37.0299,108.4733],[-150,-386.8966,155.2616]])


# C
ft = np.array([[-3.7628,-49.3078,297.9267],[106.6266,-341.9893,283.17]])


# D
at = np.array([[-3.7628,-49.3078,297.9267],[-125,-341.9893,263.0367]])


# E
d = np.array([[-4.5,-86.1697,264.1451],[-4.5,-270.4791,95.2367]])


# F
tr = np.array([[52.3129,-39.3953,144.9733], [81.6177,-417.12,183.0733]])


cP = np.array([0,0,0])

# force applied *in G's (long, lat, normal)
fap = np.array([0,1.7,0]) 


# calculation of forces w/ weight distribution

wtLat = mass*fap[1]*g*cm/tw
wtLong = mass*fap[0]*g*cm/tw
fap[2] = wd*mass*9.81/2+(wtLat+wtLong)/2
fzPercent = fap[2]/(mass*g)
fap[1] = fap[1]*g*mass*fzPercent
fap[0] = fap[0]*g*mass*fzPercent

# In[4]:


# auxillary values

# force application coords
A = fb[0]
B = fb[0]
C = ft[0]
D = at[0]
E = d[0]
F = tr[0]

# moment from contact patch
mom = np.cross(cP,fap)

# total applied force/moment vector
fapTot = np.resize(np.array([fap,mom]),(6,1))

# delta x for each tube
dxA = fb[1][0]-fb[0][0]
dxB = ab[1][0]-ab[0][0]
dxC = ft[1][0]-ft[0][0]
dxD = at[1][0]-at[0][0]
dxE = d[1][0]-d[0][0]
dxF = tr[1][0]-tr[0][0]

# delta y for each tube
dyA = fb[1][1]-fb[0][1]
dyB = ab[1][1]-ab[0][1]
dyC = ft[1][1]-ft[0][1]
dyD = at[1][1]-at[0][1]
dyE = d[1][1]-d[0][1]
dyF = tr[1][1]-tr[0][1]

# delta z for each tube
dzA = fb[1][2]-fb[0][2]
dzB = ab[1][2]-ab[0][2]
dzC = ft[1][2]-ft[0][2]
dzD = at[1][2]-at[0][2]
dzE = d[1][2]-d[0][2]
dzF = tr[1][2]-tr[0][2]

# length of each tube
lA = np.sqrt(dxA**2+dyA**2+dzA**2)
lB = np.sqrt(dxB**2+dyB**2+dzB**2)
lC = np.sqrt(dxC**2+dyC**2+dzC**2)
lD = np.sqrt(dxD**2+dyD**2+dzD**2)
lE = np.sqrt(dxE**2+dyE**2+dzE**2)
lF = np.sqrt(dxF**2+dyF**2+dzF**2)

# Force vector for each a-arm
forceA = np.array([dxA/lA, dyA/lA, dzA/lA])
forceB = np.array([dxB/lB, dyB/lB, dzB/lB])
forceC = np.array([dxC/lC, dyC/lC, dzC/lC])
forceD = np.array([dxD/lD, dyD/lD, dzD/lD])
forceE = np.array([dxE/lE, dyE/lE, dzE/lE])
forceF = np.array([dxF/lF, dyF/lF, dzF/lF])

# Moment arms w.r.t origin
mA = np.cross(A, forceA)
mB = np.cross(B, forceB)
mC = np.cross(C, forceC)
mD = np.cross(D, forceD)
mE = np.cross(E, forceE)
mF = np.cross(F, forceF)

# Calculation of instant centers
p1 = fb[1][2]
p2 = ft[1][2]

m1 = dzA/dyA
m2 = dzC/dyC

t = (p1-p2)/(m2-m1)

ic = np.array([fb[1][0],fb[1][1]+t,fb[1][2]+m2*t])
ric = (ic-np.array([0,600,0]))*np.array([0,-1,1])

# Tie rod vector
vec = tr[0]-ric
print(vec)

# plot elements

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
t = np.linspace(0,1,100)
X = fb[0][0] + t*dxA
Y = fb[0][1] + t*dyA
Z = fb[0][2] + t*dzA
ax.plot3D(X, Y, Z, 'green')
X = ab[0][0] + t*dxB
Y = ab[0][1] + t*dyB
Z = ab[0][2] + t*dzB
ax.plot3D(X, Y, Z, 'green')
X = at[0][0] + t*dxC
Y = at[0][1] + t*dyC
Z = at[0][2] + t*dzC
ax.plot3D(X, Y, Z, 'blue')
X = ft[0][0] + t*dxD
Y = ft[0][1] + t*dyD
Z = ft[0][2] + t*dzD
ax.plot3D(X, Y, Z, 'blue')
X = d[0][0] + t*(d[1][0]-d[0][0])
Y = d[0][1] + t*(d[1][1]-d[0][1])
Z = d[0][2] + t*(d[1][2]-d[0][2])
ax.plot3D(X, Y, Z, 'red')
X = tr[0][0] + t*(tr[1][0]-tr[0][0])
Y = tr[0][1] + t*(tr[1][1]-tr[0][1])
Z = tr[0][2] + t*(tr[1][2]-tr[0][2])
ax.plot3D(X, Y, Z, 'purple')

ax.set_xlabel('$X$')
ax.set_ylabel('$Y$')
ax.set_zlabel('$Z$')

# concatenation of forces into matrix

forceArray = np.array([[dxA/lA],[dxB/lB],[dxC/lC],[dxD/lD],[dxE/lE],[dxF/lF],                         [dyA/lA],[dyB/lB],[dyC/lC],[dyD/lD],[dyE/lE],[dyF/lF],                         [dzA/lA],[dzB/lB],[dzC/lC],[dzD/lD],[dzE/lE],[dzF/lF],                         [mA[0]], [mB[0]], [mC[0]], [mD[0]], [mE[0]], [mF[0]],                         [mA[1]], [mB[1]], [mC[1]], [mD[1]], [mE[1]], [mF[1]],                         [mA[2]], [mB[2]], [mC[2]], [mD[2]], [mE[2]], [mF[2]]])
forceMatrix = np.reshape(forceArray, (6,6))

# solving for forces in elements

forces = np.dot(inv(forceMatrix), fapTot)
print(forces)

LBJ = np.array([[forces[0][0]*dxA/lA+forces[1][0]*dxB/lB+forces[4][0]*dxE/lE],                [forces[0][0]*dyA/lA+forces[1][0]*dyB/lB+forces[4][0]*dyE/lE],                [forces[0][0]*dzA/lA+forces[1][0]*dzB/lB+forces[4][0]*dzE/lE]])
print(LBJ)

another issue i am having, is that the forces here seem consistently higher than expected, based on other research papers and efforts. Why?