IPython Notebook

 Untitled0 Last saved: Jun 23 9:11 PM

download lunar's raw data from http://paste.ubuntu.com/5793631/ and save it as data.txt

run this command to turn hh:mm:ss into total seconds in one colum, and altitude in the other:

awk '{split($3, t, ":"); sum = 3600*t[1] + 60*t[2] + t[3]}; {print sum" "$6}' data.txt > processed.txt

then some functions we'll need:

In [27]:
    y = np.matrix(z).T 
# some helper functions
 
def kf(x, P, z, R, Q, F, H):
    
    # measurement step
    y = np.matrix(z).T - H * x
    S = H * P * H.T + R  # residual convariance
    K = P * H.T * S.I    # Kalman gain
    x = x + K*y
    I = np.matrix(np.eye(F.shape[0])) # identity matrix
    P = (I - K*H)*P
 
    # prediction step
    x = F*x
    P = F*P*F.T + Q
    
    return x, P
 
e = 2.71828183
pi = 3.1415927
SEA_LEVEL_PRESSURE = 101325.
SEA_LEVEL_DENSITY = 1.225
RADIUS_OF_EARTH = 6371009.
 
def alt2pressure(alt):
    """Return the atmospheric pressure in Pascals given an altitude in meters.
    Uses the NASA model
    """
 
    if alt < 11000:
        temperature = 15.04 - 0.00649 * alt
        return 101290 * (((temperature + 273.1) / 288.08)**5.256)
    elif alt < 25000:
        temperature = -56.46
        return 22560 * e**(1.73 - 0.000157 * alt)
    else:
        temperature = -131.21 + 0.00299 * alt
        return 2488 * (((temperature + 273.1)/216.6)**-11.388)
 
def alt2density(alt):
    """Return the atmospheric pressure in Pascals given an altitude in meters.
    Uses the NASA model
    """
 
    if alt < 11000:
        temperature = 15.04 - 0.00649 * alt
        return alt2pressure(alt)/287.0/(temperature + 273.1)
    elif alt < 25000:
        temperature = -56.46
        return alt2pressure(alt)/287.0/(temperature + 273.1)
    else:
        temperature = -131.21 + 0.00299 * alt
        return alt2pressure(alt)/287.0/(temperature + 273.1)

now we can process the data

In [28]:
plot(times, alts) # sanity check 
 
data = zip(*genfromtxt('processed.txt'))
times = data[0]
alts = data[1]
plot(times, alts) # sanity check
Out[28]:
[<matplotlib.lines.Line2D at 0x2f0da50>]
In [42]:
plot(times, vels) 
 
vels = zeros(size(alts))
 
x = matrix([[alts[0]],[(alts[0]-alts[1])/(times[1]-times[0])]])
P = matrix(''' 100. 0.; 0. 20''')
R = 20
Q = matrix(''' 100. 0. ; 0. 50''')
H = matrix(''' 1. 0.''')
 
for i in range(size(vels)-1):
    F = matrix([[1., times[i] - times[i+1]], [0., 1.]]) # sampling time isn't constant
    x, P = kf(x, P, alts[i], R, Q, F, H)
    vels[i] = x[1]
 
plot(times, vels)
Out[42]:
[<matplotlib.lines.Line2D at 0x411ee90>]

now we can calculte Cd from the velocities

In [46]:
savefig('plot.png') 
 
def drag_co(m, density, vel, area):
    return (2*m*9.81) / (density*vel*vel*area)
 
mass = 0.490 # 490g
area = 0.362 # meters square. per irc chat with lunar
 
 
cd = zeros(size(vels))
 
for i in range(size(vels)):
    cd[i] = drag_co(mass, alt2density(alts[i]), vels[i], area)
    
plot(alts, cd)
savefig('plot.png')
In [ ]: