SVV_Structures/integrals.py at master · ftagliaca/SVV_Structures · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
import numpy as np
from aileronProperties import Aileron
from aero_loads import AerodynamicLoad

A320 = Aileron(0.547, 2.771, 0.153, 1.281, 2.681, 28.0, 22.5, 1.1, 2.9, 1.2, 1.5, 2.0, 17, 1.103, 1.642, 26, 91.7)
AELoading = AerodynamicLoad(A320, "data/aerodynamicloada320.dat")

def MapAeroLoading(filename):
    return np.genfromtxt(filename, delimiter=",")

z,x,AeroLoading = AELoading.interpolate_predefined_grid()

def LinearInterpolatePos(Q1, Q2, x_0, x_1, x):
    #print("Oh come on!", x)
    return Q1 + (Q2-Q1)/(x_1 - x_0)*(x-x_0)

def integrate_1d(x, y, x_f):
    if abs(x_f) > abs(x[-1]):   #if x_f is outside of the range covered by input, we return the total integral of what we can integrate over
        x_f = x[-1]
    if abs(x_f) < abs(x[0]):    #if x_f is lower than the lowest x_value, return 0
        return 0
    if len(x) < 2:
        return 0
    total = 0
    i = 1
    while abs(x[i]) < abs(x_f):
        total += (x[i] - x[i-1])*(y[i]+y[i-1])/2
        i += 1
    i -= 1
    if abs(x_f) > abs(x[i]):
        total += (x_f - x[i])*(LinearInterpolatePos(y[i], y[i+1], x[i], x[i+1], x_f) + y[i])/2
    return total

def integrate_1d_list(x, y, x_f):
    int_list = []
    x_new = []
    i = 1
    if abs(x_f) > abs(x[-1]):
        x_f = x[-1]
    if abs(x_f) < abs(x[0]):
        return [0], [x_f]
    if len(x) < 2:
        return [0], [x_f]
    int_list.append(0)
    x_new.append(x[0])
    while abs(x[i]) < abs(x_f):
        int_list.append(integrate_1d(x, y, x[i]))
        x_new.append(x[i])
        i += 1
    i -= 1
    if abs(x_f) > abs(x[i]):
        int_list.append(integrate_1d(x, y, x_f))
        x_new.append(x_f)

    return int_list, x_new


def integrate_1d_tau(x, y, x_f, x_sc):
    #inputs: x; an array containing all the x-locationx of the points. y; an array containing all the y values of the points. x_i; the value of x to where we integrate.
    #outputs int_x; the total integrated value until x_i.

    #assert len(x) == len(y) #both arrays must have the same size

    if abs(x_f) > abs(x[-1]):   #if x_f is outside of the range covered by input, we return the total integral of what we can integrate over
        x_f = x[-1]
    if abs(x_f) < abs(x[0]):    #if x_f is lower than the lowest x_value, return 0
        return 0
    if len(x) < 2:
        return 0
    total = 0
    i = 1
    while abs(x[i]) < abs(x_f):
        total += (x[i] - x[i-1])*(y[i]+y[i-1])/2 * ((x[i] + x[i-1])*0.5 - x_sc)
        #print(((x[i] + x[i-1])*0.5 - x_sc))
        #print(x[i-1], x[i], total)
        i += 1
    i -= 1
    if abs(x_f) > abs(x[i]):
        total += (x_f - x[i])*(LinearInterpolatePos(y[i], y[i+1], x[i], x[i+1], x_f) + y[i])/2* ((x[i] + x_f)*0.5 - x_sc)
        #print((x[i] + x_f)*0.5 - x_sc)
        #print(x[i], x_f, total)
    return total


def integrate_1d_list_tau(x, y, x_f, x_sc):
    #inputs: x, y; lists containing the locations and values of all data points. x_f; the maximum location until which we integrate
    #outputs: x_new; a list containing all the data locatations up to x_f, and x_f if that is not already in the list. int_list; a list containing the integrated values at each location in x_new
    int_list = []
    x_new = []
    i = 1
    if x_f > x[-1]:
        x_f = x[-1]
    if x_f < x[0]:
        return [0], [x_f]
    if len(x) < 2:
        return [0], [x_f]
    int_list.append(0)
    x_new.append(x[0])
    while x[i] < x_f:
        int_list.append(integrate_1d_tau(x, y, x[i], x_sc))
        x_new.append(x[i])
        i += 1
    i -= 1
    if x_f > x[i]:
        int_list.append(integrate_1d_tau(x, y, x_f, x_sc))
        x_new.append(x_f)

    return int_list, x_new


#Define w_bar
def make_w_bar(y = AeroLoading):
    w_bar = []
    for i in range(len(x)):
        w_bar = w_bar + [integrate_1d(z, y[:,i], z[-1])]
    return w_bar

#Define tau
def make_tau(x_sc, y = AeroLoading):
    tau = []
    for i in range(len(x)):
        tau = [integrate_1d_tau(z, y[:,i], z[-1], x_sc)] + tau
    return tau

w_bar = make_w_bar()

def Integral(x_f, p = 2):
    ret_list_2 = w_bar
    x_list_2 = x
    for _ in range(p-2):
        ret_list_2, x_list_2 = integrate_1d_list(x, w_bar, x_f)
    return integrate_1d(x_list_2, ret_list_2, x_f)


def IntegralShear(x_f, z_sc, p = 2):
    tau_list_2 = make_tau(z_sc)
    x_list_2 = x
    for _ in range(p-2):
        tau_list_2, x_list_2 = integrate_1d_list(x_list_2, tau_list_2, x_f)

    return integrate_1d(x_list_2, tau_list_2, x_f)