Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
integrate qu8k accelerometer data
import scipy as sp
from scipy import fft, ifft
from scipy.integrate import cumtrapz
from scipy.integrate import simps
import csv
from matplotlib import rc
fig_width_pt = 469.75
inch_per_pt = 1.0 / 72.27
figw = fig_width_pt * inch_per_pt
figh = (620.43 / 2.1) * inch_per_pt
fig_size = [figw, figh]
rc('figure', figsize=fig_size)
import matplotlib.pylab as plt
def fft_diff(y, L):
i = 1.0j
n = len(y)
k = sp.append(sp.array(range(0, n/2)), 0.0)
k = (2.0 * sp.pi / L) * sp.append(k , sp.array(range(-n/2+1,0)))
return(ifft(i * k * fft(y)).real)
def fft_int(y, x, L):
i = 1.0j
n = len(y)
k = sp.append(sp.array(range(0, n/2)), 0.0)
k = (2.0 * sp.pi / L) * sp.append(k , sp.array(range(-n/2+1,0)))
fy = fft(y)
y_int = -i * fy / k
y_int[0] = 0.0
y_int[n/2] = 0.0
return(ifft(y_int).real + (x*fy[0]/n).real)
def c_T(T):
gamma = 1.4
R = 1716.59
return(sp.sqrt(gamma*R*T))
def c_alt(h):
T0 = 518.67 # R
theta = sp.ones(h.shape)
theta[h<278385.8] = 1.237723 - h[h<278385.8] / 472687.0
theta[h<232939.6] = 1.434843 - h[h<232939.6] / 337634.0
theta[h<167322.8] = 0.939268
theta[h<154199.5] = 0.482561 + h[h<154199.5] / 337634.0
theta[h<104986.9] = 0.682457 + h[h<104986.9] / 945374.0
theta[h<65616.8] = 0.751865
theta[h<36089.2] = 1.0 - h[h<36089.2] / 145442.0
return(c_T(theta*T0))
reader = csv.reader(open("qu8k_ascent.csv","r"))
skiplines = 2
header = []
data = []
for i,row in enumerate(reader):
if(i<skiplines):
header.append(row)
else:
data.append([float(col) for col in row[0:len(row)-1]])
traj = sp.array(data)
traj = traj[0:traj.shape[0]-7,:]
traj[:,0] = traj[:,0] + 0.04
# acceleration in the frequency domain:
a_f = fft(traj[:,2])
v_fft = fft_int(traj[:,2], traj[:,0], max(traj[:,0])-min(traj[:,0]))
# set the initial condition to v=0, arbitrary constant of integration
v_fft = v_fft - v_fft[0]
d_fft = fft_int(v_fft, traj[:,0], max(traj[:,0])-min(traj[:,0]))
d_fft = d_fft - d_fft[0] + 3900.0 # black rock is at 3900 ft
mach_fft = v_fft / c_alt(d_fft)
v_trapz = sp.append(0.0, cumtrapz(traj[:,2], traj[:,0]))
fig = plt.figure()
plt.title(r'"Qu8k" at Black Rock, NV 30 Sept 2011')
ax1 = fig.add_subplot(111)
ax1.plot(traj[:,0], traj[:,2]/1e3, 'r-', label="acceleration")
#ax1.plot(traj[:,0], traj[:,3], label="velocity, rectangular")
#ax1.plot(traj[:,0], v_trapz, label="velocity, trapezoidal")
ax1.plot(traj[:,0], v_fft/1e3, 'g-', label="velocity")
ax1.plot(traj[:,0], mach_fft, 'm-', label="Mach")
xmin,xmax,ymin,ymax = ax1.axis()
ax1.set_xlim([0,xmax])
ax1.legend(loc=9)
ax1.set_xlabel(r'time $(s)$')
ax1.set_ylabel(r'$kft/s$ or $kft/s/s$ or $M$')
ax2 = ax1.twinx()
ax2.plot(traj[:,0], d_fft/1e3, 'b-', label='altitude')
ax2.set_ylabel(r'altitude $(kft)$', color='b')
ax1.text(20,0.1,'accelerometer data:\n ddeville.com/derek/Qu8k.html')
plt.savefig("qu8k_trajectory.pdf")
plt.figure()
plt.title(r'difference from FFT-based integration')
plt.semilogy(traj[:,0], abs(traj[:,3]-v_fft), label=r'1st order (rectangular)')
plt.semilogy(traj[:,0], abs(v_trapz-v_fft), label=r'2nd order (trapezoidal)')
plt.legend(loc=0)
plt.ylabel(r'velocity magnitude $(ft/s)$')
plt.xlabel(r'time $(s)$')
plt.savefig("qu8k_trajectory_error.pdf")
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.