bbcal.py 20.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
#!/usr/bin/env python2.7

# FIXME: check to see if file exists before grab data out
# FIXME: change the order that the Thermistors from the bbcal.in file get read, also change the sign


# FIXME: MeanSigmaPred doesn't come out right
# FIXME: Mean Predicted NBB Uncertainty doesn't come put right

10
import sys
William Roberts's avatar
William Roberts committed
11
import logging
12
13
import os
from datetime import datetime
14
import numbers
15

16
import math
17
from cal_val import AERI_lib
William Roberts's avatar
William Roberts committed
18
import numpy as np
19
20

logger = logging.getLogger(__name__)
21
22


23
24
25
26
27
28
29
30
def format_time(time):
    try:
        return '{:02}:{:02.0f}'.format(math.floor(time), time * 60 % 60)
    except ValueError:
        return str(time)



William Roberts's avatar
William Roberts committed
31
32
33
def interpolate(xin, yin, xout):
    # linear interpolation for use with the emissivity
    lenxin = len(xin)
34

William Roberts's avatar
William Roberts committed
35
    i1 = np.searchsorted(xin, xout)
William Roberts's avatar
William Roberts committed
36
37
38
39
40
41
42
43
44
45
    if i1 == 0:
        i1 = 1
    if i1 == lenxin:
        i1 = lenxin - 1

    x0 = xin[i1 - 1]
    x1 = xin[i1]
    y0 = yin[i1 - 1]
    y1 = yin[i1]

William Roberts's avatar
William Roberts committed
46
47
    fit = np.polyfit((x0, x1), (y0, y1), 1)
    val = np.polyval(fit, xout)
William Roberts's avatar
William Roberts committed
48
    return val
49
50


51
def process(directory, records, abb_breakout=None, output_path=''):
52
53
54
55
56
57
58
59
    try:
        NBBtemp = float(str.strip(os.environ.get('NBBTEMP')))
    except TypeError:
        NBBtemp = None
    try:
        IBBtemp = float(str.strip(os.environ.get('IBBTEMP')))
    except TypeError:
        IBBtemp = None
60
    return process_helper(1, directory, records, NBBtemp, IBBtemp, abb_breakout, output_path), process_helper(2, directory, records, NBBtemp, IBBtemp, abb_breakout, output_path)
Marco Kurzynski's avatar
Marco Kurzynski committed
61

62
def process_helper(channel_number, directory, records,  NBBtemp=None, IBBtemp=None, abb_breakout=None, output_path=''):
63
64
65
    if isinstance(records, numbers.Number):
        records = [records]

William Roberts's avatar
William Roberts committed
66
67
    # Maximun number of points
    MAXPTS = 8192
68

William Roberts's avatar
William Roberts committed
69
70
    # Maximum number of ranges
    MAXRNGS = 10
71

72
    if channel_number != 1 and channel_number != 2:
Marco Kurzynski's avatar
Marco Kurzynski committed
73
        logger.error("Invalid channel_number, must be 1 or 2")
74
        exit(1)
75

76
77
78
    try:
        LOGFILE = str.strip(os.environ.get(f'C{channel_number}_LOGFILE'))
    except TypeError:
79
        LOGFILE = f'bbcal3p_{os.path.basename(directory)[2:]}.lo{channel_number}'
80
    LOGFILE = os.path.join(output_path, LOGFILE)
81

William Roberts's avatar
William Roberts committed
82
    # Open logfile
83
    lf = open(LOGFILE, 'w')
William Roberts's avatar
William Roberts committed
84
85
    lf.write(" ----------------  bbcal3p.py -----------------------\n")
    lf.write(" -----------------------------------------------------\n")
86

William Roberts's avatar
William Roberts committed
87
    # Input File
88
    INFILE = os.path.join(directory, f"{os.path.basename(directory)[2:8]}C{channel_number}.RNC")
William Roberts's avatar
William Roberts committed
89
90
    lf.write("INPUT FILE: %s\n" % INFILE)

91
    lf.write("INPUT RECORD RANGE: {0}\n".format(records))
William Roberts's avatar
William Roberts committed
92
93

    # NBB & IBB Scene IDs
94
    IDNBB = str.strip(os.environ.get('IDNBB'))
