# 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 boundary_conditions():
#     def __init__(self, concentrated_load=np.array([[0, 0, 0], ]), load_x=[1],
#                        distributed_load = 0):
#         self.concentrated_load = concentrated_load
#         self.concentrated_x = load_x
#         self.concentrated_n = len(concentrated_load)
#         self.distributed_load = None
# 
#         self.distributed_load = distributed_load
# 
#         if len(concentrated_load) != len(load_x):
#             raise Exception('load values and x lists have to match')


from aeropy.geometry.parametric import CoordinateSystem
from aeropy.structural.stable_solution import (properties, boundary_conditions)
from aeropy.structural.shell import shell
from aeropy.xfoil_module import output_reader
from scipy.integrate import quad
from scipy.optimize import minimize
from optimization_tools.DOE import DOE
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pickle

def find_deflection(y, l, P):
    def _to_integrate(y):
        num = P/E/I*(l*y-y**2/2)
        den = np.sqrt(1-P**2/E**2/I**2*(l*y-y**2/2)**2)
        if np.isnan(den):
            return 100
        else:
            return num/den
    
    x = []
    for y_i in y:
        x_i = current_L = quad(_to_integrate, 0, y_i)[0]
        x.append(x_i)
    return x
    
np.set_printoptions(precision=5)

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

# Beam properties
bp = properties()
bc = boundary_conditions(concentrated_load=np.array([[0.0, 0.0, -1.0], ]))
EB_solution = bc.concentrated_load[0][2]/(6*bp.young*bp.inertia) * \
    np.array([0, 0, 3, -1])

curve_parent = CoordinateSystem.polynomial(D=[0, 0, 0, 0], chord = 1, color = 'b')
curve_child = CoordinateSystem.polynomial(D=EB_solution, chord = 1, color  ='0.5')


# 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 boundary_conditions():
#     def __init__(self, concentrated_load=np.array([[0, 0, 0], ]), load_x=[1],
#                        distributed_load = 0):
#         self.concentrated_load = concentrated_load
#         self.concentrated_x = load_x
#         self.concentrated_n = len(concentrated_load)
#         self.distributed_load = None
# 
#         self.distributed_load = distributed_load
# 
#         if len(concentrated_load) != len(load_x):
#             raise Exception('load values and x lists have to match')


# Path: examples/structural/beam_theory/beam.py
from aeropy.geometry.parametric import CoordinateSystem
from aeropy.structural.stable_solution import (properties, boundary_conditions)
from aeropy.structural.shell import shell
from aeropy.xfoil_module import output_reader
from scipy.integrate import quad
from scipy.optimize import minimize
from optimization_tools.DOE import DOE
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pickle


def find_deflection(y, l, P):
    def _to_integrate(y):
        num = P/E/I*(l*y-y**2/2)
        den = np.sqrt(1-P**2/E**2/I**2*(l*y-y**2/2)**2)
        if np.isnan(den):
            return 100
        else:
            return num/den
    
    x = []
    for y_i in y:
        x_i = current_L = quad(_to_integrate, 0, y_i)[0]
        x.append(x_i)
    return x
    
np.set_printoptions(precision=5)

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

# Beam properties
bp = properties()
bc = boundary_conditions(concentrated_load=np.array([[0.0, 0.0, -1.0], ]))
EB_solution = bc.concentrated_load[0][2]/(6*bp.young*bp.inertia) * \
    np.array([0, 0, 3, -1])

curve_parent = CoordinateSystem.polynomial(D=[0, 0, 0, 0], chord = 1, color = 'b')
curve_child = CoordinateSystem.polynomial(D=EB_solution, chord = 1, color  ='0.5')


beam = shell(curve_parent, curve_child, bp, bc)