#-----[            ]-----
#-----[ Example 18 ]-----
#-----[            ]-----

#-----[ 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,S,F=sp.symbols('E I S F')
a,wB=sp.symbols('a w_B')
x1,x2,x3=sp.symbols('x1 x2 x3')

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

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

#-----[ ... 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
u3=wB*x3/2/a

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]})
u3_=u3

#-----[ 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))
print("\nu3=")
sp.pprint(sp.expand(u3_).collect(wB))
print("\nw1'=")
sp.pprint(sp.expand(w1_.diff(x1)).collect(wB).collect(F))
print("\nw2'=")
sp.pprint(sp.expand(w1_.diff(x1)).collect(wB).collect(F))
print("\nu3'=")
sp.pprint(sp.expand(u3_.diff(x3)).collect(wB))
print("\nw1''=")
sp.pprint(sp.expand(w1_.diff(x1,2)).collect(wB).collect(F))
print("\nw2''=")
sp.pprint(sp.expand(w2_.diff(x2,2)).collect(wB).collect(F))
print("\nu3''=")
sp.pprint(sp.expand(u3_.diff(x3,2)).collect(wB))

#-----[ system of elgebraic equations for ]-----
#-----[ Unit-displacement-method ...      ]-----

eqn=-F*w2_.diff(wB).subs(x2,a) \
+E*I*sp.integrate(w1_.diff(x1,2)*w1_.diff(x1,x1,wB),(x1,0,a)) \
+E*I*sp.integrate(w2_.diff(x2,2)*w2_.diff(x2,x2,wB),(x2,0,a)) \
+E*S*sp.integrate(u3_.diff(x3)*u3_.diff(x3,wB),(x3,0,2*a))

print("\nUnit-displcament-method algebraic equations:")
sp.pprint(sp.expand(eqn))

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

sol=sp.solve(eqn,wB)

print("wB=")
sp.pprint(sol[0])

import matplotlib.pyplot as plt
import numpy as np

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

x=np.linspace(0,1,100);
xx=np.linspace(0,2,100);
y1=[w1_.diff(wB).subs({a:1,x1:ii}).evalf() for ii in x]
y2=[w2_.diff(wB).subs({a:1,x2:ii}).evalf() for ii in x]
y3=[u3_.diff(wB).subs({a:1,x3:ii}).evalf() for ii in xx]
y4=[w1_.diff(x1,wB).subs({a:1,x1:ii}).evalf() for ii in x]
y5=[w2_.diff(x2,wB).subs({a:1,x2:ii}).evalf() for ii in x]
y6=[u3_.diff(x3,wB).subs({a:1,x3:ii}).evalf() for ii in xx]
y7=[w1_.diff(x1,x1,wB).subs({a:1,x1:ii}).evalf() for ii in x]
y8=[w2_.diff(x2,x2,wB).subs({a:1,x2:ii}).evalf() for ii in x]
y9=[u3_.diff(x3,x3,wB).subs({a:1,x3:ii}).evalf() for ii in xx]

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

fig,([ax1,ax2,ax3])=plt.subplots(nrows=1,ncols=3,sharey=True,figsize=(18,6))

ax1.set_title("$w_1(x)$")
ax2.set_title("$w_2(x)$")
ax3.set_title("$w_3(x)$")

ax1.plot(x,y1,label='$\phi_{w_B}$')
ax1.plot(x,y4,label='$\phi_{w_B}^\prime$')
ax1.plot(x,y7,label='$\phi_{w_B}^{\prime\prime}$')
ax2.plot(x,y2,label='$\phi_{w_B}$')
ax2.plot(x,y5,label='$\phi_{w_B}^\prime$')
ax2.plot(x,y8,label='$\phi_{w_B}^{\prime\prime}$')
ax3.plot(xx,y3,label='$\phi_{w_B}$')
ax3.plot(xx,y6,label='$\phi_{w_B}^\prime$')
ax3.plot(xx,y9,label='$\phi_{w_B}^{\prime\prime}$')

ax1.set_xlabel("$x/a$")
ax2.set_xlabel("$x/a$")
ax3.set_xlabel("$x/a$")

ax1.legend()
ax2.legend()
ax3.legend()

plt.savefig('example18.png')

print("\ngraphs of fundamental functions appearing in dummy-displacemnts functions")
print("were plotted in the figgure >example18.png<")