William Roberts's avatar
William Roberts committed
95
    lf.write("NBB identifier: %s\n" % IDNBB)
96
    IDIBB = str.strip(os.environ.get('IDIBB'))
William Roberts's avatar
William Roberts committed
97
98
99
    lf.write("IBB identifier: %s\n" % IDIBB)

    # NBB & IBB Emissivity Files
100
    NBEMFL = str.strip(os.environ.get('NBEMFL'))
William Roberts's avatar
William Roberts committed
101
    lf.write("NBB EMISSIVITY FILE: %s\n" % NBEMFL)
102
    IBEMFL = str.strip(os.environ.get('IBEMFL'))
William Roberts's avatar
William Roberts committed
103
104
105
    lf.write("IBB EMISSIVITY FILE: %s\n" % IBEMFL)

    # Read NBB & IBB Emissivity Files
106
107
108
    d = os.path.dirname(sys.modules[__name__].__file__)
    emfp  = open(os.path.join(d, "emissivity/"+NBEMFL), 'r')
    imfp  = open(os.path.join(d, "emissivity/"+IBEMFL), 'r')
William Roberts's avatar
William Roberts committed
109
110
111
112
113
114
115
116
117
118
119
120

    # Read first line for Emissivity Files
    lf.write("FIRST LINE OF NBB EMISSIVITY-SPECTRUM FILE:\n")
    lf.write(emfp.readline())
    lf.write("FIRST LINE OF IBB EMISSIVITY-SPECTRUM FILE:\n")
    lf.write(imfp.readline())

    # Read the rest of the NBB & IBB Emissivity Files
    pairs = [x.split() for x in emfp.readlines()]
    ipairs = [x.split() for x in imfp.readlines()]

    # Emissivity Wavenamer and Values Arrays
William Roberts's avatar
William Roberts committed
121
122
    EMISSX = np.array(pairs)[:, 0].squeeze()
    EMISSY = np.array(pairs)[:, 1].squeeze()
William Roberts's avatar
William Roberts committed
123
124

    # Emissivity Wavenamer and Values Arrays
William Roberts's avatar
William Roberts committed
125
126
    IEMISSX = np.array(ipairs)[:, 0].squeeze()
    IEMISSY = np.array(ipairs)[:, 1].squeeze()
William Roberts's avatar
William Roberts committed
127
128
129
130
131
132

    # Close Emissivity Files
    emfp.close()
    imfp.close()

    # BB Cavity Factor
133
    CAVFAC2 = float(str.strip(os.environ.get('CAVFAC2')))
William Roberts's avatar
William Roberts committed
134
135
136
137
    lf.write("BB CAVITY FACTOR: %f\n" % CAVFAC2)
    RCF = 1.0 / CAVFAC2

    # NBB Apex, Bottom and Top weightings
138
139
140
    WTNBA = float(str.split(str.strip(os.environ.get('WTNB')))[0])
    WTNBB = float(str.split(str.strip(os.environ.get('WTNB')))[1])
    WTNBT = float(str.split(str.strip(os.environ.get('WTNB')))[2])
William Roberts's avatar
William Roberts committed
141
142
143
144
145

    # Write Weighting to Logfile
    lf.write("NBB Weights: %1.5f %1.5f %1.5f\n" % (WTNBA, WTNBB, WTNBT))

    # NBB Apex, Bottom and Top offsets
146
147
148
    TOSNBA = float(str.split(str.strip(os.environ.get('TOSNB')))[0])
    TOSNBB = float(str.split(str.strip(os.environ.get('TOSNB')))[1])
    TOSNBT = float(str.split(str.strip(os.environ.get('TOSNB')))[2])
William Roberts's avatar
William Roberts committed
149
150
151
152
153

    # Write Offsets to Logfile
    lf.write("NBB Offsets: %1.5f %1.5f %1.5f\n" % (TOSNBA, TOSNBB, TOSNBT))

    # IBB Apex, Bottom and Top weightings
154
155
156
    WTIBA = float(str.split(str.strip(os.environ.get('WTIB')))[0])
    WTIBB = float(str.split(str.strip(os.environ.get('WTIB')))[1])
    WTIBT = float(str.split(str.strip(os.environ.get('WTIB')))[2])
