#-----[            ]-----
#-----[ Example 16 ]-----
#-----[            ]-----

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

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

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

#-----[ ... 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=1/E/I*(sp.Rational(1,24)*q*x1**4+sp.Rational(1,6)*c1*x1**3+sp.Rational(1,2)*c2*x1**2+c3*x1+c4)
w2=1/E/I*(sp.Rational(1,24)*q*x2**4+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(q/E/I))
print("\nw2=")
sp.pprint(sp.expand(w2_).collect(wB).collect(q/E/I))
print("\nw1'=")
sp.pprint(sp.expand(w1_.diff(x1)).collect(wB).collect(q/E/I))
print("\nw2'=")
sp.pprint(sp.expand(w2_.diff(x2)).collect(wB).collect(q/E/I))
print("\nw1''=")
sp.pprint(sp.expand(w1_.diff(x1,2)).collect(wB).collect(q/E/I))
print("\nw2''=")
sp.pprint(sp.expand(w2_.diff(x2,2)).collect(wB).collect(q/E/I))

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

eqn=sp.integrate(q*w1_.diff(wB),(x1,0,a))+sp.integrate(q*w2_.diff(wB),(x2,0,a))-k*wB \
-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))

print("\nUnit-displcament-method algebraic equations:")
sp.pprint(sp.expand(eqn).collect(wB).collect(q/E/I))

#-----[ ... 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);
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=[w1_.diff(x1,wB).subs({a:1,x1:ii}).evalf() for ii in x]
y4=[w2_.diff(x2,wB).subs({a:1,x2:ii}).evalf() for ii in x]
y5=[w1_.diff(x1,x1,wB).subs({a:1,x1:ii}).evalf() for ii in x]
y6=[w2_.diff(x2,x2,wB).subs({a:1,x2:ii}).evalf() for ii in x]

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

fig,([ax1,ax2])=plt.subplots(nrows=1,ncols=2,sharey=True)

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

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

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

ax1.legend()
ax2.legend()

plt.savefig('example16.png')

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