I'm trying to build a sounding rocket to reach the 100km as my senior project next year. To do this I need an equation to approximate the atmospheric density at any given altitude in the rocket's ascent, so that I can better understand the engine performance requirements of the rocket. I got my atmospheric data from https://www.engineeringtoolbox.com/standard-atmosphere-d_604.html
Here is some python code with outputs:
## Atmospheric density curve fit from Engineering Toolbox US standard Atmosphere
## Design by whisperingShiba
from scipy import optimize
import numpy as np
import matplotlib.pyplot as plt
def fxn(x, a, b):
return a*np.exp(-b*x)
def polyFxn(h, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9):
return ((((((((a9*h + a8)*h + a7)*h + a6)*h + a5)*h + a4)*h + a3)*h + a2)*h + a1)*h + a0
altitudes = np.array([0,5000,10000,15000,20000,25000,30000,35000,40000,45000,50000,60000,70000,80000,90000,100000,150000,200000,250000])
density1 = np.array([23.8,20.48,17.56,14.96,12.67,10.66,8.91,7.38,5.87,4.62,3.64,2.26,1.39,.86,.56,.33,.037,.0053,.00065])
altitudeLinSpace = np.linspace(0,500000,500000)
#Popt returns a array containing constants a,b,c... etc for function 'fxn'
popt, pcov = optimize.curve_fit(fxn, altitudes, density1, p0=(1, 1e-5))
popt2, pcov2 = optimize.curve_fit(polyFxn, altitudes, density1)
print(popt)
print(popt2)
#Plots data
plt.grid(True)
plt.ylim((0,25))
plt.ylabel('Density of atmospher in slug/ft^3 * 10^-4')
plt.xlabel('Altitude, feet')
#Data from engineering toolbox
plt.plot(altitudes, density1)
#curve_fit for negative exponential
plt.plot(altitudeLinSpace, fxn(altitudeLinSpace, popt[0], popt[1]))
plt.plot(altitudeLinSpace, polyFxn(altitudeLinSpace, popt2[0],popt2[1],popt2[2],popt2[3],popt2[4],popt2[5],popt2[6],popt2[7],popt2[8],popt2[9]))
plt.show()
Here's the output of the python script, where blue is the data, orange is the negative exponential, and green is the 9th order polynomial fit:
As can be seen in the image, the 9th order polynomial fit is really good for the first part, but deviates massively past 100000 feet. The negative exponential fit has the opposite problem, having good performance as altitude approaches infinity, but having poor performance from 50000 ft to 100000 feet.
How could I best fit this data? I want to avoid super long polynomials if possible (I only have 19 data points, and the curve is already trying too hard to 'hit' all the points). I have already put this data into excel and its even worse than python's curve fit. Any advice would be greatly appreciated.
Thanks to Claude for his work. Look at this beautiful graph:
Best Answer
I should not use a strict polynomial regression.
Looking at the data $(h_i,\rho_i)$, I first had a look at the plot of $\log \left(\frac{\rho_i}{\rho_1}\right)$ as a function of $h$. This is quite nice and using a model $$\log \left(\frac{\rho}{\rho_1}\right)=\sum_{k=1}^n a_k \left(\frac {h}{10^4}\right)^k$$ I obtained a very good fit $(R^2=0.999977)$ using $n=4$. For this model, the parameters are all significant as shown below $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a_1 & -0.2457039080 & 0.006144797 & \{-0.25909227,-0.23231555\} \\ a_2 & -0.0351850842 & 0.001530225 & \{-0.03851916,-0.03185101\} \\ a_3 & +0.0021044026 & 0.000110197 & \{+0.00186430,+0.00234450\} \\ a_4 & -0.0000390562 & 0.000002385 & \{-0.00004425,-0.00003386\} \end{array}$$ Using $n>4$ leads to non-significant parameters.
So, the model would be $$\rho=23.77\, \exp \left(\sum_{k=1}^4 a_k \left(\frac {h}{10^4}\right)^k \right)$$
Below are compared the experimental values and the predicted values from this model. $$\left( \begin{array}{ccc} h & \text{measured} & \text{predicted} \\ 0 & 23.77 & 23.770 \\ 5000 & 20.48 & 20.843 \\ 10000 & 17.56 & 17.986 \\ 15000 & 14.96 & 15.296 \\ 20000 & 12.67 & 12.839 \\ 25000 & 10.66 & 10.651 \\ 30000 & 8.91 & 8.743 \\ 35000 & 7.38 & 7.112 \\ 40000 & 5.87 & 5.739 \\ 45000 & 4.62 & 4.600 \\ 50000 & 3.64 & 3.665 \\ 60000 & 2.26 & 2.297 \\ 70000 & 1.39 & 1.423 \\ 80000 & 0.86 & 0.877 \\ 90000 & 0.56 & 0.541 \\ 100000 & 0.33 & 0.335 \\ 150000 & 0.037 & 0.0366 \\ 200000 & 0.0053 & 0.00534 \\ 250000 & 0.00065 & 0.000649 \end{array} \right)$$