from pint import UnitRegistry
import sympy
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import sys
%matplotlib inline
from IPython.display import display

"""
Explanation: Semi-Monocoque Theory: corrective solutions
End of explanation
"""


from Section import Section

"""
Explanation: Import Section class, which contains all calculations
End of explanation
"""


ureg = UnitRegistry()
sympy.init_printing()

"""
Explanation: Initialization of sympy symbolic tool and pint for dimension analysis (not really implemented rn as not directly compatible with sympy)
End of explanation
"""


A, A0, t, t0, a, b, h, L, E, G = sympy.symbols('A A_0 t t_0 a b h L E G', positive=True)

"""
Explanation: Define sympy parameters used for geometric description of sections
End of explanation
"""


values = [(A, 150 * ureg.millimeter**2),(A0, 250  * ureg.millimeter**2),(a, 400 * ureg.millimeter), \
          (b, 20 * ureg.millimeter),(h, 150 * ureg.millimeter),(L, 2000 * ureg.millimeter), \
          (t, 0.8 *ureg.millimeter),(E, 72e3 * ureg.MPa), (G, 28e3 * ureg.MPa)]
datav = [(v[0],v[1].magnitude) for v in values]

"""
Explanation: We also define numerical values for each symbol in order to plot scaled section and perform calculations
End of explanation
"""


stringers = {1:[(a,h),A],
             2:[(sympy.Rational(1,2)*a,h),A],
             3:[(sympy.Integer(0),h),A],
             4:[(sympy.Integer(0),sympy.Integer(0)),A],
             5:[(sympy.Rational(1,2)*a,sympy.Integer(0)),A],
             6:[(a,sympy.Integer(0)),A]}

panels = {(1,2):t,
          (3,4):t,
          (4,5):t,
          (5,6):t,
          (6,1):t,
          (5,2):t}

"""
Explanation: First example: Simple rectangular symmetric section
Define graph describing the section:
1) stringers are nodes with parameters:
- x coordinate
- y coordinate
- Area
2) panels are oriented edges with parameters:
- thickness
- lenght which is automatically calculated
End of explanation
"""


S1 = Section(stringers, panels)

S1.cycles

"""
Explanation: Define section and perform first calculations
End of explanation
"""


start_pos={ii: [float(S1.g.node[ii]['ip'][i].subs(datav)) for i in range(2)] for ii in S1.g.nodes() }

plt.figure(figsize=(12,8),dpi=300)
nx.draw(S1.g,with_labels=True, arrows= True, pos=start_pos)
plt.arrow(0,0,20,0)
plt.arrow(0,0,0,20)
#plt.text(0,0, 'CG', fontsize=24)
plt.axis('equal')
plt.title("Section in starting reference Frame",fontsize=16);

"""
Explanation: Plot of S1 section in original reference frame
Define a dictionary of coordinates used by Networkx to plot section as a Directed graph.
Note that arrows are actually just thicker stubs
End of explanation
"""


positions={ii: [float(S1.g.node[ii]['pos'][i].subs(datav)) for i in range(2)] for ii in S1.g.nodes() }

x_ct, y_ct = S1.ct.subs(datav)

plt.figure(figsize=(12,8),dpi=300)
nx.draw(S1.g,with_labels=True, pos=positions)
plt.plot([0],[0],'o',ms=12,label='CG')
plt.plot([x_ct],[y_ct],'^',ms=12, label='SC')
#plt.text(0,0, 'CG', fontsize=24)
#plt.text(x_ct,y_ct, 'SC', fontsize=24)
plt.legend(loc='lower right', shadow=True)
plt.axis('equal')
plt.title("Section in pricipal reference Frame",fontsize=16);

"""
Explanation: Plot of S1 section in inertial reference Frame
Section is plotted wrt center of gravity and rotated (if necessary) so that x and y are principal axes.
Center of Gravity and Shear Center are drawn
End of explanation
"""


Tx, Ty, Nz, Mx, My, Mz, F, ry, ry, mz = sympy.symbols('T_x T_y N_z M_x M_y M_z F r_y r_x m_z')

S1.set_loads(_Tx=0, _Ty=0, _Nz=0, _Mx=0, _My=0, _Mz=Mz)

"""
Explanation: Standard Solution
End of explanation
"""


S1.compute_stringer_actions()
S1.compute_panel_fluxes();

S1.N

S1.q

S1.set_loads(_Tx=0, _Ty=Ty, _Nz=0, _Mx=Mx, _My=0, _Mz=0)

S1.compute_stringer_actions()
S1.compute_panel_fluxes();

S1.N

S1.q

"""
Explanation: Loads
End of explanation
"""


S1.compute_L()

S1.L

"""
Explanation: Autosolutions
Compute L matrix: with 6 nodes we expect 3 dofs, one with symmetric load and one with antisymmetric load
End of explanation
"""


S1.compute_H()

S1.H

"""
Explanation: Compute H matrix
End of explanation
"""


S1.compute_KM(A,a/2,t,datav)

S1.Ktilde

S1.Mtilde

"""
Explanation: Compute $\tilde{K}$ and $\tilde{M}$ as:
$$\tilde{K} = L^T \cdot \left[ \frac{A}{A_0} \right] \cdot L$$
$$\tilde{M} = H^T \cdot \left[ \frac{l}{l_0}\frac{t_0}{t} \right] \cdot L$$
End of explanation
"""


sol_data = (sympy.N(S1.Ktilde.subs(datav).inv())*(sympy.N(S1.Mtilde.subs(datav)))).eigenvects()

"""
Explanation: Compute eigenvalues and eigenvectors as:
$$\left| \mathbf{I} \cdot \beta^2 - \mathbf{\tilde{K}}^{-1} \cdot \mathbf{\tilde{M}} \right| = 0$$
We substitute some numerical values to simplify the expressions
End of explanation
"""


β2 = [sol[0] for sol in sol_data]
β2

"""
Explanation: Eigenvalues correspond to $\beta^2$
End of explanation
"""


X = [sol[2][0] for sol in sol_data]
X

"""
Explanation: Eigenvectors are orthogonal as expected
End of explanation
"""


λ = [sympy.N(sympy.sqrt(E*A*b/(G*t)*βi).subs(datav)) for βi in β2]
λ

"""
Explanation: From $\beta_i^2$ we compute:
$$\lambda_i = \sqrt{\frac{E A_0 l_0}{G t_0} \beta_i^2}$$
substuting numerical values
End of explanation
"""


