PSAS Launch 8 was a very successfull roll control flight. Since we had good data from that flight let's use it to reconstruct coefficient of lift and see how well our model performs.
In [1]:
%matplotlib inline
import numpy as np
import utils
from matplotlib import rcParams
rcParams.update({'font.size': 14})
from math import sin, cos, radians, degrees, exp, sqrt, fabs
import numpy
columns = numpy.loadtxt('../data/Launch-8/opal/Opal_launch_rotated.csv', delimiter=',', unpack=True)
opal_time = columns[0]
opal_accel = columns[3]
opal_roll = columns[6]
opal_mag_x = columns[7]
opal_mag_y = columns[8]
opal_roll = opal_roll*57.2957795
In [2]:
fig = utils.Plot(r"Launch 8 Roll Rate", "Time [$s$]", "Rate [${}^0/s$]")
fig.plot(opal_time, opal_roll)
fig.show()
In [3]:
fig = utils.Plot(r"Launch 8 Magnetometer", "Time [$s$]", "Magnetometer Values [$T$]")
fig.plot(opal_time, opal_mag_x, label=r"$y$")
fig.plot(opal_time, opal_mag_y, color=utils.blue, label=r"$z$")
fig.legend(loc=3)
fig.show()
In [4]:
# Trick! Get the orgin to be zero by setting the mean to be zero
opal_mag_x -= opal_mag_x.mean()
opal_mag_y -= opal_mag_y.mean()
fig = utils.Plot(r"Launch 8 Magnetometer: $[x,y]$ vector", "Time [$s$]", "Magnetometer [$T$]")
fig.plot(opal_mag_x, opal_mag_y)
fig.plot(opal_mag_x[:10], opal_mag_y[:10], color='#00ff00') # mark the begining, look for a small green dot
fig.show()
In [5]:
# use numpy to build [y,z] vector from the two arrays
mag = np.dstack((opal_mag_x, opal_mag_y))[0]
# reference angle, should be a mean, but good enough.
begin = mag[0]
# storage
angles = [] # roll angle
rollup = 0 # every time we make one revolution, add 360 to this
last_angle = 0 # for checking if we've wrapped back to -pi
for v1 in mag:
# compute angle between now and the reference angle
s = np.cross(v1, begin)
c = np.dot(v1, begin)
angle = np.arctan2(s, c)
# store
angles.append(np.degrees(angle)+rollup)
last_angle = angle
# how many times did we spin?
print "The rocket made", rollup/360, "revolutions during ascent"
fig = utils.Plot(r"Launch 8 Magetometer Derived Roll Angle", "Time [$s$]", "Total Angle [${}^0$]")
fig.plot(opal_time, angles)
fig.show()
In [6]:
# differentiate
sample_rate = 128
angle_rate = np.diff(angles)
angle_rate = angle_rate*sample_rate # correct for time units using the sample rate of the sensor
# low pass the hell out of the result
from scipy.signal import firwin, lfilter
freq_cuttoff = 9.0
filter_taps = 4096
window = firwin(numtaps=filter_taps, cutoff=freq_cuttoff, nyq=sample_rate/2)
angle_rate_padded = np.concatenate((angle_rate, np.zeros(filter_taps / 2)))
angle_rate_filtered = lfilter(window, 1.0, angle_rate_padded)[filter_taps / 2:]
fig = utils.Plot(r"Launch 8 Magetometer Derived Roll Rate", "Time [$s$]", "Angle [${}^0/s$]")
fig.plot(opal_time[1:], angle_rate, color='black', alpha=0.1, label="Roll Rate")
fig.plot(opal_time[1:], angle_rate_filtered, alpha=1, lw=1.6, label="Filtered Roll Rate")
fig.plot(opal_time, opal_roll, color=utils.blue, label="Roll Gryo")
fig.ylim([-1000, 1000])
fig.legend()
fig.show()
In [7]:
columns = numpy.loadtxt('../data/Launch-8/roll/raw/output_trim.csv', delimiter=',', unpack=True, usecols=(0,1,2,3))
roll_sample_rate = 1000
roll_time = columns[0]
roll_accel = columns[1]
roll_rate = columns[2]
fin_angle = columns[3]
roll_time = roll_time - 4540.7
roll_rate = roll_rate - roll_rate[0:10].mean()
roll_rate = roll_rate / 42.9
roll_accel -= roll_accel[0:10].mean()
roll_accel *= -0.0068
fig = utils.Plot(r"Launch 8 Roll roll", "Time [$s$]", "Roll [${}^0/s$]")
fig.plot(roll_time, roll_rate, color=utils.red, label=r"roll")
fig.plot(opal_time, opal_roll, color=utils.blue, label=r"opal")
fig.legend()
fig.show()
In [8]:
freq_cuttoff = 17.0
filter_taps = 2**12
window = firwin(numtaps=filter_taps, cutoff=freq_cuttoff, nyq=roll_sample_rate/2)
roll_rate_padded = np.concatenate((roll_rate, np.zeros(filter_taps / 2)))
roll_rate_filtered = lfilter(window, 1.0, roll_rate_padded)[filter_taps / 2:]
fig = utils.Plot(r"Launch 8 Roll Rate Filtered", "Time [$s$]", "Rate [${}^0/s$]")
fig.plot(roll_time, roll_rate_filtered)
fig.show()
In [9]:
roll_rate_diff = np.diff(roll_rate_filtered)
roll_rate_diff = roll_rate_diff*roll_sample_rate
fig = utils.Plot(r"Launch 8 Roll Rate Filtered", "Time [$s$]", "Rate [${}^0/s$]")
fig.plot(roll_time[1:], roll_rate_diff)
fig.show()
In [10]:
fin_angle -= fin_angle[0:10].mean()
fin_angle /= (fin_angle.max()/15.0)
fig = utils.Plot(r"Launch 8 Fin Angle", "Time [$s$]", "Angle [${}^0$]")
fig.plot(roll_time, fin_angle)
fig.show()
In [11]:
fig = utils.Plot(r"Launch 8 Acceleration", "Time [$s$]", "Acceleration [$m/s/s$]")
fig.plot(roll_time, roll_accel, color=utils.red, label=r"roll")
fig.plot(opal_time, opal_accel, color=utils.blue, label=r"opal")
fig.legend()
fig.show()
In [12]:
from scipy.integrate import simps
velocity = [0]
altitude = [0]
for i in range(1, len(roll_accel)):
altitude.append(simps(velocity, roll_time[:i]))
velocity.append(simps(roll_accel[:i], roll_time[:i]))
In [13]:
fig = utils.Plot(r"Launch 8 Velocity", "Time [$s$]", "Velocity [$m/s$]")
fig.plot(roll_time, velocity)
fig.show()
In [14]:
fig = utils.Plot(r"Launch 8 Altitude", "Time [$s$]", "Altitude [$m$]")
fig.plot(roll_time, altitude)
fig.show()
In [15]:
# Define PSAS:
k_p = 2.45
k_v = 3.2
Cl_base = 3.2
fin_area = 1.13e-3
def C_L(a, v):
"""Find C_L for a given speed and angle of attack
:param float a: Angle of attack, alpha in degrees
:param float v: velocity, v in m/s
"""
# math is in radians
a = radians(a)
# Subsonic case
def _subsonic():
sina = sin(a)
cosa = cos(a)
cl = k_p*sina*cosa*cosa
cl += k_v*cosa*sina*sina
return cl
# Supersonic case
def _supersonic():
cl = a*Cl_base
return cl
if v <= 265:
return _subsonic()
elif v < 330:
# Intepolate between super and subsonic
y0 = _subsonic()
y1 = _supersonic()
x0 = 265
x1 = 330
cl = y0 + (y1-y0)*(v - x0)/(x1-x0)
return cl
else:
return _supersonic()
def lift(a, v, alt):
"""Compute the lift of one fin at an angle of
attack, velocity, and altitdue
:param float a: Angle of attack in degrees
:param float v: velocity in m/s
:param float alt: altitude MSL in meters
:returns float lift in Newtons:
"""
# get density of atmosphere with quick exponential model
rho = 1.2250 * exp((-9.80665 * 0.0289644 * alt)/(8.31432*288.15))
l = 0.5*C_L(a, v)*rho*v*v*fin_area
return l
def anglar_accel(a, x, v, t):
"""Compute angular accelation for a single point
:param a: Fin alpha (degrees)
:param x: altitude (meters, MSL)
:param v: air velocity (m/s)
:param t: time since launch (seconds)
"""
fin_arm = 0.085
I_0 = 0.086
I_1 = 0.077
# compute I
if t <= 0:
I = I_0
elif t < 5.6:
I = I_0 + (I_1-I_0)*t/5.6
else:
I = I_1
return degrees((4*lift(a, v, x)*fin_arm)/I)
#return 4*lift(a, v, x)*fin_arm
model_accel = []
for i, t in enumerate(roll_time):
aa = anglar_accel(fin_angle[i], 1390+altitude[i], velocity[i], t)
model_accel.append(aa)
fig = utils.Plot(r"Launch 8 Modeled Angular Accel", "Time [$s$]", "Accel [${}^0/s/s$]")
fig.plot(roll_time, model_accel)
fig.plot(roll_time[1:], roll_rate_diff, color=utils.blue)
fig.show()