Plot TAI-weights from BIPM Circular-T data file

tai_weight_plot.py

"""
    Read a clock-weight file from the BIPM Circular-T website
    plot some figres
    
    AW2021-03-28, update 2021-03-29
"""

import os
import datetime

#import mjdutil
import matplotlib.pyplot as plt
                
# BIPM clock types
# https://webtai.bipm.org/database/clock.html
type_dict = { 18: "18 MICROSEMI Cs 4000",
              35: "35 MICROSEMI 5071A HIGH PERFORMANCE TUBE.",
              36: "36 MICROSEMI 5071A LOW/STANDARD PERFORMANCE TUBE",
              40: "40 UNSPECIFIED HYDROGEN MASER",
              41: "41 HYDROGEN MASER",
              92: "92 GROUND-STATE HYPERFINE TRANSITION OF 133 Cs",
              93: "93 GROUND-STATE HYPERFINETRANSITION OF 87 Rb", 
            }

def read_weight(directory=os.getcwd(), prefix="w", dt=datetime.datetime(2017,12,01)):
    fname = directory + "%s%02d.%02d"%(prefix, dt.year-2000, dt.month) # filename for year and month file e.g. "w17.12"
    with open(fname) as f:
        mjd_found=False
        mjds=[]
        weights_type={} # dict for weight per clock type
        weights_lab={} # dict for weight per lab
        n_clocks = 0
        for line in f:
            line=line.split()
            if len(line) and line[0] == "LAB." and not mjd_found: # this repeats many times!
                for n in range(6):
                    mjds.append( int(line[2+n]) )
                mjd_found=True
                #print "FOUND mjds",mjds

            clk_type = 0
            try:
                clk_type=int(line[1])
            except:
                clk_type=0
                
            if len(line) >= 9 and clk_type in type_dict.keys():
                #clk_type = int(line[1])
                clk_id = int(line[2])
                lab = line[0]
      
                n=5 # read the last MJD date column
                w = line[3+n]
                if w=="*****" or w=="*********":
                    w=0
                else:
                    w = float(w)
                #print clk_type, clk_type, w
                n_clocks += 1
                if clk_type not in weights_type:
                    weights_type[clk_type]=0
                weights_type[clk_type] += w
                
                if lab not in weights_lab:
                    weights_lab[lab]=0
                weights_lab[lab] += w
    mjd0 = mjds[-1]
    return mjd0, weights_type, weights_lab, n_clocks



if __name__ == "__main__":
    directory = os.getcwd() + "/ClockData/" # where we find the w-file
    print "looking for w-file in: ", directory
    dt = datetime.datetime(2021,01,1) # year and month (day doesn't matter)
    m, w, l, n = read_weight(directory, 'w', dt)
    
    print "Weight-file for ", dt.year, ".", dt.month
    print "number of clocks ", n
    print "MJD of last column ", m
    print "Number of types: ", len(w)
    print "Sum of weights per type: ", sum(w.values())
    print "Type dict ", w
    print "Number of labs ", len(l)
    print "Sum of weight per lab: ", sum(l.values())

    # Plot figures
    fig1, ax1 = plt.subplots()
    plt.pie( w.values(), labels= [type_dict[ k ] for k in w.keys()], autopct='%1.1f%%' , rotatelabels=True)
    ax1.axis('equal') 
    plt.title('TAI weight per clock type for %d.%d' % (dt.year, dt.month))

    fig1, ax1 = plt.subplots()
    plt.pie( l.values(), labels = l.keys(), autopct='%1.1f%%', rotatelabels=True )
    ax1.axis('equal') 
    plt.title('TAI weight per laboratory for %d.%d' % (dt.year, dt.month))    
    
    plt.show()

"""
    Output for w21.01
    
    Weight-file for  2021 . 1
    number of clocks  422
    MJD of last column  59244
    Number of types:  7
    Sum of weights per type:  99.984
    Type dict  {35: 7.026999999999994, 36: 0.2800000000000001, 40: 75.53500000000004, 41: 12.853, 18: 0.008, 92: 0.221, 93: 4.06}
    Number of labs  66
    Sum of weight per lab:  99.984


"""