# Path: aeropy/structural/beam.py
# class beam_chen():
#     def __init__(self, geometry, properties, load, s):
#         self.g = geometry
#         self.p = properties
#         self.l = load
#         self.s = s
#         self.length = self.g.arclength()
# 
#     def G(self, x):
#         c = 1/self.p.young/self.p.inertia
#         # print(x)
#         # return c*(quad(self.M, 0, x)[0])
#         c = self.g.chord
#         return self.l.concentrated_load[0][-1]/self.p.young/self.p.inertia*(c*x-x**2/2)
# 
#     def M(self, x):
#         c = self.g.chord
#         for i in range(len(self.l.concentrated_s)):
#             concentrated_x_i = self._x(self.l.concentrated_s[i])
#             return self.l.concentrated_load[0][-1]*(concentrated_x_i*-x)
# 
#     def s_to_x(self):
#         self.x = np.zeros(len(self.s))
#         for i in range(len(self.s)):
#             self.x[i] = self._x(self.s[i])
# 
#     def _x(self, s):
#         def _to_minimize(l):
#             def _to_integrate(x):
#                 den = np.sqrt(1-self.G(x)**2)
#                 if np.isnan(den):
#                     return 100
#                 else:
#                     return 1/den
#             l = l[0]
#             current_L = quad(_to_integrate, 0, l)[0]
#             return abs(s-current_L)
#         return minimize(_to_minimize, s, method = 'Nelder-Mead',).x[0]
# 
#     def find_deflection(self):
#         def _to_integrate(x):
#             G = self.G(x)
#             den = np.sqrt(1-G**2)
#             if np.isnan(den):
#                 return 100
#             else:
#                 return G/den
# 
#         self.y = []
#         for x_i in self.x:
#             y_i = current_L = quad(_to_integrate, 0, x_i)[0]
#             self.y.append(y_i)
# 
# Path: aeropy/structural/stable_solution.py
# class properties():
#     def __init__(self, young=70e9, poisson=.3, dimensions=[0.01, 0.01],
#                  crosssection='square'):
#         """Stores calculates material or structural properties. If 'square',
#         dimensions = [width, height]"""
#         self.young = young
#         self.poisson = poisson
#         self.dimensions = dimensions
#         self.length = 1
# 
#         if crosssection == 'square':
#             self.area = dimensions[0]*dimensions[1]
#             self.inertia = (dimensions[0]*dimensions[1]**3)/12.
#         if young < 0 or poisson < 0 or poisson > 1 or self.area <= 0:
#             raise Exception('Material properties need to make sense')
# 
# class loads():
#     def __init__(self, concentrated_load=np.array([[0, -1], ]), load_s=[1],
#                        distributed_function = 0):
#         self.concentrated_load = concentrated_load
#         self.concentrated_s = load_s
#         self.distributed_function = distributed_function


import matplotlib.pyplot as plt
import numpy as np
import pickle
from aeropy.structural.beam import beam_chen
from aeropy.structural.stable_solution import properties, loads
from aeropy.geometry.parametric import CoordinateSystem




g = CoordinateSystem.polynomial(D=[0, 0, 0, 0], chord = 1, color = 'b')
p = properties()
l = loads()
s = np.linspace(0, 1, 10)

b = beam_chen(g, p, l, s)
b.s_to_x()
b.find_deflection()

# Results from Abaqus
abaqus_data = pickle.load(open('neutral_line.p', 'rb'))

# Euler-Bernoulle
EB_solution = l.concentrated_load[0][-1]/(6*p.young*p.inertia) * \
    np.array([0, 0, 3, -1])

curve_EB = CoordinateSystem.polynomial(D=EB_solution, chord = 1, color  ='0.5')

curve_EB.calculate_x1(b.s)

# Path: aeropy/structural/beam.py
# class beam_chen():
#     def __init__(self, geometry, properties, load, s):
#         self.g = geometry
#         self.p = properties
#         self.l = load
#         self.s = s
#         self.length = self.g.arclength()
# 
#     def G(self, x):
#         c = 1/self.p.young/self.p.inertia
#         # print(x)
#         # return c*(quad(self.M, 0, x)[0])
#         c = self.g.chord
#         return self.l.concentrated_load[0][-1]/self.p.young/self.p.inertia*(c*x-x**2/2)
# 
#     def M(self, x):
#         c = self.g.chord
#         for i in range(len(self.l.concentrated_s)):
#             concentrated_x_i = self._x(self.l.concentrated_s[i])
#             return self.l.concentrated_load[0][-1]*(concentrated_x_i*-x)
# 
#     def s_to_x(self):
#         self.x = np.zeros(len(self.s))
#         for i in range(len(self.s)):
#             self.x[i] = self._x(self.s[i])
# 
#     def _x(self, s):
#         def _to_minimize(l):
#             def _to_integrate(x):
#                 den = np.sqrt(1-self.G(x)**2)
#                 if np.isnan(den):
#                     return 100
#                 else:
#                     return 1/den
#             l = l[0]
#             current_L = quad(_to_integrate, 0, l)[0]
#             return abs(s-current_L)
#         return minimize(_to_minimize, s, method = 'Nelder-Mead',).x[0]
# 
#     def find_deflection(self):
#         def _to_integrate(x):
#             G = self.G(x)
#             den = np.sqrt(1-G**2)
#             if np.isnan(den):
#                 return 100
#             else:
#                 return G/den
# 
#         self.y = []
#         for x_i in self.x:
#             y_i = current_L = quad(_to_integrate, 0, x_i)[0]
#             self.y.append(y_i)
# 
# Path: aeropy/structural/stable_solution.py
# class properties():
#     def __init__(self, young=70e9, poisson=.3, dimensions=[0.01, 0.01],
#                  crosssection='square'):
#         """Stores calculates material or structural properties. If 'square',
#         dimensions = [width, height]"""
#         self.young = young
#         self.poisson = poisson
#         self.dimensions = dimensions
#         self.length = 1
# 
#         if crosssection == 'square':
#             self.area = dimensions[0]*dimensions[1]
#             self.inertia = (dimensions[0]*dimensions[1]**3)/12.
#         if young < 0 or poisson < 0 or poisson > 1 or self.area <= 0:
#             raise Exception('Material properties need to make sense')
# 
# class loads():
#     def __init__(self, concentrated_load=np.array([[0, -1], ]), load_s=[1],
#                        distributed_function = 0):
#         self.concentrated_load = concentrated_load
#         self.concentrated_s = load_s
#         self.distributed_function = distributed_function


# Path: examples/structural/beam_theory/beam_chen.py
import matplotlib.pyplot as plt
import numpy as np
import pickle
from aeropy.structural.beam import beam_chen
from aeropy.structural.stable_solution import properties, loads
from aeropy.geometry.parametric import CoordinateSystem





g = CoordinateSystem.polynomial(D=[0, 0, 0, 0], chord = 1, color = 'b')
p = properties()
l = loads()
s = np.linspace(0, 1, 10)

b = beam_chen(g, p, l, s)
b.s_to_x()
b.find_deflection()

# Results from Abaqus
abaqus_data = pickle.load(open('neutral_line.p', 'rb'))

# Euler-Bernoulle
EB_solution = l.concentrated_load[0][-1]/(6*p.young*p.inertia) * \
    np.array([0, 0, 3, -1])

curve_EB = CoordinateSystem.polynomial(D=EB_solution, chord = 1, color  ='0.5')

curve_EB.calculate_x1(b.s)

eulerBernoulle = curve_EB.r(b.s)