-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathhelperfuncs.py
106 lines (79 loc) · 3.53 KB
/
helperfuncs.py
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
import math as m
sign = lambda x: m.copysign(1, x)
def machangle(mach):
return m.asin(1/mach)
def calcv(gamma, mach1, mach2):
k = m.sqrt((gamma + 1)/(gamma - 1)) # dunno what actual significance of this is
alpha1 = machangle(mach1) # first mach angle
alpha2 = machangle(mach2) # second mach angle
v1 = k * m.atan(1/(m.tan(alpha1) * k)) - (m.radians(90) - alpha1) # first expansion angle
v2 = k * m.atan(1/(m.tan(alpha2) * k)) - (m.radians(90) - alpha2) # second expansion angle
return m.degrees(v2 - v1) # additional degrees of expansion needed to go from mach1 to mach2
def binarysearch(lowerbound, upperbound, target, steps, func):
# note: only applies to monotonically increasing functions
val = (upperbound + lowerbound)/2
for i in range(steps):
if func(val) > target:
upperbound = val
val = (upperbound + lowerbound) / 2
if func(val) < target:
lowerbound = val
val = (upperbound + lowerbound) / 2
if func(val) == target:
return val
return val
# calcmach has been verified against the table, at least for gamma of 1.4
def calcmach(gamma, mach1, angle, steps=30):
if angle == 0:
return 1
f = lambda mach2: calcv(gamma, mach1, mach2)
return binarysearch(1, 100, angle, steps, f)
def calcshockemitangle(gamma, angle, mach1=-1.0, v1=-1.0):
if mach1 == -1:
mach1 = calcmach(gamma, 1, v1)
mach2 = calcmach(gamma, mach1, abs(angle))
alpha1 = m.degrees(machangle( mach1)) # first mach angle
alpha2 = m.degrees(machangle( mach2)) # second mach angle
abar = (alpha1 + alpha2)/2
return abar + .5 * abs(angle)
def calcvmax(gamma):
k = m.sqrt((gamma + 1) / (gamma - 1))
return m.degrees(m.pi/2 * (k - 1))
def calcshockpropangle(gamma, v1, theta, turningangle):
mach1 = calcmach(gamma, 1, v1)
mach2 = calcmach(gamma, 1, v1 + abs(turningangle))
alpha1 = m.degrees(machangle( mach1)) # first mach angle
alpha2 = m.degrees(machangle( mach2)) # second mach angle
abar = -sign(turningangle) * (alpha1 + alpha2) / 2
return abar + theta - turningangle/2
def alphadiff(gamma, v1, v2):
mach1 = calcmach(gamma, 1, v1)
mach2 = calcmach(gamma, 1, v2)
alpha1 = m.degrees(machangle( mach1)) # first mach angle
alpha2 = m.degrees(machangle( mach2)) # second mach angle
return alpha1 - alpha2
def calcarearatio(gamma, mach):
return 1/mach * ((2 + (gamma - 1) * mach ** 2)/(gamma + 1)) ** (.5 * ((gamma + 1)/(gamma - 1)))
def calcmachfromarearatio(gamma, ratio, steps=20):
f = lambda mach: calcarearatio(gamma, mach)
return binarysearch(1.1, 1000, ratio, steps, f)
def shockangle(gamma, v, theta, turningangle):
mach1 = calcmach(gamma, 1, v)
alpha1 = m.degrees(machangle(mach1))
return -sign(turningangle) * alpha1 + theta
def shockprop(gamma, v, theta, turningangle):
angle1 = shockangle(gamma, v, theta, turningangle)
angle2 = shockangle(gamma, v + abs(turningangle), theta + turningangle, turningangle)
return (angle1 + angle2)/2
# GAMMAS
# air: 1.4
# rocket exhaust: 1.25 is a pretty good guess
if __name__ == "__main__":
print(calcv(1.25, 1, 5))
# should differ by 8.1 degrees, because one is measuring from the x axis and one is measuring from the (post-turn) wall
# return the same value
# probably why I had that problem
# Desired Mach Resulting Mach desired expansion angle resulting exp angle
# 5 9.6 72.82 134
# 2 2.72
#Next up: a cone, or