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 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 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 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).

Pegasus III landing spot – predicted from ascent data and drag co-efficient (1Km squares)

Using drag coefficient data from Pegasus 1 - on Pegasus 3 ascent data (and allowing for payload weight differences) it's landing spot was estimated within 1.2Km of the actual. Not bad considering Peg 3 had most of its balloon still attached (affecting descent drag) and Peg 1 was a clean burst.

The algorithm works like this:

Infinite loop:

• Wait approx. 1 minute
• 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 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 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.
• Add these scaled latitude and longitude differences to a cumulative total.
• Add the cumulative totals to the current position to get the predicted landing spot - re-translating the longitude given the predicted latitude.

A very simple formula to calculate air density (kg/cu m) is

   1.205 * EXP(-height / 7990.6) 

where height is in meters. The (better) NASA algorithem for calcualting air density is given in the pseudocode below.

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 done with the formulas:

   LonAdj = Longitude * COS(Latitude)

and

   Longitude = LonAdj / COS(Latitude)

This algorithm assumes that the earth is a sphere (rather than the ovaloid that it actually is).

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_alt;
	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_alt);
	while(1)
	{
		sleep(1);
		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_alt		= this_alt - last_alt;
		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_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) * 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 * DEG2RAD)); // Translate longtitude back relative to predicted latitude 
		last_time = this_time;
		last_lat = this_lat;
		last_lon = this_lon;
		last_alt = this_alt;
	}
}