William Roberts's avatar
William Roberts committed
157
158
159
160
161

    # Write Weighting to Logfile
    lf.write("IBB Weights: %1.5f %1.5f %1.5f\n" % (WTIBA, WTIBB, WTIBT))

    # IBB Apex, Bottom and Top offsets
162
163
164
    TOSIBA = float(str.split(str.strip(os.environ.get('TOSIB')))[0])
    TOSIBB = float(str.split(str.strip(os.environ.get('TOSIB')))[1])
    TOSIBT = float(str.split(str.strip(os.environ.get('TOSIB')))[2])
William Roberts's avatar
William Roberts committed
165
166
167
168
169

    # Write Offsets to Logfile
    lf.write("IBB Offsets: %1.5f %1.5f %1.5f\n" % (TOSIBA, TOSIBB, TOSIBT))

    # MLT file
170
171
172
    try:
        MLTFILE = str.strip(os.environ.get(f'C{channel_number}_MLTFILE'))
    except TypeError:
173
        MLTFILE = f'bbcal_{os.path.basename(directory)[2:]}.mlt{channel_number}'
174
    MLTFILE = os.path.join(output_path, MLTFILE)
William Roberts's avatar
William Roberts committed
175
176
177
178
179

    # Write MLT file to Log file
    lf.write("ASCII Output (MLT) file: % s\n" % MLTFILE)

    # Summary File
180
181
182
    try:
        SUMFILE = str.strip(os.environ.get(f'C{channel_number}_SUMFILE'))
    except TypeError:
183
        SUMFILE = f'bbcal_{os.path.basename(directory)[2:]}.sum{channel_number}'
184
    SUMFILE = os.path.join(output_path, SUMFILE)
William Roberts's avatar
William Roberts committed
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
    lf.write("TEXT OUTPUT (SUMMARY) file: %s\n" % SUMFILE)

    # Array of Metadata required
    IHWDS = []
    # HWDS 0 = secene Mirror Position (Dummy postion)
    IHWDS.append("sceneMirrorPosition")
    # HWDS 1 = secene Mirror Position
    IHWDS.append("sceneMirrorPosition")
    # HWDS 2 = missingDataFlag
    IHWDS.append("missingDataFlag")
    # HWDS 3 = IBBapexTemp
    IHWDS.append("IBBapexTemp")
    # HWDS 4 = IBBbottomTemp
    IHWDS.append("IBBbottomTemp")
    # HWDS 5 = IBBtopTemp
    IHWDS.append("IBBtopTemp")

    # Using breakout cable we will be looking in
    # ABB apex and ABB bottom for 3rd Body temperatures

    if abb_breakout != None:
        # HWDS 6 = NBBapexTemp
        IHWDS.append("ABBapexTemp")
        # HWDS 7 = NBBbottomTemp
        IHWDS.append("ABBbottomTemp")

    # normally look for values in NBB apex and bottom
    else:
        # HWDS 6 = NBBapexTemp
        IHWDS.append("NBBapexTemp")
        # HWDS 7 = NBBbottomTemp
        IHWDS.append("NBBbottomTemp")

    # HWDS 8 = NBBtopTemp
    IHWDS.append("NBBtopTemp")
    # HWDS 9 = BBsupportStructureTemp
    IHWDS.append("BBsupportStructureTemp")
    # HWDS 10 = timeHHMMSS
223
224
225
    IHWDS.append("Time")
    # HWDS 11 = instrumentUnitNumber
    IHWDS.append("instrumentUnitNumber")
William Roberts's avatar
William Roberts committed
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251

    # Read Required Metadata
    HWDS = AERI_lib.getMetaValues(INFILE, IHWDS)

    # back out with ABB coefficients and process with NBB coefficients
    if abb_breakout != None:
        old_abb_aa = abb_breakout["old_abb_aa"]
        old_abb_ab = abb_breakout["old_abb_ab"]
        old_abb_ac = abb_breakout["old_abb_ac"]

        old_abb_ba = abb_breakout["old_abb_ba"]
        old_abb_bb = abb_breakout["old_abb_bb"]
        old_abb_bc = abb_breakout["old_abb_bc"]

        new_nbb_aa = abb_breakout["new_nbb_aa"]
        new_nbb_ab = abb_breakout["new_nbb_ab"]
        new_nbb_ac = abb_breakout["new_nbb_ac"]

        new_nbb_ba = abb_breakout["new_nbb_ba"]
        new_nbb_bb = abb_breakout["new_nbb_bb"]
        new_nbb_bc = abb_breakout["new_nbb_bc"]

        # Covert temperature to resitance using old coefficients
        ar_all = []
        for old_apex_temp in HWDS[6]:
            alpha_a = (old_abb_aa - (1 / old_apex_temp)) / old_abb_ac
