contrails/predict.py at master · LBaritone/contrails · 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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot
import numpy as np
import re

fig = pyplot.figure()
ax = Axes3D(fig)

# Readings that do not result in a contrail prediction
rtemp = []
rpress = []
rrh = []

# Readings that result in a contrail prediction
ptemp = []
ppress = []
prh = []

# Physical constants
t0 = 273.16e0                # Absolute T
rd = 287.05e4                # Gas constant for dry air
rh2o = 461.5e4               # Gas constant for H2O
cp = 1.005e7                 # Specific heat capacity of air
bk = 1.381e-16               # Boltzman's constant
ba = 6.1115e0                # Constants for ice saturation density
bb = 23.036e0                # (Buck [J. Atmos. Sci., 20, 1527, 1981])
bc = 279.82e0
bd = 333.7
bal = 6.1121e0               # liquid saturation
bbl = 18.729e0
bcl = 257.87e0
bdl = 227.3e0

# Constants for saturation vapor pressure from Tabazadeh et al. [1997]
c1 = 18.452406985e0
c2 = -3505.1578807e0
c3 = -330918.55082e0
c4 = 12725068.262e0

delh = 42.e10                # Heat liberated per gram fuel
wexh = 1.25e0                # Water vapor mixing ratio in exhaust
effic = 0.3e0                # fraction of combustion heat converted to propulsion

total = 0
contrails = 0

outfile = open("preds/2017prediction87418.txt", "w")
with open("data/2017weather87418.txt", "r") as infile:
    for line in infile:
        data = line.split()
        if (len(data) > 10) and \
           (not data[0].startswith('--')) and \
           (not data[0].startswith('PRES')) and \
           (not data[0].startswith('hPa')) and \
           (not re.search('[^\.0-9]', data[0])) and \
           (not re.search('[^\.0-9]', data[1])) :

          total += 1                        # total valid readings
          pamb = float(data[0]) * 1.e3      # convert to bar
          rhamb = float(data[5]) / 100.


          #  To find threshold temperature, begin at high ambient temperature and
          #  steadily decrease the temperature.  At each step, calculate the maximum
          #  plume relative humidity (with respect to water) [RH].  When the peak RH
          #  exceeds 1, then the threshold temperature has been reached.
          #  Prime slplume

          tls = -31.9
          slplume = [0.1]

          itamb = 1
          idelt = 1

          while (slplume[0] < 1) :
            tls = tls - 0.1
            tt = tls + t0

            # Ambient air density
            rnair = (pamb / tt) / bk

            # Liquid sat. H2O number density
            rnsatl = 1000. * np.exp(c1 + c2/tt + c3/(tt**2) + c4/(tt**3)) / bk / tt

            # Ice sat. H2O number density
            fexp = np.exp((bb - tls/bd) * tls / (tls + bc))
            rnsati = 1.e3 * ba * fexp / bk / (tls + t0)

            # Ambient RHI
            rhiamb = rhamb * rnsatl/rnsati

            # Ambient water vapor number density
            # if the abmient RH > 1, then use RH = 1
            if (rhamb >= 0.9999) :
              rnwamb = 0.9999 * rnsatl
            else :
              rnwamb = rnsatl * rhamb

            # Find the peak plume  saturation ratio by starting with a very large
            # deltaT (= Tplume-tls) and slowly decrementing deltaT until the peak
            # saturation is found.  i.e., start with conditions very near the
            # engine exit and move downstream.
            deltaT = 100.

            # Needs to be a list to access inside while loop
            ##########################################
            ############ plume = contrail ############
            ##########################################
            slplume_prev = [0.]

            while (slplume[0] > slplume_prev[0]) :
              slplume_prev[0] = slplume[0]

              deltaT = deltaT / 1.3
              tplume = tls + deltaT
              tt = tplume + t0

              # Liquid sat. H2O number density
              rnsatl = 1000. * np.exp(c1 + c2/tt + c3/(tt**2) + c4/(tt**3)) / bk / tt

              slplume[0] = (rnwamb/rnsatl + rnair * wexh * cp * rh2o * deltaT /
                                            (rnsatl * delh * (1. - effic) * rd))

          if (rhamb >= 0.9999) :
            print "Warning! RH > 100%;  RH = " + str(rhamb * 100)

          # if tamb < Tls then contrail will form
          tamb = float(data[2])
          # outfile.write("tamb: " + str(tamb) + str((tamb - tls) < 0.) + " Tls: " + str(tls) + "\n")
          if tamb < tls :
            contrails += 1
            output = ("Pressure: " + str(pamb/1e3) + "mbar"+ " RHI: " +
                      str(rhiamb * 100) + "% Tls: " + str(tls) + "C" +
                      " Tamb: " + str(tamb) + "C ---> PREDICTED"+ "\n")
            outfile.write(output)

            # Configure prediction point for graph
            ptemp.append(tamb)
            ppress.append(pamb/1e3)
            prh.append(rhamb * 100)
          else :
            # Configure reading point for graph
            rtemp.append(tamb)
            rpress.append(pamb/1e3)
            rrh.append(rhamb * 100)

outfile.write("\nTotal Lines: " + str(total))
outfile.write("\nContrail Predictions: " + str(contrails))


ax.scatter(rtemp, rpress, rrh, c='b')
ax.scatter(ptemp, ppress, prh, c='r')
ax.set_xlabel('Ambient Temperature (C)')
ax.set_ylabel('Ambient Pressure (mbar)')
ax.set_zlabel('Ambient RH')

pyplot.show()


  #