Differences

This shows you the differences between two versions of the page.

--- ideas:landing_spot_prediction [2006/08/22 19:11] – rocketboy
+++ ideas:landing_spot_prediction [2008/07/19 23:33] (current) – external edit 127.0.0.1
@@ Line 3: / Line 3: @@
 I wanted to write an algorithm that could continuously calculate the expected payload landing spot – to end the flight early if it looked like the payload might drift outside of pre-defined bounds. The algorithm would continuously compare the predicted spot this with an internal map of the UK outline and operate a cut-down device if it looked likely that it might drift out to sea.  This would allow the balloon payload to run autonomously, reduce the possibility of it being lost and maximise on balloon flight duration.
-I reasoned that the assent data from an onboard GPS would be a good way of measuring the wind speed, direction at each level, and together data about the drag coefficient of the parachute and payload weight – the decent path could then be predicted and hence a landing spot.
+I reasoned that the ascent data from an onboard GPS would be a good way of measuring the wind speed, direction at each level, and together data about the drag coefficient of the parachute and payload weight – the decent path could then be predicted and hence a landing spot.
-I wanted an algorithm that would be easy to implement in a small micro (PIC) and would not use a lot of processor power, program or memory space.  The algorithm chosen estimates the decent drift in incremental steps – as the balloon ascends – it predicts the decent drift each minute of the assent and adds this to a cumulative total for the flight so far. The predicted landing spot is just the current position plus the calculated cumulative decent drift.
+I wanted an algorithm that would be easy to implement in a small micro (PIC) and would not use a lot of processor power, program or memory space.  The algorithm chosen estimates the decent drift in incremental steps – as the balloon ascends – it predicts the descent drift each minute of the ascent and adds this to a cumulative total for the flight so far. The predicted landing spot is just the current position plus the calculated cumulative decent drift.
-I tested the algorithm using data form the Pegasus III balloon flight. – I used the decent data (height & duration) to back calculate the balloon/payload drag coefficient/weight.  Then I used these to calculate the landing spot using just the assent data (time, latitude, longitude and height) up to apogee.  The result was much more accurate than I expected – predicting the landing spot to **within a few hundred meters**. Given the coarseness of discrete integration being used, the assumptions about co-ordinate maths, simplicity of the atmopheric model and expected changes in wind directions over the duration of the flight it is pretty amazing (IMO).
+I tested the algorithm using data form the Pegasus III balloon flight. – I used the descent data (height & duration) to back calculate the balloon/payload drag coefficient/weight.  Then I used these to calculate the landing spot using just the ascent data (time, latitude, longitude and height) up to apogee.  The result was much more accurate than I expected – predicting the landing spot to **within a few hundred meters**. Given the coarseness of discrete integration being used, the assumptions about co-ordinate maths, simplicity of the atmopheric model and expected changes in wind directions over the duration of the flight it is pretty amazing (IMO).
 {{:ideas:predict.jpg|:ideas:predict.jpg}}
@@ Line 21: / Line 21: @@
   * Compare the current time, latitude, longitude and height with the previous sample and calculate the differences of each.
   * Calculate the expected air density at the mid-point height between the current and previous sample
-  * Calculate the decent rate that would apply given the air density, drag coefficient of the parachute/payload and its weight.
+  * Calculate the descent rate that would apply given the air density, drag coefficient of the parachute/payload and its weight.
-  * Calculate how long it would take to fall from the current height to the previous height given that decent rate.
+  * Calculate how long it would take to fall from the current height to the previous height given that descent rate.
   * Scale the difference in latitude and longitude by the ratio of the time between the previous and current sample points and the calculated decent time.
   * Translate the difference in longitude using the midpoint of the current and previous sample latitude.
@@ Line 28: / Line 28: @@
   * Add the cumulative totals to the current position to get the predicted landing spot - re-translating the longitude given the predicted latitude.
-The formula to calculate air density (kg/cu m) is
+A very simple formula to calculate air density (kg/cu m) is
 .205 * EXP(-height / 7990.6)
-where height is in meters.
+where height is in meters.   The (better) NASA algorithem for calcualting air density is given in the pseudocode below.
-The formula to calculate decent rate(m/sec)  is
+The formula to calculate descent rate(m/sec)  is
      SQRT((weight * 9.81)/(0.5 * density * cd * area))
 where cd is the drag co-efficient, area is the effective area, and weight is the weight (kg) of the parachute/payload combination.
