/*-------------------------------------------------------------------------
 *  Rocket density vs altitude calculator - version 2.0
 *-------------------------------------------------------------------------
 * This program provides a fast executing approximation of air density 
 * (kg/m**3) at various altitudes. Two tables are provided, a low altitude 
 * table (0..19800 meters) and a high altitude table (20000+ meters).
 *
 * The low altitude table simply indexes into the density chart. A 
 * temperature correction is applied below 11000 meters. After 11000m
 * ground temperature is considered unimportant. 
 * 
 * Note that the temperature correction is an approximation, accurate to 
 * within a few percent of theoretical. (plenty accurate for RASP) To 
 * disable the temperature correction (but not altitude), change the 
 * #define TEMP_CORRECTION below. 
 *
 * The high altitude table provides density data at altitudes in excess
 * of 20000 meters. The density is extrapolated from the two nearest 
 * data points. This routine is much slower than the algorithm above, but 
 * model rocket applications should not require 20000+ meter data points. 
 *
 * Usage: double air_density(altitude, temperature)
 *        double altitude - altitude in meters
 *        double temperature - temperature in degK at ground zero 
 *
 * Author: Mark C. Spiegl
 *
 * Copyright:
 * This program may be freely distributed, copied or modified provided
 * it is not used for profit. No warranty provided.
 *
 * Revisions:
 * 02/94   v1.0   Original
 * 02/94   v2.0   Add temperature function 
 *
 *------------------------------------------------------------------------
 */

#include <math.h> 

#define DELTA_ALT        100
#define OFFSET_TO_ZERO   300
#define DT_DH            0.0065   /* degK per meter */
#define SITE_ALTITUDE    0.0      /* meters (hardcoded for fast execution) */ 
#define TEMP0            288.1    /* Temp of std atmosphere at sea level */ 
#define MAX_INDEX        (sizeof(density_chart)/sizeof(double)-1)

/* 
 * 0- do not do temperature correction   
 * 1- do temperature correction 
 */
#define TEMP_CORRECTION  1 


/*********************************************************************
* This table defines air density at various altitudes. The altitude is
* measured in meters and the density in kg/m**3
*********************************************************************/

