summaryrefslogtreecommitdiff
path: root/src/util/atmosphere.5c
blob: 9b5107f09396ea2a99a502d11c1a81bb8d2af723 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
#!/usr/bin/nickle -f
/*
 * Pressure Sensor Model, version 1.1
 *
 * written by Holly Grimes
 *
 * Uses the International Standard Atmosphere as described in
 *   "A Quick Derivation relating altitude to air pressure" (version 1.03)
 *    from the Portland State Aerospace Society, except that the atmosphere
 *    is divided into layers with each layer having a different lapse rate.
 *
 * Lapse rate data for each layer was obtained from Wikipedia on Sept. 1, 2007
 *    at site <http://en.wikipedia.org/wiki/International_Standard_Atmosphere
 *
 * Height measurements use the local tangent plane.  The postive z-direction is up.
 *
 * All measurements are given in SI units (Kelvin, Pascal, meter, meters/second^2).
 *   The lapse rate is given in Kelvin/meter, the gas constant for air is given
 *   in Joules/(kilogram-Kelvin).
 */

const real GRAVITATIONAL_ACCELERATION = -9.80665;
const real AIR_GAS_CONSTANT = 287.053;
const int NUMBER_OF_LAYERS = 7;
const real MAXIMUM_ALTITUDE = 84852;
const real MINIMUM_PRESSURE = 0.3734;
const real LAYER0_BASE_TEMPERATURE = 288.15;
const real LAYER0_BASE_PRESSURE = 101325;

/* lapse rate and base altitude for each layer in the atmosphere */
const real[NUMBER_OF_LAYERS] lapse_rate = {
	-0.0065, 0.0, 0.001, 0.0028, 0.0, -0.0028, -0.002
};
const int[NUMBER_OF_LAYERS] base_altitude = {
	0, 11000, 20000, 32000, 47000, 51000, 71000
};


/* outputs atmospheric pressure associated with the given altitude. altitudes
   are measured with respect to the mean sea level */
real altitude_to_pressure(real altitude) {

   real base_temperature = LAYER0_BASE_TEMPERATURE;
   real base_pressure = LAYER0_BASE_PRESSURE;

   real pressure;
   real base; /* base for function to determine pressure */
   real exponent; /* exponent for function to determine pressure */
   int layer_number; /* identifies layer in the atmosphere */
   int delta_z; /* difference between two altitudes */

   if (altitude > MAXIMUM_ALTITUDE) /* FIX ME: use sensor data to improve model */
      return 0;

   /* calculate the base temperature and pressure for the atmospheric layer
      associated with the inputted altitude */
   for(layer_number = 0; layer_number < NUMBER_OF_LAYERS - 1 && altitude > base_altitude[layer_number + 1]; layer_number++) {
      delta_z = base_altitude[layer_number + 1] - base_altitude[layer_number];
      if (lapse_rate[layer_number] == 0.0) {
         exponent = GRAVITATIONAL_ACCELERATION * delta_z
              / AIR_GAS_CONSTANT / base_temperature;
         base_pressure *= exp(exponent);
      }
      else {
         base = (lapse_rate[layer_number] * delta_z / base_temperature) + 1.0;
         exponent = GRAVITATIONAL_ACCELERATION /
              (AIR_GAS_CONSTANT * lapse_rate[layer_number]);
         base_pressure *= pow(base, exponent);
      }
      base_temperature += delta_z * lapse_rate[layer_number];
   }

   /* calculate the pressure at the inputted altitude */
   delta_z = altitude - base_altitude[layer_number];
   if (lapse_rate[layer_number] == 0.0) {
      exponent = GRAVITATIONAL_ACCELERATION * delta_z
           / AIR_GAS_CONSTANT / base_temperature;
      pressure = base_pressure * exp(exponent);
   }
   else {
      base = (lapse_rate[layer_number] * delta_z / base_temperature) + 1.0;
      exponent = GRAVITATIONAL_ACCELERATION /
           (AIR_GAS_CONSTANT * lapse_rate[layer_number]);
      pressure = base_pressure * pow(base, exponent);
   }

   return pressure;
}


/* outputs the altitude associated with the given pressure. the altitude
   returned is measured with respect to the mean sea level */
real pressure_to_altitude(real pressure) {

   real next_base_temperature = LAYER0_BASE_TEMPERATURE;
   real next_base_pressure = LAYER0_BASE_PRESSURE;

   real altitude;
   real base_pressure;
   real base_temperature;
   real base; /* base for function to determine base pressure of next layer */
   real exponent; /* exponent for function to determine base pressure
                             of next layer */
   real coefficient;
   int layer_number; /* identifies layer in the atmosphere */
   int delta_z; /* difference between two altitudes */

   if (pressure < 0)  /* illegal pressure */
      return -1;
   if (pressure < MINIMUM_PRESSURE) /* FIX ME: use sensor data to improve model */
      return MAXIMUM_ALTITUDE;

   /* calculate the base temperature and pressure for the atmospheric layer
      associated with the inputted pressure. */
   layer_number = -1;
   do {
      layer_number++;
      base_pressure = next_base_pressure;
      base_temperature = next_base_temperature;
      delta_z = base_altitude[layer_number + 1] - base_altitude[layer_number];
      if (lapse_rate[layer_number] == 0.0) {
         exponent = GRAVITATIONAL_ACCELERATION * delta_z
              / AIR_GAS_CONSTANT / base_temperature;
         next_base_pressure *= exp(exponent);
      }
      else {
         base = (lapse_rate[layer_number] * delta_z / base_temperature) + 1.0;
         exponent = GRAVITATIONAL_ACCELERATION /
              (AIR_GAS_CONSTANT * lapse_rate[layer_number]);
         next_base_pressure *= pow(base, exponent);
      }
      next_base_temperature += delta_z * lapse_rate[layer_number];
   }
   while(layer_number < NUMBER_OF_LAYERS - 1 && pressure < next_base_pressure);

   /* calculate the altitude associated with the inputted pressure */
   if (lapse_rate[layer_number] == 0.0) {
      coefficient = (AIR_GAS_CONSTANT / GRAVITATIONAL_ACCELERATION)
                                                    * base_temperature;
      altitude = base_altitude[layer_number]
                    + coefficient * log(pressure / base_pressure);
   }
   else {
      base = pressure / base_pressure;
      exponent = AIR_GAS_CONSTANT * lapse_rate[layer_number]
                                       / GRAVITATIONAL_ACCELERATION;
      coefficient = base_temperature / lapse_rate[layer_number];
      altitude = base_altitude[layer_number]
                      + coefficient * (pow(base, exponent) - 1);
   }

   return altitude;
}