252
253
254
            beta_a = np.sqrt((np.pow((old_abb_ab / (3 * old_abb_ac)), 3) + (np.pow(alpha_a, 2) / 4)))
            ar_all.append(np.pow(np.e, (
                    np.pow((beta_a - (alpha_a / 2)), (1.0 / 3)) - np.pow((beta_a + (alpha_a / 2)), (1.0 / 3)))))
William Roberts's avatar
William Roberts committed
255
256
257
258

        br_all = []
        for old_bot_temp in HWDS[7]:
            alpha_b = (old_abb_ba - (1 / old_bot_temp)) / old_abb_bc
259
260
261
            beta_b = np.sqrt((np.pow((old_abb_bb / (3 * old_abb_bc)), 3) + (np.pow(alpha_b, 2) / 4)))
            br_all.append(np.pow(np.e, (
                    np.pow((beta_b - (alpha_b / 2)), (1.0 / 3)) - np.pow((beta_b + (alpha_b / 2)), (1.0 / 3)))))
William Roberts's avatar
William Roberts committed
262
263
264
265
266

        # Calculate new temperature
        nbb_apex_temps = []
        for ar in ar_all:
            nbb_apex_temps.append(
267
                1 / (new_nbb_aa + new_nbb_ab * (np.log(ar)) + new_nbb_ac * np.pow(np.log(ar), 3)))
William Roberts's avatar
William Roberts committed
268
269
270
        nbb_bot_temps = []
        for br in br_all:
            nbb_bot_temps.append(
271
                1 / (new_nbb_ba + new_nbb_bb * (np.log(br)) + new_nbb_bc * np.pow(np.log(br), 3)))
William Roberts's avatar
William Roberts committed
272

William Roberts's avatar
William Roberts committed
273
274
275
        HWDS[6] = np.array(nbb_apex_temps)
        HWDS[7] = np.array(nbb_bot_temps)
        HWDS[8] = np.array(nbb_bot_temps)
William Roberts's avatar
William Roberts committed
276

William Roberts's avatar
William Roberts committed
277
278
279
    HWDS[6][np.isnan(HWDS[6])] = 0
    HWDS[7][np.isnan(HWDS[7])] = 0
    HWDS[8][np.isnan(HWDS[8])] = 0
William Roberts's avatar
William Roberts committed
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297

    # Record the first and last times
    TIME1 = HWDS[10][0]
    TIME2 = HWDS[10][-1]

    # Apply offset, weight and average NBB thermistors
    AVGNBT = (((HWDS[6] + TOSNBA) * WTNBA) + ((HWDS[7] + TOSNBB) * WTNBB) + ((HWDS[8] + TOSNBT) * WTNBT))

    # Apply offset, weight and average IBB thermistors
    AVGIBT = (((HWDS[3] + TOSIBA) * WTIBA) + ((HWDS[4] + TOSIBB) * WTIBB) + ((HWDS[5] + TOSIBT) * WTIBT))

    # Overwrite the NBB & IBB Temps with values from Blue Rack
    if NBBtemp != None:
        AVGNBT = [float(NBBtemp) for I in AVGNBT]
    if IBBtemp != None:
        AVGIBT = [float(IBBtemp) for I in AVGIBT]

    # From first record to last record
298
    for NR in records:
