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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
/*
* Copyright © 2019 Keith Packard <keithp@keithp.com>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*/
#include <math.h>
#include "ao-atmosphere.h"
#define GRAVITY 9.80665
/*
* 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).
*/
#define GRAVITATIONAL_ACCELERATION (-GRAVITY)
#define AIR_GAS_CONSTANT 287.053
#define NUMBER_OF_LAYERS 7
#define MAXIMUM_ALTITUDE 84852.0
#define MINIMUM_PRESSURE 0.3734
#define LAYER0_BASE_TEMPERATURE 288.15
#define LAYER0_BASE_PRESSURE 101325
/* lapse rate and base altitude for each layer in the atmosphere */
static const double lapse_rate[] = {
-0.0065, 0.0, 0.001, 0.0028, 0.0, -0.0028, -0.002
};
static const double 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
*/
double
ao_altitude_to_pressure(double altitude)
{
double base_temperature = LAYER0_BASE_TEMPERATURE;
double base_pressure = LAYER0_BASE_PRESSURE;
double pressure;
double base; /* base for function to determine pressure */
double exponent; /* exponent for function to determine pressure */
int layer_number; /* identifies layer in the atmosphere */
double 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 */
double
ao_pressure_to_altitude(double pressure)
{
double next_base_temperature = LAYER0_BASE_TEMPERATURE;
double next_base_pressure = LAYER0_BASE_PRESSURE;
double altitude;
double base_pressure;
double base_temperature;
double base; /* base for function to determine base pressure of next layer */
double exponent; /* exponent for function to determine base pressure
of next layer */
double 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;
}
|
