GravityWaveAnalysis/pres2alt_complex.m at main · arthurMedforth/GravityWaveAnalysis · GitHub

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


function z = pres2alt_complex(pres)

z = zeros(size(pres));


% These are the levels from the wiki page for the BAROMETRIC FORMULA


% define pressure "base" levels:

p0 = 1013.25; % (hPa) (sea level)
p1 = 226.3210;
p2 = 54.7489;
p3 = 8.6802;
p4 = 1.1091;
p5 = 0.6694;
p6 = 0.0396;

% and their corresponding altitude levels (km)

z0 = 0;
z1 = 11;
z2 = 20;
z3 = 32;
z4 = 47;
z5 = 51;
z6 = 71;


% now calculate scale heights...

H0 = - (z1-z0) / (log(p1/p0));
H1 = - (z2-z1) / (log(p2/p1));
H2 = - (z3-z2) / (log(p3/p2));
H3 = - (z4-z3) / (log(p4/p3));
H4 = - (z5-z4) / (log(p5/p4));
H5 = - (z6-z5) / (log(p6/p5));

H6 = H5;


% and now calculate the altitudes from pres2alt.m...


for i = 1:length(pres),

    p = pres(i);


    %added by Corwin to make NaNsafe
    if ~isfinite(p);
      z(i) = NaN;
      continue
    end


%  %  if p >= p0, % Quick clause to cope with above normal surface pressure.
%  %              % Instead of saying things are underground, set them to zero
%  %              % and move on.
%  %
%  %      z(i) = 0;
%  %
%  %      continue
%  %
%  %  end


if p <= p0 && p > p1, % LAYER 0

    ref_pres = p0; ref_alt = z0;

    H = H0;

end

if p <= p1 && p > p2, % LAYER 1

    ref_pres = p1; ref_alt = z1;

    H = H1;

end

if p <= p2 && p > p3, % LAYER 2

    ref_pres = p2; ref_alt = z2;

    H = H2;

end

if p <= p3 && p > p4, % LAYER 3

    ref_pres = p3; ref_alt = z3;

    H = H3;

end

if p <= p4 && p > p5, % LAYER 4

    ref_pres = p4; ref_alt = z4;

    H = H4;

end

if p <= p5 && p > p6, % LAYER 5

    ref_pres = p5; ref_alt = z5;

    H = H5;

end

if p <= p6, % LAYER 6 AND ABOVE

    ref_pres = p6; ref_alt = z6;

    H = H6;

end


z(i) = ( -H * log(p/ref_pres) ) + ref_alt;


end % end cycling through all of pres's element(s)...