' - line = line + " No Valid data for last 2 hours." - line = line + '
\n' - return line -# -#--- there are good data -# - p56 = p2dat[ind] - p56 = p56[p56 >0] - p130 = p3dat[ind] - p130 = p130[p130 > 0] - if len(p56) > 0: - p56a = numpy.mean(p56) - p56m = min(p56) - else: - p56a = 0 - p56m = 0 + if len(dat) > 0: + # There are nominal quality data in the last 2 hours, + # so create a string, ace_table, with data for html display. + # However, some channels may still contain bad data + # indicated by -1.e5 values. These are going to be + # displayed but filtered out before computing the stats. - if len(p130) > 0: - p130a = numpy.mean(p130) - p130m = min(p56) - else: - p130a = 0 - p130m = 0 + # Add scaled P5 and P6 + dat['p3_scaled_p5'] = dat['p5'] * P3_P5_SCALE + dat['p3_scaled_p6'] = dat['p6'] * P3_P6_SCALE - e38 = de1[ind2] - e38 = e38[e38 > 0] - e175 = de4[ind2] - e175 = e175[e175 >0] + # ace_table with data for html display - if len(e38) > 0: - e38a = numpy.mean(e38) - e38m = min(e38) - else: - e38a = 0 - e38m = 0 - if len(e175) > 0: - e175a = numpy.mean(e175) - e175m = min(e175) - else: - e175a = 0 - e175m = 0 -# -#--- compute fluence of 2 hours -# - e38f = e38a * 7200. - e175f = e175a * 7200. - p56f = p56a * 7200. - p130f = p130a * 7200. - p337f = p337a * 7200. - p761f = p761a * 7200. - p1073f = p1073a * 7200. - p5_p3f = p5_p3a * 7200. - p6_p3f = p6_p3a * 7200. - p7_p3f = p7_p3a * 7200. -# -#--- compute spectral indcies -# - if p337a > 0: - p3_p5 = p130a / p337a - else: - p3_p5 = 0 - if p761a > 0: - p3_p6 = p130a / p761a - else: - p3_p6 = 0 - if p761a > 0: - p5_p6 = p337a / p761a - else: - p5_p6 = 0 - if p1073a > 0: - p6_p7 = p761a / p1073a - else: - p6_p7 = 0 - if p1073a > 0: - p5_p7 = p337a / p1073a - else: - p337a = 0 -# -#--- violation check -# - if p130f > 360000000: - val = "%.4e" % p130f - ace_violation_protons(val) + cols = ("yr", "mo", "da", "hhmm", "de1", "de4", + "p2", "p3", "p3_scaled_p5", "p3_scaled_p6", + "p5", "p6", "p7") + row_format = ' '.join(["{:s} " * 4, "{:11.3f} " * 9, "\n"]) - p5_p6_lim = 1.0 #----????? this is what given in the original, but use below - p5_p6_lim = 1.0e10 + for row in dat: + ace_table = ace_table + row_format.format(*row[cols]) - if (p5_p6 > p5_p6_lim) or (p5_p6 < 1): - speci = "%12.1f" % p5_p6 - speci_lim = "%8.1f" % p5_p6_lim + # Add a header for the table with stats - ace_invalid_spec(speci, speci_lim) -# -#--- crate a summary table -# - hfile = templ_dir + 'header2' - line = line + "\n" + read_file(hfile) + with open(f"{args.ace_dir}/Scripts/Template/header2") as fp: + header = fp.read() - line = line + "%7s %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f\n"\ - % ("AVERAGE ", e38a, e175a, p56a, p130a, p5_p3a, p6_p3a, p337a, p761a, p1073a) + ace_table = ace_table + "\n" + header - line = line + "%7s %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f\n"\ - % ("MINIMUM ", e38m, e175m, p56m, p130m, p5_p3m, p6_p3m, p337m, p761m, p1073m) + # Calculate electron/proton stats (mean, min, 2h fluence, spectra) - line = line + "%7s %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n\n"\ - % ("FLUENCE ", e38f, e175f, p56f, p130f, p5_p3f, p6_p3f, p337f, p761f, p1073f) + stats = calculate_stats(dat) - line = line + "%7s %11s %11.3f %11s %11.3f %11s %11.3f %11s %11.3f \n\n"\ - % ("SPECTRA ", "p3/p5", p3_p5, "p3/p6", p3_p6, "p5/p6", p5_p6, "p6/p7", p6_p7) + # Add table with stats for the display - line = line + "%62s %4.1f\n"\ - % ("* This P3 channel is currently scaled from P5 data. P3* = P5 X ", p5_p3_scale) + for stat, name, fmt in zip((stats['medians'], stats['mins'], stats['fluences_2h']), + ('AVERAGE', 'MINIMUM', 'FLUENCE'), + ('{:11.3f}', '{:11.3f}', '{:11.4e}')): + vals = stat.values() + row_format = ' '.join([f"{name.ljust(15)}", f"{fmt} " * len(vals), "\n"]) + ace_table = ace_table + row_format.format(*list(vals)) - line = line + "%62s %4.1f\n"\ - % ("** This P3 channel is currently scaled from P6 data. P3** = P6 X ", p6_p3_scale) + # Add info on spectral indexes - line = line + "%62s %4.1f\n"\ - % ("*** This P3 channel (not shown) is currently scaled from P7 data. P3*** = P7 X ", p7_p3_scale) + ace_table = ace_table + "\n" + "%7s %11s %11.3f %11s %11.3f %11s %11.3f %11s %11.3f \n\n"\ + % ("SPECTRA ", "p3/p5", stats['spectra']['p3_p5'], + "p3/p6", stats['spectra']['p3_p6'], + "p5/p6", stats['spectra']['p5_p6'], + "p6/p7", stats['spectra']['p6_p7']) - return line + ace_table = ace_table + "%62s %4.1f\n"\ + % ("* This P3 channel is currently scaled from P5 data. P3* = P5 X ", + P3_P5_SCALE) -#--------------------------------------------------------------------------------------------------- -#-- ace_violation_protons: send out a warning email -- -#--------------------------------------------------------------------------------------------------- + ace_table = ace_table + "%62s %4.1f\n"\ + % ("** This P3 channel is currently scaled from P6 data. P3** = P6 X ", + P3_P6_SCALE) -def ace_violation_protons(val): + ace_table = ace_table + "%62s %4.1f\n"\ + % ("*** This P3 channel (not shown) is currently scaled from P7 data. P3*** = P7 X ", + P3_P7_SCALE) + + else: + # No valid data in the last 2h + ace_table = """ ++ No valid data for last 2 hours. +
+""" + + # download images and reverse the color: + download_img(img=HTTP_EPAM, + html_dir=args.html_dir, + ace_dir=args.ace_dir) + download_img(img=HTTP_MAG, + html_dir=args.html_dir, + ace_dir=args.ace_dir) + + # read in headers, footers and image html from template files + # and insert ace_table string with current data and stats + contents = [] + files = ['header', 'header1', 'image2', 'footer'] + for file in files: + with open(f"{args.ace_dir}/Scripts/Template/{file}") as fp: + contents.append(fp.read()) + + contents.insert(2, ace_table) + + # update ACE html page + with open(f"{args.html_dir}/ACE/ace.html", "w") as fo: + fo.write("\n".join(contents)) + + +def read_2h_ace_data_from_12h_archive(fname): """ - send out a warning email - input: val --- value - output: email sent out + Reads 12h ACE archive and keeps only the last 2h + of nominal quality (both the e_status and p_status + flags equal 0). + + :param fname: ACE 12h local MTA archive + :returns: astropy Table with ACE data """ - line = 'A Radiation violation of P3 (130KeV) has been observed by ACE\n' - line = line + 'Observed = ' + val + '\n' - line = line + '(limit = fluence of 3.6e8 particles/cm2-ster-MeV within 2 hours)\n' - line = line + 'see http://cxc.harvard.edu/mta/ace.html\n' -# -#--- add current spacecraft info -# - line = line + curr_state() -# -#--- read alert case -# - if os.path.isfile('/data/mta4/www/Snapshot/.scs107alert'): - line = line + 'The ACIS on-call person should review the data and call a telecon if necessary.\n' - line = line + 'This message sent to sot_ace_alert\n' - send_mail('ACE_p3', line, 'sot_ace_alert@cfa.harvard.edu') -# -#--- yellow alert case -# - else: - line = line + 'This message sent to sot_yellow_alert\n' - send_mail('ACE_p3', line, 'sot_yellow_alert@cfa.harvard.edu') -#--------------------------------------------------------------------------------------------------- -#-- ace_invalid_spec: sending out a warning email -- -#--------------------------------------------------------------------------------------------------- + # TODO: update scripts that create the 12h archive to + # include the column names + names = ("yr", "mo", "da", "hhmm", "mjd", "mjdsec", + "e_status", "de1", "de4", + "p_status", "p2", "p3", "p5", "p6", "p7", "anis_index") -def ace_invalid_spec(speci, speci_lim): - """ - sending out a warning email - input: speci --- the value reported - speci_lim --- the limit of the vluae - output: email sent out - """ -# -#--- check whether the mail is recently sent out -# - out = '/tmp/mta/prot_spec_violate' - if os.path.isfile(out): - cmd = 'date >> ' + out - os.system(cmd) -# -#--- if not, send the warning email -# - else: - line = " A spectral index violation of P5/P6 has been observed by ACE, " - line = line + "indicating a possibly invalid P5 channel.\n" - line = line + 'Observed = ' + speci - line = line + '(limit = ' + speci_lim + ')\n' - line = line + 'see http://cxc.harvard.edu/mta/ace.html\n' - line = line + 'This message sent to mtadude\n' - - send_mail('ACE_p5/p6', line, 'mtadude@cfa.harvard.edu') -# -#--- open a file indicating that the mail was sent -# - with open(out, 'w') as fo: - fo.write(line) + formats = ('S4', 'S2', 'S2', 'S4', 'f8', 'f8', + 'i1', 'f8', 'f8', + 'i1', 'f8', 'f8', 'f8', 'f8', 'f8', 'i1') -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- + dat = np.loadtxt(fname, dtype={'names': names, 'formats' : formats}) + dat = Table(dat) -def send_mail(subject, content, address): + # remove duplicate lines; TODO check why they appear + dat = unique(dat, keys=("yr", "mo", "da", "hhmm")) - fo = open(zspace, 'w') - fo.write(content) - fo.close() + fits_fmt = [f"{yr}-{mo.zfill(2)}-{da}T{hhmm[:2]}:{hhmm[-2:]}:00.000" + for yr, mo, da, hhmm in zip(dat['yr'], dat['mo'], + dat['da'], dat['hhmm'])] - cmd = 'cat ' + zspace + '|mailx -s "' + subject + '" ' + address - os.system(cmd) - - cmd = 'rm ' + zspace - os.system(cmd) + tt = DateTime(np.array(fits_fmt), format='fits').date -#--------------------------------------------------------------------------------------------------- -#-- curr_state: extract some satellite related information -- -#--------------------------------------------------------------------------------------------------- + # keep only the most recent 2h of nominal quality data + # subtract 2.25h to account for cronjob processing every 5 min + now_minus_2h = DateTime(DateTime() - 2.25 / 24).date + + # status: 0 = nominal, 4,6,7,8 = bad data, unable to process, 9 = no data + ok = (dat['e_status'] == 0) & (dat['p_status'] == 0) & (tt > now_minus_2h) + dat = dat[ok] + + return dat -def curr_state(): + +def calculate_stats(dat): """ - extract some satellite related information - input: none - output: line --- lines which display the satellite information + Calculate statistics: medians, mins, 2h fluences in each channel, + and spectral indexes. + :param dat: astropy Table with ACE data + :returns stats: dictionary with values being dictionaries of stats + for ACE channels """ -# -#--- read CRM file to get some information -# - sfile = crm_dir + 'CRMsummary.dat' - data = mcf.read_data_file(sfile) - - dirct = 'NA' - dist = 'NA' - if len(data) < 1: - dline = '' - else: - dline = data[0] - for ent in data: - mc = re.search('Geocentric Distance', ent) - if mc is not None: - atemp = re.split(':', ent) - btemp = re.split('\s+', atemp[1].strip()) - dist = btemp[0] - if btemp[1] == 'A': - dirct = 'Ascending' - else: - dirct = 'Descending' - - line = 'Currently ' + dirct + ' through ' + dist + 'km with ' + dline + '\n\n' -# -#-- get DSN contact info from ephem data -# - ctime = time.strftime("%Y:%j:%H:%M:%S", time.gmtime()) - start = Chandra.Time.DateTime(ctime).secs - stop = start + 3.0 * 86400. - - ifile = ephem_dir + 'dsn_summary.dat' - data = mcf.read_data_file(ifile) - data = data[2:] - line = line + data[0] + '\n' + data[1] + '\n' - for ent in data: - atemp = re.split('\s+', ent) - year = atemp[10] - yday = atemp[11] - stime = convert_to_stime(year, yday) - if stime >= start and stime < stop: - line = line + ent + '\n' - - line = line + '\n' - - return line -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- + stats = {'medians': {}, + 'mins': {}, + 'fluences_2h': {}, + 'spectra': {'p3_p5': -1e5, 'p3_p6': -1e5, + 'p5_p6': -1e5, 'p6_p7': -1e5}} + + cols = ("de1", "de4", "p2", "p3", "p3_scaled_p5", "p3_scaled_p5", + "p5", "p6", "p7") + + for col in cols: + # Filter out negative (bad) values + ok = dat[col] > 0 + if len(dat[col][ok]) > 2: + # At least 2 measurements of nominal quality + # in the last 2 hours + stats['medians'][col] = np.median(dat[col][ok]) + stats['mins'][col] = np.min(dat[col][ok]) + stats['fluences_2h'][col] = np.median(dat[col][ok]) * 7200 + else: + stats['medians'][col] = -1e5 + stats['mins'][col] = -1e5 + stats['fluences_2h'][col] = -1e5 -def convert_to_stime(year, yday): + # Spectral indexes + meds = stats['medians'] + if all(x > 0 for x in [meds['p3'], meds['p5']]): + stats['spectra']['p3_p5'] = meds['p3'] / meds['p5'] - atemp = re.split('\.', yday) - frac = float(atemp[1]) - val = 24. * frac - hh = int(val) - diff = val - hh - val = 60 * diff - mm = int(val) - diff = val - mm - ss = int(60 *diff) + if all(x > 0 for x in [meds['p3'], meds['p6']]): + stats['spectra']['p3_p6'] = meds['p3'] / meds['p6'] - htime = str(year) + ':' + atemp[0] + ':' + str(hh) + ':' + str(mm) + ':' + str(ss) - - stime = Chandra.Time.DateTime(htime).secs + if all(x > 0 for x in [meds['p5'], meds['p6']]): + stats['spectra']['p5_p6'] = meds['p5'] / meds['p6'] - return stime + if all(x > 0 for x in [meds['p6'], meds['p7']]): + stats['spectra']['p6_p7'] = meds['p6'] / meds['p7'] + + return stats -#--------------------------------------------------------------------------------------------------- -#-- download_img: down load an image from web site -- -#--------------------------------------------------------------------------------------------------- -def download_img(file, chg=1): +def download_img(img, html_dir, ace_dir, chg=1): """ - down load an image from web site - input: file --- image file address - chg --- if >0, reverse the color - output: