Skip to content

Instantly share code, notes, and snippets.

@jstults
Created November 17, 2011 16:45
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save jstults/1373698 to your computer and use it in GitHub Desktop.
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