William Roberts's avatar
William Roberts committed
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
        # FIXME: If you can't read from the file, exit loop
        # Check BAD DATA FLAG
        if HWDS[2][NR - 1] != 0:
            lf.write("BAD DATA FLAG FOUND, REC.NO.: %d\n" % NR)
        else:
            lf.write("FROM INPUT RECORD NUMBER: %d\n" % NR)
            # write uncorrected and weighted NBB temps
            lf.write("\tNBB raw temps (apex, bottom, top): %3.3f, %3.3f, %3.3f\n" % (
                HWDS[6][NR - 1], HWDS[7][NR - 1], HWDS[8][NR - 1]))
            # Write NBB Average temperature (weighted and offset already applied)
            lf.write("\tNBB AVERAGE temperature: %3.3f\n" % (AVGIBT[NR - 1]))

            # write uncorrected and weighted IBB temps
            lf.write("\tIBB raw temps (apex, bottom, top): %3.3f, %3.3f, %3.3f\n" % (
                HWDS[3][NR - 1], HWDS[4][NR - 1], HWDS[5][NR - 1]))
            # Write IBB Average temperature (weighted and offset already applied)
            lf.write("\tIBB AVERAGE temperature: %3.3f\n" % (AVGIBT[NR - 1]))
    lf.close()

    # convert from paint emissivity to cavity emissivity
    EMISSY = [float(emi) for emi in EMISSY]
    EMISSY = [(pe / (pe + (1 - pe) / CAVFAC2)) for pe in EMISSY]
    IEMISSY = [float(emi) for emi in IEMISSY]
    IEMISSY = [(pe / (pe + (1 - pe) / CAVFAC2)) for pe in IEMISSY]

    # interpolate "NBB" & "IBB" EMISSIVITIES TO AERI Wavenumbers
    interx = AERI_lib.get_wavenumbers(INFILE, 1)

    EMISSX = [float(emi) for emi in EMISSX]
    IEMISSX = [float(emi) for emi in IEMISSX]

330
331
    BBEN = np.array([interpolate(EMISSX, EMISSY, y) for y in interx])
    iBBEN = np.array([interpolate(IEMISSX, IEMISSY, y) for y in interx])
William Roberts's avatar
William Roberts committed
332
333

    # set up arrays
William Roberts's avatar
William Roberts committed
334
335
336
337
338
339
340
341
342
    AVRADO = np.zeros(len(interx))
    SDRADO = np.zeros(len(interx))
    AVRADP = np.zeros(len(interx))
    SDRADP = np.zeros(len(interx))
    AVGOBT = np.zeros(len(interx))
    AVGPBT = np.zeros(len(interx))
    DIFFBT = np.zeros(len(interx))
    UNCBTO = np.zeros(len(interx))
    UNCBTP = np.zeros(len(interx))
William Roberts's avatar
William Roberts committed
343
    # set up arrays
William Roberts's avatar
William Roberts committed
344
345
346
347
348
349
350
351
352
    iAVRADO = np.zeros(len(interx))
    iSDRADO = np.zeros(len(interx))
    iAVRADP = np.zeros(len(interx))
    iSDRADP = np.zeros(len(interx))
    iAVGOBT = np.zeros(len(interx))
    iAVGPBT = np.zeros(len(interx))
    iDIFFBT = np.zeros(len(interx))
    iUNCBTO = np.zeros(len(interx))
    iUNCBTP = np.zeros(len(interx))
William Roberts's avatar
William Roberts committed
353
354
355
356
357
358
359
360

    # INITIALIZE SAMPLE COUNT ...
    NSMPLS = 0
    NPVALS = 0
    iNSMPLS = 0
    iNPVALS = 0

    # From first record to last record
361
    for NR in records:
William Roberts's avatar
William Roberts committed
362
363
364
365
366
        SMP = chr(int(HWDS[1][NR - 1]))
        if SMP != IDNBB and SMP != IDIBB:
            logger.info("SCENE-MIR.POS.: %s" % SMP)
            logger.info("RECORD NUMBER: %d" % NR)
        elif SMP == IDNBB:
367
368
369
370
371
372
373
374
375
376
377
378
379
            RAD = np.array(AERI_lib.get_rec_rad(INFILE, NR))

            # Accumulate SUMS (AND & SUMS OF SQUARES) OF OVSERVED BB RADIANCES ...
            AVRADO = AVRADO + RAD
            SDRADO = SDRADO + RAD ** 2

            # ACCUMULATE SUMS (& SUMS OF SQUARES) OF PREDICTED BB RADIANCES ...
            # REFRAD: REFlected RaDiance ....
            REFRAD = np.vectorize(AERI_lib.bt2rad)(interx, HWDS[9][NR - 1])
            RADNBB = BBEN * np.vectorize(AERI_lib.bt2rad)(interx, AVGNBT[NR - 1]) + (1.0 - BBEN) * REFRAD
            AVRADP = AVRADP + RADNBB
            SDRADP = SDRADP + RADNBB ** 2