-Longitude adjustment is dome with the formulas:
+Longitude adjustment is done with the formulas:
      LonAdj = Longitude * COS(Latitude)
 and
@@ Line 43: / Line 43: @@
 ===== Co-Ordinate Maths =====
 This algorithm assumes that the earth is a sphere (rather than the ovaloid that it actually is).
 ====== Pseudo Code ======
+<code c>double weight, cd, area;
+// NASA air temperature/pressure/denisity model
+double air_density(double alt)
+{
+     if (alt < 11000.0)
+     {    // below 11Km - Troposphere
+          temp = 15.04 - (0.00649 * alt);
+          pres = 101.29 * power((temp + 273.1) / 288.08, 5.256);
+     }
+     else
+     {
+         if (alt < 25000.0)
+         { // between 11Km and 25Km - lower Stratosphere
+              temp = -56.46;
+              pres = 22.65 * exp(1.73 - ( 0.000157 * alt));
+         }
+         else
+         { // above 25Km - upper Stratosphere
+               temp = -131.21 + (0.00299 * alt);
+               pres = 2.488 * power((temp + 273.1) / 216.6, -11.388);
+         }
+     }
+     return(pres / (0.2869 * (temp + 273.1)));
+}
+void predict()
-	double this_time, this_lat, this_lon, this_height;
+{
-	double last_time, last_lat, last_lon, last_height;
+	double this_time, this_lat, this_lon, this_alt;
-	double delta_time, delat_lat, delta_lon, delta_height; // differences between readings
+	double last_time, last_lat, last_lon, last_alt;
+	double delta_time, delta_lat, delta_lon, delta_alt; // differences between readings
 	double sum_lon, sum_lat; // Cumlative (scaled and translated) totals of lat/lon deltas
 	double density, descent_rate, descent_time, scale; // transient variables
 	double predict_lat, predict_lon; // currently predicted landing spot
-	get_gps(&last_time,&last_lat,&last_lon,&last_height);
+	get_gps(&last_time,&last_lat,&last_lon,&last_alt);
 	while(1)
 	{
 		sleep(1);
-		get_gps(&this_time,&this_lat,&this_lon,&this_height);
+		get_gps(&this_time,&this_lat,&this_lon,&this_alt);
 		delta_time		= this_time - last_time;
 		delta_lon		= this_lon - last_lon;
 		delta_lat		= this_lat - last_lat;
-		delta_height	= this_height - last_height;
+		delta_alt		= this_alt - last_alt;
-		density = 1.205 * EXP(-((this_height + last_height)/2) / 7990.6); // model for air density
+		density = air_density((this_alt + last_alt)/2); // calculate average air density between altitudes
-		descent_rate = SQRT((weight * 9.81)/(0.5 * density * cd * area)); // calcualte descent rate given air density, weight, drag ...
+		descent_rate = sqrt((weight * 9.81)/(0.5 * density * cd * area)); // calcualte descent rate given air density, weight, drag ...
-		descent_time = delta_height / descent_rate; // time taken to fall through difference in heights
+		descent_time = delta_alt / descent_rate; // time taken to fall through difference in heights
 		scale = descent_time / delta_time;			// ratio of descent/ascent times
 		delta_lat *= scale;							// scale latitude and longtitude deltas by ratio of times
 		delta_lon *= scale;
-		delta_lon *= COS((this_lat + last_Lat)/2);	// Translate longtitude relative to equator
+		delta_lon *= cos(((this_lat + last_lat)/2) * DEG2RAD);	// Translate longtitude relative to equator
 		sum_lat += delta_lat;						// Cumlative total of latitude
 		sum_lon += delta_lon;						// and (translated) longtitude
 		predict_lat = this_lat + sum_lat;			// Calculate predicted landing spot
-		predict_lon = this_lon + (sum_lon / COS(predict_lat)); // Translate longtitude back relative to predicted latitude
+		predict_lon = this_lon + (sum_lon / cos(predict_lat * DEG2RAD)); // Translate longtitude back relative to predicted latitude
 		last_time = this_time;
 		last_lat = this_lat;
 		last_lon = this_lon;
-		last_height = this_height;
+		last_alt = this_alt;
 	}
+}</code>