sf.sagews

from sympy.physics.units import *
from sympy import *
from numpy import linspace, array
from matplotlib import pyplot as plt
from sys import exit

print("\n--- User input: ---------------------")

symbolic = True
symbolic = False

if symbolic:
    # quantity = symbol:
    a, q = var("a, q")
else:
    # Derived units:
    N = kg*m/s/s
    # factors:
    a_val, q_val = 1, 1
    # a_val, q_val = 1, -2
    #
    # quantity =  factor times unit:
    a, q = a_val * m, q_val * N/m

print("\n--- a: -------------------------------")
x = var("x")

# Resultant:
hx, Rx = x/2, q*x

# Section forces:
Q, M = - Rx, - hx*Rx

print("\n--- b: -------------------------------")

print("\n--- c: -------------------------------")

# make unitless, but keep quantities:
xia = linspace(float(0), float(1), 40)
xia = [xiai *m for xiai in xia]

xi = var("xi")

if symbolic:
    Q = Q / (q*a)
    M = M / (q*a*a)
else:
    Q = Q / N
    M = M / (N*m)

[Q,M] = [X.expand().subs(x/a,xi) for X in [Q,M]]
# Make array for plotting:
xia = linspace(float(0),float(1),float(20))
def fill_array(f,x,xia):
    if f.is_constant(x):
        res = full_like(xia,f)
    else:
        res = lambdify(x,f,"numpy")
        res = array(res(xia))
    return res

[Q,M] = [fill_array(X,xi,xia) for X in [Q,M]]

if symbolic:
    tmp = plt.plot(xia, array(Q), "b-",  label= "$Q /  qa$")
    tmp = plt.plot(xia, array(M), "r--", label= "$M /  qa^2$")
    tmp = plt.xlabel(r"$x / a$")
else:
    tmp = plt.plot(xia, array(Q), "b-",  label= "$Q /  \mathrm{N}$")
    tmp = plt.plot(xia, array(M), "r--", label= "$M / \mathrm{Nm}$")
    tmp = plt.xlabel(r"$x / \mathrm{m}$")

# annotation:
plt.grid()
tmp = plt.legend(loc="best")

# save, show:
plt.savefig("plot_temp.svg", transparent=True)
plt.show()

--- User input: --------------------- --- a: ------------------------------- --- b: ------------------------------- --- c: -------------------------------