#-----[            ]-----
#-----[ Example 20 ]-----
#-----[            ]-----

#-----[ import library >sympy<  ]-----
#-----[ and setup best printing ]-----

import sympy as sp
sp.init_printing()

#-----[ symbol declaration ]-----

c1,c2,c3,c4=sp.symbols('c1 c2 c3 c4')
d1,d2,d3,d4=sp.symbols('d1 d2 d3 d4')
E,I=sp.symbols('E I')
a,F=sp.symbols('a F')
x1,x2=sp.symbols('x1 x2')

#-----[ algebraic equations for constants ]-----
#-----[ of the equations of the elastic   ]-----
#-----[ curve ...                         ]-----

eqn1=c3
eqn2=c4
eqn3=sp.Rational(1,6)*c1*a**3+sp.Rational(1,2)*c2*a**2+c3*a+c4-d4
eqn4=sp.Rational(1,2)*c1*a**2+c2*a+c3-d3
eqn5=c1*a+c2-d2
eqn6=c1-d1+1/E/I-S/I/2/a*(sp.Rational(1,6)*c1*a**3+sp.Rational(1,2)*c2*a**2+c3*a+c4)
eqn7=d1*a+d2
eqn8=d1

#-----[ ... and their solution ]-----

sol=sp.solve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8],[c1,c2,c3,c4,d1,d2,d3,d4])

print("c1=");sp.pprint(sol[c1])
print("c2=");sp.pprint(sol[c2])
print("c3=");sp.pprint(sol[c3])
print("c4=");sp.pprint(sol[c4])
print("d1=");sp.pprint(sol[d1])
print("d2=");sp.pprint(sol[d2])
print("d3=");sp.pprint(sol[d3])
print("d4=");sp.pprint(sol[d4])

#-----[ elastic curve for section I and II ]-----

w1=sp.Rational(1,6)*c1*x1**3+sp.Rational(1,2)*c2*x1**2+c3*x1+c4
w2=sp.Rational(1,6)*d1*x2**3+sp.Rational(1,2)*d2*x2**2+d3*x2+d4

w1_=w1.subs({c1:sol[c1],c2:sol[c2],c3:sol[c3],c4:sol[c4]})
w2_=w2.subs({d1:sol[d1],d2:sol[d2],d3:sol[d3],d4:sol[d4]})

#-----[ elastic curves expressed using ]-----
#-----[ dummy-displacement-functions   ]-----

print("\nw1=")
sp.pprint(sp.expand(w1_).collect(wB).collect(F))
print("\nw2=")
sp.pprint(sp.expand(w2_).collect(wB).collect(F))

#-----[ Betti's theorem ]-----

wB=F*w2_.subs(x2,a)

print("wB=")
sp.pprint(sp.simplify(wB))

#-----[ Plot of funcdamental solution ]-----

import matplotlib.pyplot as plt
import numpy as np

#-----[ make data ]-----

x=np.linspace(0,1,100);
xx=np.linspace(1,2,100);

y1=[w1_.subs({a:1,S:-1,I:1,E:1,x1:ii}).evalf() for ii in x]
y2=[w2_.subs({a:1,S:-1,I:1,E:1,x2:ii}).evalf() for ii in x]

#-----[ plot ]-----

fig,ax=plt.subplots()

ax.set_title("Deflextion of beam under load $F_B=1$")
ax.invert_yaxis()

ax.plot(x,y1)
ax.plot(xx,y2)

ax.set_xlabel("$x/a$")
ax.set_ylabel("$w(5Fa^3)^{-1}E(Sa^2+6I)$")

plt.savefig('example20.png')

print("\ngraph of fundamental solution was made in figure >example20.png<")