William Roberts's avatar
William Roberts committed
380
381
382
            NSMPLS = NSMPLS + 1
            NPVALS = NPVALS + 1
        elif SMP == IDIBB:
383
384
385
386
387
388
389
390
391
392
393
394
            RAD = np.array(AERI_lib.get_rec_rad(INFILE, NR))

            # Accumulate SUMS (AND & SUMS OF SQUARES) OF OVSERVED BB RADIANCES ...
            iAVRADO = iAVRADO + RAD
            iSDRADO = iSDRADO + RAD ** 2

            # ACCUMULATE SUMS (& SUMS OF SQUARES) OF PREDICTED BB RADIANCES ...
            # REFRAD: REFlected RaDiance ....
            REFRAD = np.vectorize(AERI_lib.bt2rad)(interx, HWDS[9][NR - 1])
            RADNBB = iBBEN * np.vectorize(AERI_lib.bt2rad)(interx, AVGIBT[NR - 1]) + (1.0 - iBBEN) * REFRAD
            iAVRADP = iAVRADP + RADNBB
            iSDRADP = iSDRADP + RADNBB ** 2
William Roberts's avatar
William Roberts committed
395
396
397

            iNSMPLS = iNSMPLS + 1
            iNPVALS = iNPVALS + 1
William Roberts's avatar
William Roberts committed
398

399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
    # COMPUTE & OUTPUT SPECTRUM OF MEAN OBSERVED B.T. ...
    AVRADO = AVRADO / NSMPLS
    AVGOBT = np.vectorize(AERI_lib.rad2bt)(interx, AVRADO)
    iAVRADO = iAVRADO / iNSMPLS
    iAVGOBT = np.vectorize(AERI_lib.rad2bt)(interx, iAVRADO)

    # COMPUTE & OUTPUT SPECTRUM OF MEAN PREDICTED B.T. ...
    AVRADP = AVRADP / NPVALS
    AVGPBT = np.vectorize(AERI_lib.rad2bt)(interx, AVRADP)
    iAVRADP = iAVRADP / iNPVALS
    iAVGPBT = np.vectorize(AERI_lib.rad2bt)(interx, iAVRADP)

    # COMPUTE & OUTPUT SPECTRUM OF DIFFERNCE (OBS.-PRED.) B.T. ...
    DIFFBT = AVGOBT - AVGPBT
    iDIFFBT = iAVGOBT - iAVGPBT

    # COMPUTE & OUTPUT SPECTRUM OF UNCERTAINTY IN OBSERVED B.T. ...
    SDRADO = SDRADO / NSMPLS
    SDRADO = abs(SDRADO - AVRADO ** 2)
    SDRADO = np.sqrt(SDRADO)
    ARG = NSMPLS
    AVPLUS = AVRADO + SDRADO / np.sqrt(ARG)
    UNCBTO = np.vectorize(AERI_lib.rad2bt)(interx, AVPLUS) - np.vectorize(AERI_lib.rad2bt)(interx, AVRADO)
    iSDRADO = iSDRADO / iNSMPLS
    iSDRADO = abs(iSDRADO - iAVRADO ** 2)
    iSDRADO = np.sqrt(iSDRADO)
    ARG = iNSMPLS
    AVPLUS = iAVRADO + iSDRADO / np.sqrt(ARG)
    iUNCBTO = np.vectorize(AERI_lib.rad2bt)(interx, AVPLUS) - np.vectorize(AERI_lib.rad2bt)(interx, iAVRADO)

    # COMPUTE & OUTPUT SPECTRUM OF UNCERTAINTY IN PREDICTED B.T. ...
    SDRADP = SDRADP / NPVALS
    SDRADP = abs(SDRADP - AVRADP ** 2)
    SDRADP = np.sqrt(SDRADP)
    ARG = NPVALS
    AVPLUS = AVRADP + SDRADP / np.sqrt(ARG)
    UNCBTP = np.vectorize(AERI_lib.rad2bt)(interx, AVPLUS) - np.vectorize(AERI_lib.rad2bt)(interx, AVRADP)
    iSDRADP = iSDRADP / iNPVALS
    iSDRADP = abs(iSDRADP - iAVRADP ** 2)
    iSDRADP = np.sqrt(iSDRADP)
    ARG = iNPVALS
    AVPLUS = iAVRADP + iSDRADP / np.sqrt(ARG)
    iUNCBTP = np.vectorize(AERI_lib.rad2bt)(interx, AVPLUS) - np.vectorize(AERI_lib.rad2bt)(interx, iAVRADP)