static double density_chart[] = { 
   1.26110,                     /* -300 */
   1.24927,                     /* -200 */
   1.23744,                     /* -100 */
   1.22560,                     /* 0 */
   1.21376,                     /* 100 */
   1.20209,                     /* 200 */
   1.19059,                     /* 300 */
   1.17909,                     /* 400 */
   1.16773,                     /* 500 */
   1.15657,                     /* 600 */
   1.14542,                     /* 700 */
   1.13426,                     /* 800 */
   1.12310,                     /* 900 */
   1.11223,                     /* 1000 */
   1.10141,                     /* 1100 */
   1.09058,                     /* 1200 */
   1.07976,                     /* 1300 */
   1.06903,                     /* 1400 */
   1.05855,                     /* 1500 */
   1.04807,                     /* 1600 */
   1.03767,                     /* 1700 */
   1.02752,                     /* 1800 */
   1.01713,                     /* 1900 */
   1.00678,                     /* 2000 */
   0.99697,                     /* 2100 */
   0.98693,                     /* 2200 */
   0.97684,                     /* 2300 */
   0.96703,                     /* 2400 */
   0.95722,                     /* 2500 */
   0.94745,                     /* 2600 */
   0.93798,                     /* 2700 */
   0.92851,                     /* 2800 */
   0.91904,                     /* 2900 */
   0.90957,                     /* 3000 */
   0.90010,                     /* 3100 */
   0.89063,                     /* 3200 */
   0.88150,                     /* 3300 */
   0.87253,                     /* 3400 */
   0.86374,                     /* 3500 */
   0.85462,                     /* 3600 */
   0.84563,                     /* 3700 */
   0.83684,                     /* 3800 */
   0.82865,                     /* 3900 */
   0.82016,                     /* 4000 */
   0.81103,                     /* 4100 */
   0.80219,                     /* 4200 */
   0.79362,                     /* 4300 */
   0.78550,                     /* 4400 */
   0.77711,                     /* 4500 */
   0.76875,                     /* 4600 */
   0.76063,                     /* 4700 */
   0.75251,                     /* 4800 */
   0.74440,                     /* 4900 */
   0.73629,                     /* 5000 */
   0.72841,                     /* 5100 */
   0.72063,                     /* 5200 */
   0.71285,                     /* 5300 */
   0.70507,                     /* 5400 */
   0.69729,                     /* 5500 */
   0.68952,                     /* 5600 */
   0.68195,                     /* 5700 */
   0.67451,                     /* 5800 */
   0.66706,                     /* 5900 */
   0.65982,                     /* 6000 */
   0.65270,                     /* 6100 */
   0.64526,                     /* 6200 */
   0.63800,                     /* 6300 */
   0.63090,                     /* 6400 */
   0.62379,                     /* 6500 */
   0.61685,                     /* 6600 */
   0.61009,                     /* 6700 */
   0.60332,                     /* 6800 */
   0.59656,                     /* 6900 */
   0.58979,                     /* 7000 */
   0.58303,                     /* 7100 */
   0.57627,                     /* 7200 */
   0.56951,                     /* 7300 */
   0.56303,                     /* 7400 */
   0.55660,                     /* 7500 */
   0.55017,                     /* 7600 */
   0.54375,                     /* 7700 */
   0.53742,                     /* 7800 */
   0.53133,                     /* 7900 */
   0.52524,                     /* 8000 */
   0.51915,                     /* 8100 */
   0.51307,                     /* 8200 */
   0.50698,                     /* 8300 */
   0.50096,                     /* 8400 */
   0.49521,                     /* 8500 */
   0.48924,                     /* 8600 */
   0.48319,                     /* 8700 */
   0.47744,                     /* 8800 */
   0.47190,                     /* 8900 */
   0.46646,                     /* 9000 */
   0.46071,                     /* 9100 */
   0.45515,                     /* 9200 */
   0.44975,                     /* 9300 */
   0.44433,                     /* 9400 */
   0.43892,                     /* 9500 */
   0.43350,                     /* 9600 */
   0.42810,                     /* 9700 */
   0.42285,                     /* 9800 */
   0.41777,                     /* 9900 */
   0.41270,                     /* 10000 */
   0.40763,                     /* 10100 */
   0.40255,                     /* 10200 */
   0.39748,                     /* 10300 */
   0.39253,                     /* 10400 */
   0.38780,                     /* 10500 */
   0.38278,                     /* 10600 */
   0.37760,                     /* 10700 */
   0.37218,                     /* 10800 */
   0.36677,                     /* 10900 */
   0.36131,                     /* 11000 */
   0.35574,                     /* 11100 */
   0.35016,                     /* 11200 */
   0.34466,                     /* 11300 */
   0.33941,                     /* 11400 */
   0.33417,                     /* 11500 */
   0.32896,                     /* 11600 */
   0.32389,                     /* 11700 */
   0.31881,                     /* 11800 */
   0.31378,                     /* 11900 */
   0.30905,                     /* 12000 */
   0.30431,                     /* 12100 */
   0.29959,                     /* 12200 */
   0.29503,                     /* 12300 */
   0.29046,                     /* 12400 */
   0.28590,                     /* 12500 */
   0.28150,                     /* 12600 */
   0.27711,                     /* 12700 */
   0.27271,                     /* 12800 */
   0.26848,                     /* 12900 */
   0.26426,                     /* 13000 */
   0.26003,                     /* 13100 */
   0.25612,                     /* 13200 */
   0.25223,                     /* 13300 */
   0.24834,                     /* 13400 */
   0.24444,                     /* 13500 */
   0.24056,                     /* 13600 */
   0.23667,                     /* 13700 */
   0.23307,                     /* 13800 */
   0.22952,                     /* 13900 */
   0.22597,                     /* 14000 */
   0.22241,                     /* 14100 */
   0.21886,                     /* 14200 */
   0.21531,                     /* 14300 */
   0.21201,                     /* 14400 */
   0.20879,                     /* 14500 */
   0.20559,                     /* 14600 */
   0.20249,                     /* 14700 */
   0.19945,                     /* 14800 */
   0.19640,                     /* 14900 */
   0.19336,                     /* 15000 */
   0.19031,                     /* 15100 */
   0.18728,                     /* 15200 */
   0.18433,                     /* 15300 */
   0.18145,                     /* 15400 */
   0.17858,                     /* 15500 */
   0.17580,                     /* 15600 */
   0.17309,                     /* 15700 */
   0.17039,                     /* 15800 */
   0.16777,                     /* 15900 */
   0.16523,                     /* 16000 */
   0.16270,                     /* 16100 */
   0.16016,                     /* 16200 */
   0.15762,                     /* 16300 */
   0.15509,                     /* 16400 */
   0.15262,                     /* 16500 */
   0.15025,                     /* 16600 */
   0.14789,                     /* 16700 */
   0.14558,                     /* 16800 */
   0.14338,                     /* 16900 */
   0.14118,                     /* 17000 */
   0.13898,                     /* 17100 */
   0.13679,                     /* 17200 */
   0.13459,                     /* 17300 */
   0.13243,                     /* 17400 */
   0.13041,                     /* 17500 */
   0.12838,                     /* 17600 */
   0.12639,                     /* 17700 */
   0.12452,                     /* 17800 */
   0.12266,                     /* 17900 */
   0.12081,                     /* 18000 */
   0.11894,                     /* 18100 */
   0.11708,                     /* 18200 */
   0.11523,                     /* 18300 */
   0.11336,                     /* 18400 */
   0.11150,                     /* 18500 */
   0.10966,                     /* 18600 */
   0.10797,                     /* 18700 */
   0.10627,                     /* 18800 */
   0.10458,                     /* 18900 */
   0.10306,                     /* 19000 */
   0.10154,                     /* 19100 */
   0.10003,                     /* 19200 */
   0.09850,                     /* 19300 */
   0.09698,                     /* 19400 */
   0.09546,                     /* 19500 */
   0.09394,                     /* 19600 */
   0.09241                      /* 19700 */
}; 