William Roberts's avatar
William Roberts committed
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463

    logger.debug("Writing MLT file {0}".format(MLTFILE))
    mlt = open(MLTFILE, 'w')
    mlt.write('#MLT format\n')
    mlt.write('#Software: BBCAL3P.EXE\n')
    mlt.write('#label 1 Wavenumber\n')
    mlt.write('#label 2 Mean Observed NBB Brightness Temperature\n')
    mlt.write('#label 3 Mean Predicted NBB Brightness Temperature\n')
    mlt.write('#label 4 NBB Observed B.T. - Predicted B.T.\n')
    mlt.write('#label 5 Mean Observed NBB Uncertainty\n')
    mlt.write('#label 6 Mean Predicted NBB Uncertainty\n')
    mlt.write('#label 7 Mean Observed NBB Brightness Temperature\n')
    mlt.write('#label 8 Mean Predicted NBB Brightness Temperature\n')
    mlt.write('#label 9 NBB Observed B.T. - Predicted B.T.\n')
    mlt.write('#label 10 Mean Observed NBB Uncertainty\n')
    mlt.write('#label 11 Mean Predicted NBB Uncertainty\n')
    for i in range(len(interx)):
        mlt.write("%10.3f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f\n" % (
            interx[i], AVGOBT[i], AVGPBT[i], DIFFBT[i], UNCBTO[i], UNCBTP[i], iAVGOBT[i], iAVGPBT[i], iDIFFBT[i],
            iUNCBTO[i], iUNCBTP[i]))
    mlt.close()

Marco Kurzynski's avatar
Marco Kurzynski committed
464
465
466
467
468
469
470
    # Read number of ranges
    NRNGS = int(str.strip(os.environ.get('NRNGS')))

    # FIXME: Do something in the future if there is more then 1 range
    if NRNGS >= 2:
        logger.error("Add more code :)")

471
    WN = str.strip(os.environ.get(f'C{channel_number}_WN'))
William Roberts's avatar
William Roberts committed
472
    # Lower bound of spectral range
473
    WNLO = int(str.split(WN)[0])
William Roberts's avatar
William Roberts committed
474
    # Higher bound of spectral range
475
    WNHI = int(str.split(WN)[1])
William Roberts's avatar
William Roberts committed
476
477
478
479
480
481
482
483
484
485

    XMIN = interx[0]
    XMAX = interx[len(interx) - 1]

    DX = (XMAX - XMIN) / len(interx)

    IWNLO = int((WNLO - XMIN) / DX)
    IWNHI = int((WNHI - XMIN) / DX)

    # CALCULATE MEAN AND STD. DEV. OF DIFFERENCE (OBS-PRED) B.T. ...
486
    SMPLS = (IWNHI + 1 - IWNLO)
William Roberts's avatar
William Roberts committed
487
488
    AVUNCBTP = 0
    iAVUNCBTP = 0
489
490
491
492
493
494
    AVOMPBT = DIFFBT[IWNLO:IWNHI + 1].sum() / float(SMPLS)
    SDOMPBT = (DIFFBT ** 2)[IWNLO:IWNHI + 1].sum() / float(SMPLS)
    AVUNCBTO = UNCBTO[IWNLO:IWNHI + 1].sum() / float(SMPLS)
    iAVOMPBT = iDIFFBT[IWNLO:IWNHI + 1].sum() / float(SMPLS)
    iSDOMPBT = (iDIFFBT ** 2)[IWNLO:IWNHI + 1].sum() / float(SMPLS)
    iAVUNCBTO = iUNCBTO[IWNLO:IWNHI + 1].sum() / float(SMPLS)