/*********************************************************************
* The following table provides density data for altitudes in the
* 20000+ meter range.
*********************************************************************/

static struct high_altitude_def {
   double alt;
   double density;
}                 high_alt_table[] = {
   { 19800.0,  0.09089 },
   { 20000.0,  8.8034e-2 },
   { 25000.0,  4.0016e-2 },
   { 32000.0,  1.2721e-2 },
   { 47000.0,  1.4845e-3 },
   { 53000.0,  7.1881e-4 },
   { 75000.0,  4.339e-5 },
   { 90000.0,  3.213e-6 },
   { 126000.0, 1.566e-8 },
   { 175000.0, 2.655e-10 },
   { 300000.0, 3.279e-12 },
   { 1e10,     0.0 }
};

/*********************************************************************
* Lookup function
*********************************************************************/

double
air_density(alt, temp)
double alt;
double temp;
{
   int alt_rounded;
   int i, j;
   double ret;

   /*
    * This section handles low altitude requests. Simply index into
    * density_chart[]. Multiply density stored in chart by a temperature
    * correction.
    */
   alt_rounded = (int) alt + (OFFSET_TO_ZERO + SITE_ALTITUDE);
   i = (alt_rounded / DELTA_ALT) +
      (alt_rounded % DELTA_ALT > DELTA_ALT / 2);
   if (i <= MAX_INDEX) {

#if TEMP_CORRECTION  
      ret = density_chart[i] * (TEMP0) / (temp + SITE_ALTITUDE * DT_DH); 
#else
      ret = density_chart[i]; 
#endif 

      /*
       * This section handles high altitude requests. Extrapolate data point
       * from high_alt_table[]. Assume that ground temperature is 
       * unimportant.
       */
   } else {
      for (j = 0; high_alt_table[j].alt < alt; j++);
      i = j - 1;
      ret = high_alt_table[i].density +
         (high_alt_table[i].density - high_alt_table[j].density) *
         (alt - high_alt_table[i].alt) /
         (high_alt_table[i].alt - high_alt_table[j].alt);
   }
   return (ret);
}