William Roberts's avatar
William Roberts committed
495
496
497
498
499
500

    # don't divide by number of samples
    SDOMPBT = SDOMPBT - AVOMPBT ** 2
    iSDOMPBT = iSDOMPBT - iAVOMPBT ** 2

    # Uncertainty in the mean:
501
502
    SDOMPBT = np.sqrt(SDOMPBT) / np.sqrt(SMPLS)
    iSDOMPBT = np.sqrt(iSDOMPBT) / np.sqrt(SMPLS)
William Roberts's avatar
William Roberts committed
503
504
505

    logger.debug("Writing summary file {0}".format(SUMFILE))

506
507
    ihms = format_time(TIME1)
    jhms = format_time(TIME2)
William Roberts's avatar
William Roberts committed
508
509
    # jhms=ihms

510
    dateYYMMDD = os.path.basename(directory)[2:8]
William Roberts's avatar
William Roberts committed
511
512
513
    with open(SUMFILE, 'w') as sf:
        sf.write("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n\n")
        sf.write("Mean and Standard deviation of Observed minus Predicted BT\n\n\n")
514
515
        sf.write("Instrument: AERI-%s\n" % (int(HWDS[11][records[0] - 1])))
        sf.write('Date of measurements: 20{0}-{1}-{2}\n'.format(dateYYMMDD[:2], dateYYMMDD[2:4], dateYYMMDD[4:]))
William Roberts's avatar
William Roberts committed
516
        sf.write("Time: %s - %s\n" % (ihms, jhms))
517
        sf.write("Input Records: {0}\n\n\n".format(records))
William Roberts's avatar
William Roberts committed
518
519
520
521
        sf.write("   Vmean   DIFFmean    RSSmean  MeanSigmaObs  MeanSigmaPred  SigmaWave     V1      V2\n")

        VMEAN = 0.5 * (WNLO + WNHI)
        arg = AVUNCBTO ** 2 + (10 * AVUNCBTP) ** 2 + SDOMPBT ** 2
522
        RSSmean = np.sqrt(arg)
William Roberts's avatar
William Roberts committed
523
        arg = iAVUNCBTO ** 2 + (10 * iAVUNCBTP) ** 2 + iSDOMPBT ** 2
524
        iRSSmean = np.sqrt(arg)
William Roberts's avatar
William Roberts committed
525
526
        sf.write("NBB\n\n")
        sf.write("%8.2f %10.5f %10.5f %10.5f %13.5f %13.5f %9.2f %7.2f\n\n" % (
527
            VMEAN, AVOMPBT, RSSmean, AVUNCBTO, AVUNCBTP, SDOMPBT, WNLO, WNHI))
William Roberts's avatar
William Roberts committed
528
529
530
        sf.write("   Vmean   DIFFmean    RSSmean  MeanSigmaObs  MeanSigmaPred  SigmaWave     V1      V2\n")
        sf.write("IBB\n\n")
        sf.write("%8.2f %10.5f %10.5f %10.5f %13.5f %13.5f %9.2f %7.2f\n\n" % (
531
            VMEAN, iAVOMPBT, iRSSmean, iAVUNCBTO, iAVUNCBTP, iSDOMPBT, WNLO, WNHI))
William Roberts's avatar
William Roberts committed
532
533
534

        sf.write("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n")

535
536
537
538
    if not np.isnan(TIME1):
        TIME1 = int(TIME1)
    if not np.isnan(TIME2):
        TIME2 = int(TIME2)
William Roberts's avatar
William Roberts committed
539
    # return  wavenumbers, Mean Observed BrightNess Temp, Predicted Temp, Diff Mean, Time Range, Record Range
540
    return [interx, iAVGOBT, iAVGPBT, AVGOBT, AVGPBT, AVOMPBT, iAVOMPBT, TIME1, TIME2, records]
William Roberts's avatar
William Roberts committed
541
542
543
544
545
546
547
548
549


if __name__ == '__main__':
    doc = """
python bbcal3p.py inputFile

Outputs LOG, MLT and Summary files

"""
550

William Roberts's avatar
William Roberts committed
551
    from sys import argv
552

William Roberts's avatar
William Roberts committed
553
554
555
556
    if len(argv) == 2:
        process(argv[1])
    else:
        logger.error(doc)