# -*- coding: utf-8 -*-
"""
TA_gen2_auto12A.py:
Simulation tool of timescale generation

- Input datafiles are read from the list './clock_file_list.txt'.
- Input clock data are prepared in the '../Data_Clock".
- Output of menu 0,1,2 are saved in the './TA_RESULT'.
- Output of menu 2 is saved in the './TA_GEN'.

Created on 22 May, 2025
@author: Y. Hanado
"""

import numpy as np
import matplotlib.pyplot as plt
import datetime
import allantools
import os
import shutil

# ==============================================================
# DEF : calculation ADEV by "allantools"
def adev_calc(data_list, menu_n, x_limit):                                         

    line = []
    name = []
    tt = []                                                                        
    xx = []                                                                        
    taus = []
    adevs = []
    data_adev = []
    data_adev_list = []
    MJDspan = []                                                                   

  # search the MJD span of calclated TA   
    for nclk, line in enumerate(data_list):
        if (menu_n =='2'):
            name = line[0]
            if ('TA' in name):
                tt = line[1][0]
                xx = line[1][1]  
                t = []
                x = []
                for i in range(len(tt)):
                    if ((xx[i] < x_limit) & (np.isnan(xx[i])!=True)):                    
                        t.append(tt[i])
                        x.append(xx[i])
                TA_MJD0 = int(t[0])              
                TA_MJD1 = int(t[-1])             

    name = []
    tt = []                                
    xx = []  
  # calculate ADEV for each data of data_list
    for nclk, line in enumerate(data_list):
        name = line[0]
        tt = line[1][0]
        xx = line[1][1]
        
      # if the plot is for TA calculation.
        t = []
        x = []
        if (menu_n =='2'):
            for i in range(len(tt)):
                if (tt[i] >= TA_MJD0) and (tt[i] <= TA_MJD1):
                    t.append(tt[i])                    
                    x.append(xx[i])                    
            tt = t
            xx = x
                    
      # fix MJD span for ADEV  
        t = []
        x = []
        p_data = []
        E_flag = False
        for i in range(len(tt)):
          # Separate data to portion by before/after missing data   
            if ((xx[i] < x_limit) & (np.isnan(xx[i])!=True)):                           
                t.append(tt[i])
                x.append(xx[i])
                E_flag = True            
            else:
                if (E_flag == True):
                    MJD0 = int(t[0])              
                    MJD1 = int(t[-1])             
                    MJDspan = [MJD0, MJD1]        
                    p_data.append([t, x, MJDspan])
                    E_flag = False
                    t = []
                    x = []
                    
        if (len(t) != 0):     
            MJD0 = int(t[0])              
            MJD1 = int(t[-1])             
            MJDspan = [MJD0, MJD1]        
        # p_data is the portion data list deparated by the data missing 
            p_data.append([t, x, MJDspan])
       
      # Adev calclation for each portion data 
        x = []                                                                    
        t = []                                                                    
        for i in range(len(p_data)):                                                    
            if (i > 0):                                                           
                name = name + "*" + str(i)                                        
            t = p_data[i][0]                                                      
            x = p_data[i][1]                                                       
            MJDspan = p_data[i][2]                                                
            t_int = float(t[1] - t[0])*86400   
          # np.nan is in x?                                                         
            npx=np.isnan(x)                                                         
            if ((any(xx >= x_limit for xx in x) != True) 
                & (npx.any() != True)):                                               
              # Adev by allantools
                (taus, adevs, errors, ns) = allantools.oadev(x)
                taus = taus*t_int
                adevs = adevs/t_int
                data_adev = np.array([taus, adevs])
              # data_adev_list[] is [[CS**, [tau, ADEV]],,,]
                data_adev_list.append([name, data_adev, MJDspan])                  
            else:
                print("!!!!! " + name + "  is not available for allantools.")

    return data_adev_list

# ==============================================================
# DEF : plot ADEV results 
def adev_plot(outfile,wgtfile,all_wgt_list,MJDall,log,comment,menu_n,x_limit):                  

  # data_list[] is the data from outfile <== "all_*.dat"
    data_list = []
    data_list = read_jointdata(outfile, MJDall, x_limit)                           

  # pickup only the selected clock if menu_n=2                                    
    data_list_ForAdev = []
    if (menu_n == '2'):
        for nclk, line in enumerate(data_list):
            if (line[-3] != 0 or line[-2] == 3):                                        
                data_list_ForAdev.append(line)
        data_list = data_list_ForAdev
  
  # calculate ADEV by allantools
  # data_adev_list[] is [["UTCk-CSxx", [[tau,,,],[ADEV,,,]]],,,]
    data_adev_list = []
    data_adev_list = adev_calc(data_list, menu_n, x_limit)                         

  # Error message
    if (data_adev_list == []):
        print("!!!!! There are no clock for ADEV calcluation !")
        return
  
  # format Adev data-set for plotting   ............................               
    tau_list = data_adev_list[0][1][0]
    for nclk in range(1,len(data_adev_list)):
        tau_nclk = data_adev_list[nclk][1][0]
        if (max(tau_nclk) > max(tau_list)):
            tau_list = np.append(tau_list, max(tau_nclk))                          
    for n in range(1,len(data_adev_list)):
        tau_clk= data_adev_list[n][1][0]
        adv_clk= data_adev_list[n][1][1]
        if (max(tau_clk) < max(tau_list)):
            tau_clk = np.append(tau_clk, max(tau_list))
            adv_clk = np.append(adv_clk, np.nan)
        data_adev_list[n][1] = [tau_clk,adv_clk]
  
  # plot graph
    fig = plot_data(data_list, data_adev_list, menu_n)
    plt.show(block=False)

  # plot `all_weight.dat' 
    w_fig = ""                                                      
    if (menu_n == '2'):  
        w_fig = plot_weight(all_wgt_list)
        plt.show(block=False)

  # save? or not?
    print ('\n' + comment)                                                        
    while (save_yn := input ("????? Save figure? (y/n) = ")) not in ('y','n'):                   
        print("!!!! Input is not correct. Please try again!" + '\n')                             
    if (save_yn == "y"):
        log, fig = save_logfigdat(outfile, data_adev_list, wgtfile, 
                                  log, fig, w_fig)                                 
    return

# ==============================================================
# DEF : read initial weight and flag
def clk_init_wgt(infile_clklist):
    
  # read initial weight from 'clock_file_list.txt'                                 
    wgt_1 = []
    with open(infile_clklist, "r") as infile:
        for line in infile:
            if (line[0] != "#"):
                name = line.split()[1]                                             
                name = name.split('.')[0]
                name = name.split('-')[1]                                            
                wgt = float(line.split()[3])                                       
              # if flg=1 then weight=0      
                wgt_1.append([name, wgt])
    
    return wgt_1

# ==============================================================
# DEF: joint single column clock data 
def joint_rawdata(outfile, MJDall, data_list, log):
 
  # set MJDall[] as the first element of all_data[]
    all_data = []
    all_data.append(MJDall)     # all_data = [[56xx0,56xx1,,,]]
    output_line0 = ["fname "]                                                      
    output_line1 = ["weight"]
    output_line2 = ["flag  "]
    output_line3 = ["cid   "]                                                      
    
  # read clock data from data_list[]
    for nclk, line in enumerate(data_list):
        name = line[0]          # name is "UTCk-CSxx"
        t = line[1][0]          # t[] is list of MJD
        x = line[1][1]          # x[] is list of UTCk_CS data
        wgt = line[2]           # wgt is a value of weight (not list)
        flg = line[3]           # flg is a value (not list)
        cid = line[4]           # cid is a value of clock id (not list)                  
        
        output_line0.append(" " + name + " ")
        output_line1.append("      " + str(wgt) + "      ")
        output_line2.append("       " + str(flg) + "       ")
        output_line3.append("       " + str(cid) + "       ")                      
        
      # Missing data => np.nan                                                     
        xx=[]
        for i in range(len(MJDall)):
            aruyo = (MJDall[i] in t)
            if (aruyo == True): 
                for j in range(len(t)):
                    if (t[j] == MJDall[i]): 
                        xx.append(x[j])
                        break
            else: 
                xx.append(np.nan)
                
      # all_data[] is a joint of each column of clock data
      # all_data=[[MJD,,,],[c#1data,,],[C#2data,,,],,,]
        col_clock = []
        for i in range(len(MJDall)):
            col_clock = np.hstack((col_clock,xx[i])) # col_clock is a column of UTCk-CS data
        all_data.append(col_clock)                   # all_data is multi-columns of UTCk-CS data
        
  # save the joint data file ==> "temp.dat"
    n_lin = len(col_clock)
    n_col = len(all_data)
    output = np.transpose(all_data)
    np.savetxt("temp.dat", output.reshape((n_lin,n_col)),
               fmt='  '.join(['%i'] + ['%+.6e']*(n_col-1)))

  # save the joint data file with title & weight & flag ==> outfile
    with open(outfile, "w") as outf:
        outf.writelines(output_line0 + ['\n'])
        outf.writelines(output_line1 + ['\n'])
        outf.writelines(output_line2 + ['\n'])
        outf.writelines(output_line3 + ['\n'])                                     
        with open("temp.dat", "r") as inf:
            for line in inf:
                outf.writelines(line)

    return outfile, log

# ==============================================================
# DEF: assign the MJD range within available one  
def MJD_limit(MJDall, log):
    
  # input the param of MJD range as "MJD1 MJD2" from the console 
    chk1 = False
    while chk1 == False:
        try:
            comment = ("\n**1** Assign the MJD range for calculation in [" 
                       + str(MJDall[0]) + "," + str(MJDall[-1]) + '].\n')
            MJD1, MJD2 = map(int, input(comment + "\n??1?? Input MJD1 "
                                        "and MJD2 with spacing = ").split()) 
        except ValueError:
            print('!!! Input ValueError, try again !!!')
        else:
            print('\n**1** Selected range is : [', MJD1, ' - ', MJD2, ']')
            if MJD1 < MJDall[0]:
                print("!!!!! MJD1 should be equal or larger than ", 
                      MJDall[0], ". !!!!!")                                        
            elif MJD2 > MJDall[-1]:
                print("!!!!! MJD2 should be equal or smaller than ", 
                      MJDall[-1], ". !!!!!")                                       
            elif MJD1 > MJD2:
                print("!!!!! MJD1 should be smaller than MJD2. !!!!!")                
            else:
                chk1 = True

  # MJDall[] is [MJD1, MJD2]  
    MJDall = [MJDall[i] for i in range(len(MJDall)) 
        if (MJDall[i] >= MJD1) & (MJDall[i] <= MJD2)]
    MJDall = np.array(MJDall)

  # comment to log[]n
    comment = "**1** Assign the MJD range for calculation *****" + '\n'
    log_0 = ("::::: Selected MJD     :  [" + str(MJD1) + " - " + str(MJD2) 
             + ']\n')
    log = log + [comment + log_0]

    return MJDall, log

# ==============================================================
# DEF : maximum MJD for all selected clock file
def MJD_max(data_list):

  # search the Maximum MJD range over the data_list
    for nclk, line in enumerate(data_list):
        MJD_min_0 = int(np.amin(line[1][0]))
        MJD_max_0 = int(np.amax(line[1][0]))
        if (nclk == 0):
            MJD_min = MJD_min_0
            MJD_max = MJD_max_0
        if (MJD_min_0 <= MJD_min):
            MJD_min = MJD_min_0
        if (MJD_max_0 >= MJD_max):
            MJD_max = MJD_max_0
        
  # make MJDall[] in range of (MJD_min, MJD_max+1)
    MJDall = np.array([i for i in range(MJD_min, MJD_max+1, 1)])

    return MJDall

# ==============================================================
# DEF : plot the original & Allan deviation
def plot_data(data_list1, data_list2, menu_n):

  # set the plotting area
    fig = plt.figure(figsize=(12, 7))   

  # left box position is set as [left-x, ground-y, x-width, ywidth]
    rect = [0.05, 0.1, 0.35, 0.5] #0307
    xlabel = 'measurement time (MJD)'
    ylabel = 'time difference (sec)'
    ax1 = plt.axes(rect)
    ax1 = plot_set(ax1, xlabel, ylabel)

  # right box position is set as [left-x, ground-y, x-width, ywidth]
    rect = [0.48, 0.1, 0.35, 0.5] #0307
    xlabel = 'averaging time (sec)'
    ylabel = 'allan deviation'
    ax2 = plt.axes(rect)
    ax2 = plot_set(ax2, xlabel, ylabel)

  # set the colortable                                                            
    colors = plt.cm.tab20(np.linspace(0, 1, 20))                                  
  
  # set the fig-param in left box
    clk_param = []                                                                        
    for i in range(len(data_list1)):
        name1 = data_list1[i][0]
        data1 = data_list1[i][1]
        cwgt1 = data_list1[i][2]                                                          
        cflg1 = data_list1[i][3]                                                          
        c_color1 = colors[data_list1[i][4]]                                       
        if ("TA" in name1): c_color1 = "black"                                    
        if (cflg1 == 1):                                                          
            refclk = name1.split('-')[1].split('.')[0]                             
        clk_param.append([name1,c_color1,cwgt1, cflg1])                                   
        figx = data1[0]
        figy = data1[1]
        figy = [None if k == 1000 else k for k in figy]
        ax1.set_xlim(xmin=data1[0][0], xmax=data1[0][-1])
        ax1.plot(figx, figy, color=c_color1, marker="o", label=name1)
        ax1.legend(bbox_to_anchor=(0.28, 1.05, 0.38, 0.2), loc='center', ncol=3)

    if (menu_n == 0): 
        ax１.text(0.01, 0.01, f'${refclk}$ = reference clock', fontsize=10, 
                 color='dimgray',transform=ax1.transAxes) 
        
  # set the fig-param in right box
    for i in range(len(data_list2)):
        name2 = data_list2[i][0]
        data2 = data_list2[i][1]
        for j in range(len(clk_param)):                                                   
          # if multi ADEV in one clock                                            
            if ('*' in name2):                                                    
                name2a = name2.split('*') [0]                                     
                name2b = name2.split('*') [1]                                     
            else:                                                                 
                name2a = name2                                                    
                name2b = 0                                                        
            if (name2a == clk_param[j][0]):                                       
                c_color2 = clk_param[j][1]                                                
            if ("TA" in name2): c_color2 = "black"                                
        figx = data2[0]
        figy = data2[1]
        ax2.set_xlim(xmin=data2[0][0], xmax=data2[0][-1])
        ax2.loglog(figx, figy, color=c_color2, marker=f'${name2b}$', label=name2)                
        ax2.legend(bbox_to_anchor=(0.28, 1.05, 0.38, 0.2), loc='center', ncol=3)

    if (menu_n == '2'):                                                        
        ax2.text(0.01, 0.01, 'All ADEVs are calculated within the MJDs' 
                 +' of TA.', fontsize=9, color='red',transform=ax2.transAxes) 
      
    return fig

# ==============================================================
# DEF : set the basic plotting parameter
def plot_set(ax, xlabel, ylabel):

  # setting for the figure parameters
    ax.grid(which='both', color='gray', linestyle='dotted')
    ax.set_xlabel('x: ' + xlabel)
    ax.set_ylabel('y: ' + ylabel)
    ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0))

    return ax

# ==============================================================                   
# DEF : plot the weights 
def plot_weight(all_wgt_list):

  # set the plotting area
    fig = plt.figure(figsize=(12, 7))   

  # left box position is set as [left-x, ground-y, x-width, ywidth]
    rect = [0.05, 0.1, 0.35, 0.5] #0307
    xlabel = 'time (MJD)'
    ylabel = 'weight'
    ax1 = plt.axes(rect)
    ax1 = plot_set(ax1, xlabel, ylabel)

  # set the colortable                                             
    colors = plt.cm.tab20(np.linspace(0, 1, 20))                   
  
  # reconstract the lists from all_wgt_list
  # all_wgt_list=[[mjd1,[[clk,wgt,cid][clk,wgt,cid]..]],[mjd2,[][]..],]
    mjd = [t[0] for t in all_wgt_list]                              
    wgts = []                                                       
    for n in range(len(all_wgt_list[0][1])):                        
        xx = [x[1][n][1] for x in all_wgt_list]  #xx = [each clk's weights] 
        wgts.append(xx)                          #wgts = [[clk1's weights],,,]             
    data = [mjd] + wgts                                             
    cdat = [[c[0],c[2]] for c in all_wgt_list[0][1]]  #c[0]=clk, c[2]=cid 
  
  # set the fig-param in left box
    for i in range(len(cdat)):                                      
        name = cdat[i][0]                                           
        cid  = cdat[i][1]                                          
        c_color1 = colors[cid]                                     
        figx = data[0]
        figy = data[i+1]
        ax1.set_xlim(xmin=figx[0], xmax=figx[-1])
        ax1.plot(figx, figy, color=c_color1, marker=".", label=name)
        ax1.legend(bbox_to_anchor=(0.28, 1.01, 0.38, 0.2), loc='center', ncol=4)
        
    return fig

# ==============================================================
# DEF : read jointed all clock data
def read_jointdata(infile_all, MJDall, x_limit):                                   

  # data_org[] is the data of infile_all <== "all_*_clk.dat"
    data = []
    data_list = []
    data_org = np.genfromtxt(infile_all, dtype=str, invalid_raise=False)

  # data_org[] is the data of infile_all
    for ncol in range(len(data_org[0])):  # data_org[0] is ['MJD' 'UTCk-CS**.txt' ,,,]
        name = data_org[0][ncol]
        if (ncol == 0):                   # in the case of MJD
            wgt = 0
            flg = 0
            cid = 0                                                                
            t = [float(data_org[n][ncol]) for n in range(4, len(data_org))]
        else:                             # in the case of clock data
            wgt = float(data_org[1][ncol])
            flg = int(data_org[2][ncol])
            cid = int(data_org[3][ncol])                                           
            x = [float(data_org[n][ncol]) for n in range(4, len(data_org))]                
            data = np.array([t, x])
          # data_list[] is [[name, [MJD, UTCk-CS**], weight, flag],[],,,]
            data_list.append([name, data, wgt, flg, cid])                          
          # message of phase jump
            if (any(xx >= x_limit for xx in x) == True):                              
                print("\n" + "!!!!! " + name + "  includes outliers")

    return data_list

# ==============================================================                      
# DEF : read parameter file
def read_param_file():
    
    with open('param_list.txt', "r") as infile:
        for line in infile:
            param = line.split('#')[0]        
            param_nam = param.split('=')[0]            
            param_val = param.split('=')[1]
            param_nam = param_nam.strip()
            param_val = param_val.strip()
            if (param_nam == 't_intvl'): t_intvl = int(param_val)
            if (param_nam == 'x_limit'): x_limit = float(param_val)
            if (param_nam == 'Twgt_N') : Twgt_N  = int(param_val)
            if (param_nam == 'wgt_lmt'): wgt_lmt = float(param_val)
    
    return t_intvl, x_limit, Twgt_N, wgt_lmt

# ==============================================================
# DEF : read raw clock data
def read_rawdata(t_intvl, x_limit):                                                   

  # read the clockdata file name from 'clock_file_list.txt'        
    in_clocklist = 'clock_file_list.txt'   # list of clock data                    
    clock_list = []
    with open(in_clocklist, "r") as infile:
        for line in infile:
            if (line[0] != "#"):
                cnum = int(line.split()[0])                                        
                name = line.split()[1]                                             
                wgt = float(line.split()[3])                                       
                flg = int(line.split()[5])                                         
              # if flg=1 then weight=0                                                    
                clock_list.append([name, wgt, flg, cnum])                            

  # If too-many or no reference clock => error message
    flg1_clk_count = 0 
    for i in range(len(clock_list)):
        if (clock_list[i][-2] == 1):                                               
            flg1_clk_count = flg1_clk_count + 1
    if (flg1_clk_count == 0):
        print("!!!!! No reference clock was assigned in clock_file_list.txt => ")
        print("      Set the 1st clock as the reference clock.")
        clock_list[0][-2] = 1                                                      
        print("      ",clock_list[0],'\n')        
    if (flg1_clk_count > 1):
        print("!!!!! Too many reference clock is assigned in clock_file_list.txt!" + '\n')
        print("      Set the 1st clock as the reference clock.")
        clock_list[0][-2] = 1                                                      
        print(clock_list[0],'\n')    
  
  # initialize
    data = []
    t_raw = []
    x_raw = []
    namelist = []
    data_raw = []
    data_raw_list = []

  # read clock data file from ../Data_Clock/****.txt (<== clock_file_list.txt)
    for line in clock_list:
        clk_name = "../Data_Clock/" + line[0]
        namelist.append(line[0])
        clk_wgt = line[1]
        clk_flg = line[2]
        clk_id  = line[3]                                                          
        data = np.genfromtxt(clk_name, usecols=(0, 1), invalid_raise=False)
        t_raw = [data[n][0] for n in range(len(data))]
        x_raw = [data[n][1] for n in range(len(data))]
      # data_raw[] is [[MJD,,,],[phase,,,]]
        data_raw = np.array([t_raw, x_raw])
      # data_raw_list[] is [["UTCk-CSxx", [[MJD,,,],[phase,,,]], wgt, flg, id], [],,]
        data_raw_list.append([line[0], data_raw, clk_wgt, clk_flg, clk_id])                
        
  # comment to log[]n
    comment = ('**0** Read and plot initial clock data and Allan deviation' 
               '*****' + '\n')
    log_0 = "::::: Input clocklist   : " + in_clocklist + '\n'
    log_1 = "::::: Input clock file  : " + str(namelist)+ '\n'
    log_2 = "::::: data time interval (t_intval): " + str(t_intvl) + '(sec)\n'                   
    log_3 = "::::: Outlier limit (x_limit))     : " + str(x_limit) + '\n'                        
    log = [comment + log_0 + log_1 + log_2 + log_3]                                   

    return data_raw_list, log, comment

# ==============================================================                   
# DEF : read weight data
def read_wgt_file(wgtfile):

    with open(wgtfile, "r") as infile:
        w_clk = []
        for line in infile:
            if (line.startswith('#')):
                clk_list = line.split() 
                clk_num = len(clk_list)
                for i in range(2,clk_num):
                    w_clk.append(clk_list[i])
        data = np.genfromtxt(wgtfile, delimiter='', skip_header=1)
        w_data = data.T

    return w_data, w_clk

# ==============================================================                      
# DEF : renew the weights to be less than the maximum limit         
def re_weight_under_limit(wgt_lmt,wgt1):                                              
    
    wgt3 = wgt1                                                                       
    while any(x > wgt_lmt for x in wgt3):                                             
        wgt2 = [wgt_lmt if x > wgt_lmt else x for x in wgt3]                          
        wsum2 = sum(x for x in wgt2)                                                  
        wgt3  = [x/wsum2 for x in wgt2]                                               
    wgt3  = [x for x in wgt3]                                                         

    return wgt3                                                                       

# ==============================================================
# DEF: joint each column of clock data 
def save_jointdata(outfile, data_list, t_title, log):

  # title of the time as "Tau(s)"
    output_line0 = ["   " + str(t_title) + "  "]
    output_line1 = ["  MJDspan   "]            #***** yuko0607 

  # ...............................................................                
  # First culumn of all_data[] is values of Tau
  # get the maximum column from the ADEVs of all clocks  
    all_data = []
    maxt = []
    for nclk, line in enumerate(data_list):
        if (len(line[1][0]) > len(maxt)):
            maxt = line[1][0]    
    all_data = [maxt]
  # ...............................................................                

  # read Adev & clock data from data_list[]
    for nclk, line in enumerate(data_list):
        name = line[0]                # name = "Tau(s)" or "UTCk-CSxx"
        x = line[1][1]                # x is list of UTCk_CS data or ADEV
        MJD0 = line[2][0]                                                           
        MJD1 = line[2][1]                                                           
        output_line0.append(name + " ")
        output_line1.append(str(MJD0) + "-" + str(MJD1) + "    ")                  
      # ............................................................              
      # np.transpose needs the same length of all columns
        lenx = len(maxt) - len(x)
        if (lenx > 0):
          for i in range(lenx):                   
            x = np.append(x,np.nan)
      # ............................................................              
      # all_data[] is [[tau], [Adev of clock1], [Adev of clock2],,]
        col_clock = []
        for i in range(len(maxt)):                                                
            col_clock = np.hstack((col_clock,x[i])) 
        all_data.append(col_clock)                  
  
  # save the joint data file ==> "temp2.dat"
    n_lin = len(col_clock)
    n_col = len(all_data)
    output = np.transpose(all_data)
    np.savetxt("temp2.dat", output.reshape((n_lin,n_col)),
               fmt='  '.join(['%9d'] + ['%+.6e']*(n_col-1)))

  # save the joint data file with title ==> "all_adev.dat"
    with open(outfile, "w") as outf:
        outf.writelines(output_line0 + ['\n'] + output_line1 + ['\n'])                
        with open("temp2.dat", "r") as inf:
            for line in inf:
                outf.writelines(line)

    return outfile, log

# ==============================================================
# DEF : save data & log
def save_logfigdat(alldata, data_adev_list, wgtfile, log, fig, w_fig):                 

  # make the directory
    dt_str, dt_com, path = set_dirpath("./TA_RESULT/")                             
    try:
        os.makedirs(path)
    except FileExistsError:
        print("Directry already exists.")
  
  # set the output filename
    output_log = path + "log_" + dt_str + ".txt"
    output_fig = path + "fig_" + dt_str + ".pdf"
    output_adv = path + "all_adev_" + dt_str + ".dat"
    output_dat = path + "" + alldata.split('.')[0] + "_" + dt_str + ".dat"  
    output_wgt = path + wgtfile.split('.')[0] + "_" + dt_str + ".dat"                 
    output_wfg = path + "wgt_fig_" + dt_str + ".pdf"                               
    
  # message of the output path and filename
    print("")
    print("::::: save logfile = ", output_log)
    print("::::: save figfile = ", output_fig)
    print("::::: save alldata = ", output_dat)
    print("::::: save ADEVdat = ", output_adv)

  # save logfile to './TA_RESULT/log_mmdd_hhmm.txt'
    np.savetxt(output_log, np.transpose(log), fmt="%s", 
               header = "Data log at " + dt_com + '\n')

  # save figure to './TA_RESULT/fig_mmdd_hhmm.pdf'
    fig.savefig(output_fig, bbox_inches='tight') 
    if (w_fig != ""):                                                               
        w_fig.savefig(output_wfg, bbox_inches='tight')                             
    
  # copy the alldata(='all_*.dat') to './TA_RESULT/all_*_mmdd_hhmm.dat' 
    shutil.copyfile(alldata, output_dat)

  # save the adev result in 'all_adev.dat'
    outfile, log = save_jointdata("all_adev.dat", data_adev_list, "Tau(s)", log)
  # copy the outfile(='all_adev.dat') to './TA_RESULT/all_adev_mmdd_hhmm.dat' 
    shutil.copyfile(outfile, output_adv)

  # copy the alldata(='all_*.dat') to './TA_RESULT/all_*_mmdd_hhmm.dat' 
    if (wgtfile != ""):                                                            
        print("::::: save weights = ", output_wgt)
        print("::::: save wgt_fig = ", output_wfg)                                   
        shutil.copyfile(wgtfile, output_wgt)                                       

    return log, fig

# ==============================================================
# DEF : merge weight 
def save_weight(all_wgt_list):
    
  # all_wgt_list = [MJD, [[clk1,w1],[clk2,w2],,,  ]]
    all_MJD = []
    wlist = {}
    for i in range(len(all_wgt_list)):
        MJD = all_wgt_list[i][0]
        c_w = all_wgt_list[i][1]
      # c_w = [[clk1,w1],[clk2,w2],,,  ]              
        for j in range(len(c_w)):
            clk = c_w[j][0]
            wgt = c_w[j][1]
            wgt = '{: 5.3f}'.format(wgt)
          # wgt is merged as one-list in each clock    
            if (MJD == all_wgt_list[0][0]): 
                wlist[f"{clk}"] = []
            wlist[f"{clk}"].append(wgt)

        MJD = '{:.1f}'.format(MJD)
        all_MJD.append(MJD)
    
  # all_data = [[MJD,,],[wgt1,,],[wgt2,,],,]    
    all_w_data = [all_MJD]
    line0 = "MJD  "
    for name, lst in wlist.items():
        line0 = line0 + "   " + name
        all_w_data.append(lst)
        
    wgtfile = 'all_weight.dat'
    np.savetxt(wgtfile, np.transpose(all_w_data), header=line0, fmt='%s') 
          
    return wgtfile

# ==============================================================
# DEF : make the data directry                                                     
def set_dirpath(dirname):
    
  # get now-time
    dt_now = str(datetime.datetime.now())
    dt_mmdd = dt_now[5:7] + dt_now[8:10]
    dt_hhmmss = dt_now[11:13] + dt_now[14:16] + dt_now[17:19]
    dt_str = str(dt_mmdd)+"_"+str(dt_hhmmss)
    dt_com = (dt_now[0:4] + "." + dt_now[5:7] + "." + dt_now[8:10] 
              + "_" + dt_now[11:13] + "h" + dt_now[14:16] + "m"
              + "_" + dt_now[17:19] + "s")
    path = dirname + dt_str + "/"                                                  
    
    return dt_str, dt_com, path

# ==============================================================
# DEF :  calculate a TAx prediction (= xhat) for next day
def TA_calc_x(TAx_list, TAxhat_list, diff_list, MJD, TAgen_dir, 
              TArate_span, TAwgt_span, wgt_addinfo, x_limit, Twgt_N, wgt_lmt): 
                                                                                     
    diff_data = []
    data = []
  # diff_data[] is all the data of diff_list[] 
    for nclk in range(len(diff_list)):
        name = diff_list[nclk][0]
        name = name.split('-')[0]
        t = diff_list[nclk][1]
        x = diff_list[nclk][2]
      # diff_date is the same as MJD? 
        for i in range(len(t)):
            if (t[i] == MJD):
                data = [name, float(x[i])]
        diff_data.append(data)

  # loop for making the sum of weight                                              
    wsum = 0
    wgt_list = []                         
    for nclk in range(len(TAxhat_list)):
        TA_name = TAxhat_list[nclk][0]
        TA_t = int(TAxhat_list[nclk][1])
        TA_xhat = TAxhat_list[nclk][4]
        wgt = TAxhat_list[nclk][5]                         
        flg = TAxhat_list[nclk][6]
        cid = TAxhat_list[nclk][7]                                               
      # clk0 also joins                                                            
        if (flg == 1):                                                            
            wsum = wsum + wgt                                                     
      # search the "clk0-clk" data on this MJD from diff_data[]
        for n in range(len(diff_data)):
            df_name = diff_data[n][0]
            df_x = diff_data[n][1]
          # wsum by available clocks for TA calculation
            if (TA_name == df_name) & (TA_t == MJD):
              # judge of missing-data or outlier
                if (((TA_xhat - df_x) < x_limit )                                     
                    and (np.isnan(TA_xhat - df_x) != True)):                          
                    wsum = wsum + wgt
                else:                                   
                    wgt = 0.0                          
        wgt_list.append([TA_name, wgt, cid])
        
  # weight normaization                                                              
    for i in range(len(wgt_list)) : 
        if (wsum != 0.0):
            wgt_list[i][1] = wgt_list[i][1]/wsum
        else:
            print(MJD, " !!!!! No available clock in averaging !!!!") 

  # renew the weight under the max_limit                                              
    if (wgt_lmt != 'OFF'):                                                            
        wgt1 = [x[1] for x in wgt_list]                                               
        wgt3 = re_weight_under_limit(wgt_lmt,wgt1)                                    
        wgt_list = [[a,b,c] for (a,_,c), b in zip(wgt_list, wgt3)]                     

  # loop for reading TAx & TAxhat from TAxhat_list[]
    TA_x0 = 0.0
    TA_x0_0 = 0.0

    for nclk in range(len(TAxhat_list)):
        TA_name = TAxhat_list[nclk][0]
        TA_t = int(TAxhat_list[nclk][1])
        TA_xhat = TAxhat_list[nclk][4]
        flg = TAxhat_list[nclk][6]
      # Search clk0  
        if (flg == 1):
            TA_name0 = TA_name
      # search weight                                                              
        for i in range(len(wgt_list)):                                             
            if (TA_name == wgt_list[i][0]):                                        
                wgt_n = wgt_list[i][1]                                             
      # search the "clk0-clk" data on this MJD from diff_data[]
        for n in range(len(diff_data)):
            df_name = diff_data[n][0]
            df_x = diff_data[n][1]
          # weighted average for TAx of clock0(flg=1)
            if (TA_name == df_name) & (TA_t == MJD):
              # judge of missing-data or outlier
              # if anomalies, TA is saved as formaer value
                if (np.isnan(TA_xhat - df_x) != True):                             
                    TA_x0 = TA_x0 + wgt_n * (TA_xhat - df_x)
                    TA_x0_0 = TA_x0 
                else:
                    TA_x0 = TA_x0_0                                                
      # clk0 also joins                                                           
        if (flg == 1):                                                            
            TA_x0 = TA_x0 + wgt_n * TA_xhat                                       

    # If there is no available clock             
    if (wsum == 0):    
        print(MJD, " !!!!! No available clock for TA !!!!") 
        TA_x0 = np.nan
  
  # save the TAx of the reference clock                    
    TA_data = [[TA_name0, MJD, TA_x0]]

  # calculate TAx of other clocks
    for n in range(len(diff_data)):
        df_name = diff_data[n][0]
        df_x = diff_data[n][1]
        TA_x = float(df_x) + TA_x0
        TA_data.append([df_name, MJD, TA_x])

  # append TA_x to "/TA-GEN/TAx_xxxx.dat"
    for nclk in range(len(TA_data)):
        name = TA_data[nclk][0]
        t = TA_data[nclk][1]
        x = TA_data[nclk][2]
        t_str = '{:5d}'.format(t)
        x_str = '{:+.6e}'.format(x)
        TA_str = str(t_str) + " " + str(x_str) + "  *"
        #infile = TAgen_dir + name + "_TAx.dat"                                    
        infile = TAgen_dir + "TAx_" + name +".dat"                                 
        with open(infile, 'a') as infile:
            infile.write(TA_str + '\n')

  # append TA_x to the TAx_list[]
    TAx_list_new = []
    TA_w = []                                                                      
    for nclk in range(len(TAx_list)):     # TAx_list consists of only clock data
        TA_name = TAx_list[nclk][0]
        TA_t = TAx_list[nclk][1][0]
        TA_x = TAx_list[nclk][1][1]
        flg = TAx_list[nclk][3]
        cid = TAx_list[nclk][4]                                                    
      # actually-used-weight                                                      
        for i in range(len(wgt_list)):                                            
            if (TA_name == wgt_list[i][0]):                                       
                wgt = wgt_list[i][1]                                              
      # append the [t,x] to the TAx_list
        for n in range(len(TA_data)):
            if (TA_data[n][0] == TA_name):
                t = TA_data[n][1]
                x = TA_data[n][2]
        TA_t = np.hstack((TA_t, t))
        TA_x = np.hstack((TA_x, x))
        TA_w = np.hstack((TA_w, wgt))                                              
      # TA_list_new[] is the renewd TAx_list[]
        data_TAx = np.array([TA_t, TA_x])  # data_raw = [[MJD,,,],[phase,,,]]
        TAx_list_new.append([TA_name, data_TAx, TA_w, flg, cid])                   

  # ................................................................               
  # formatting for output file
    TAxhat_out = []                          
    for n in range(len(TAxhat_list)):        
        name = TAxhat_list[n][0]             
        t1   = TAxhat_list[n][1]             
        x0   = TAxhat_list[n][2]              
        rate = TAxhat_list[n][3]             
        xhat = TAxhat_list[n][4]             
        flg  = TAxhat_list[n][6]
        cid  = TAxhat_list[n][7]                                                                
        for i in range(len(wgt_list)):       
            if (name == wgt_list[i][0]):     
                wgt = wgt_list[i][1]
      # pick up the weights
        len_wgtlist = len(wgt_addinfo)                              
        j = 0                                                       
        wgt_adev = np.nan
        wgt_adev_tau = np.nan
        while (j < len_wgtlist):
            if (name == wgt_addinfo[j][0]):
                wgt_adev_tau = wgt_addinfo[j][1]                                   
                wgt_adev = wgt_addinfo[j][2]
                if not np.isnan(wgt_adev_tau):wgt_adev_tau0 = wgt_adev_tau 
            j += 1            
      # formatting for output file
        t1 = '{:5d}'.format(int(t1))
        x0 = '{:+.6e}'.format(x0)
        rate = '{:+.6e}'.format(rate)
        xhat = '{:+.6e}'.format(xhat)
        adv = '{:+.6e}'.format(wgt_adev)                                          
        mwgt = '{:.6f}'.format(wgt)
        flg = '{:1d}'.format(flg)
        cid = '{:.0f}'.format(float(cid))                                          
        data = (str(t1) + "  " + name + "  " 
                + str(x0) + "  " + str(rate) + "  " + str(xhat) + "   " 
                + str(adv) + "  " + str(mwgt) + "   " 
                + str(flg) + "    " + str(cid))                                      
      # TAxhat_out[] is for the outfile
        TAxhat_out.append([data])
      
  # save TAxhat_out[] ==> "./TA_GEN/CSxx_TAxhat.dat"
    title = "t0  clock     x(t0-1day)   rate(" + str(TArate_span) + "days)"
    title = title + "     xhat(t0)      ADEV(" + str(wgt_adev_tau0) + "days)"                                       
    title = title + "    weight   flg1 cid"                                            
    outfile = TAgen_dir + str(MJD) + '_TAxhat.dat'
    np.savetxt(outfile, TAxhat_out, header=title, fmt='%s')

  # comment to log[]n
    comment = '**2** make a timescale with prediction *****' + '\n'
    log_0 = "::::: saved TAxhat data : "+ TAgen_dir + "(MJD)_TAxhat.dat" + '\n'               
    log_1 = "::::: saved TAx data : "+ TAgen_dir + "TAx_(clock).dat" + '\n'                   
    log_2 = "::::: interval for rate estimation (TArate_span)      : "+str(TArate_span)+'(days)\n'                     
    log_3 = "::::: ADEV tau for clock weight (wgt_adev_tau0)       : "+str(wgt_adev_tau0)+'(days)\n'                   
    log_4 = "::::: max weight limit (wgt_lmt)                      : "+str(wgt_lmt)+'\n'                               
    log_5 = "::::: sample numbers for ADEV for weighting (Twgt_N)) : "+str(Twgt_N)+'\n'                                
    log = [comment + log_0 + log_1 + log_2 + log_3 + log_4 + log_5]                   

    return TAx_list_new, wgt_list, comment, log                                    

# ==============================================================
# DEF : calculate clock weight from ADEV
def TA_calc_wgt(TAx_list, MJD, TAwgt_span, clk_initparam, x_limit, Twgt_N, 
                wgt_lmt, t_intvl):                                                    
    
    t0 = MJD - TAwgt_span*Twgt_N                                                    
    t1 = MJD - 1
    
  # read TAx data from TAx_list[] 
    clk_invAdev = []
    wgt_addinfo = []                                                              
    wsum = 0.0
    for nclk in range(len(TAx_list)):
        name = TAx_list[nclk][0]
        TA_t = TAx_list[nclk][1][0]
        TA_x = TAx_list[nclk][1][1]

      # set the initial_weights (wgt0) on the clocks
        for n in range(len(clk_initparam)):
            if (name == clk_initparam[n][0]): wgt0 = clk_initparam[n][1]
        
      # operation targets are the clocks with wgt0!=0 
        if (wgt0 != 0.0): 
       
          # pickup data within[MJD-1 : MJD-TAwgt_span*10]
            x = []
            for i in range(len(TA_t)):
                if ((int(TA_t[i])>=t0) & (int(TA_t[i])<=t1)):
                    x.append(TA_x[i]) 
                    
          # ADEV calculation for assigned period
          # This ADEV is for the sampling interval of 86400 !!!!!!!!!!
            taus = []
            adevs = []
            adev_list = []
            npx=np.isnan(x)
          
          # if ADEV calculation is possible                         
            if (any(xx >= x_limit for xx in x) != True 
                & (npx.any() != True)):                                                 
              # ADEV calculation
                (taus, adevs, errors, ns) = allantools.oadev(x,taus="octave")
                taus = taus*t_intvl
                adevs = adevs/t_intvl
                adev_list = np.array([taus, adevs])
              # get the closest ADEV at TAwgt_span  
                sa = [] 
                for i in range(len(adev_list[0])):
                    sa.append(abs(adev_list[0][i] - TAwgt_span*t_intvl))             
                adev_idx = sa.index(min(sa))
              # wgt_adev_tau : tau of ADEV(tau) for calculating weight     
                wgt_adev = adev_list[1][adev_idx]                                 
                wgt_adev_tau = adev_list[0][adev_idx]                             
                wgt_adev_tau = int(wgt_adev_tau/t_intvl)
              # Allantools gives zero, set adev to minimum limit
                if (wgt_adev == 0.0):                                
                    print("!!!!", MJD, name, "ADEV is set to be 1.0e-28")
                    wgt_adev = 1.0e-28                               
                print("    ", MJD, name, "adev(" + str(wgt_adev_tau) 
                      + "days)=", f"{wgt_adev:.6e}")                
              # calculate 1/AVAR(tau~TAwgt_span)
                inv_avar_1e10 = 1/(wgt_adev*1e10)**2                                              
              
          # if ADEV calculation is impossible                                         
            else:
                print("!!!!", MJD, name, "ADEV impossible")
                inv_avar_1e10 = 0.0                                                   
                wgt_adev = np.nan                                                 
                wgt_adev_tau = np.nan                                             
          
          # for normalization
            wgt_addinfo.append([name, wgt_adev_tau, wgt_adev])                    
            clk_invAdev.append([name,inv_avar_1e10])                                      
            wsum = wsum + inv_avar_1e10                                                
          
  # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: yuko250303
  # normalize the weights  
    if (wsum != 0):
        clk1 = []
        wgt1 = []
        for i in range(len(clk_invAdev)):
            clk1.append(clk_invAdev[i][0])
            wgt1.append(clk_invAdev[i][1]/wsum)                                       
        clk_wgt = [clk1, wgt1]

      # renew the weight under the max_limit                                           
        if (wgt_lmt != 'OFF'):
            wgt3 = re_weight_under_limit(wgt_lmt,wgt1)                                
            clk_wgt = [clk1, wgt3]                                                    

  # save data in TAx_list 
    TAx_wgt_list = []
    for nclk in range(len(TAx_list)):
        name = TAx_list[nclk][0]
        TA_t = TAx_list[nclk][1][0]
        TA_x = TAx_list[nclk][1][1]
        flg  = TAx_list[nclk][3]
        wgt  = 0.0
        cid  = TAx_list[nclk][4]
        for i in range(len(clk_wgt[0])):
            if (name == clk_wgt[0][i]):
              wgt = clk_wgt[1][i]

      # Only weights are renewed in TAx_list           
        data_TAx = np.array([TA_t, TA_x])  # data_raw = [[MJD,,,],[phase,,,]]
        TAx_wgt_list.append([name, data_TAx, wgt, flg, cid])                       
    
    return TAx_wgt_list, wgt_addinfo                                              

# ==============================================================
# DEF : calculate the difference of reference_clock(with flg=1) and other clocks
def TA_clk_clk(infile_all, MJDall, x_limit):                                       
  
  # data_list[] is the data from infile_all <== "all_raw_clk.dat"
    data_list = []
    data_list = read_jointdata(infile_all, MJDall, x_limit)                        
    
  # find the reference clock ("clk0") with "flg=1"
    for nclk in range(len(data_list)):
        flg = data_list[nclk][3]
        if (flg == 1):
            name0 = data_list[nclk][0]
            name0 = name0.split('.')[0]
            name0 = name0.split('-')[1]                                            
            t0 = data_list[nclk][1][0]
            x0 = data_list[nclk][1][1]
            t0 = np.array([t0[n] for n in range(len(t0)) 
                           if (t0[n] >= MJDall[0]) & (t0[n] <= MJDall[-1])])
            x0 = np.array([x0[n] for n in range(len(t0)) 
                           if (t0[n] >= MJDall[0]) & (t0[n] <= MJDall[-1])])

  # calculate the difference of "clk0" and other clocks
    diffdata = [t0]
    data_diff_list = []
    diffname = "MJD   "                                                            
    for nclk in range(len(data_list)):
        flg = data_list[nclk][3]
        if (flg != 1):
            name = data_list[nclk][0]
            name = name.split('.')[0]
            name = name.split('-')[1] + "-" + name0                                
            t = data_list[nclk][1][0]
            x = data_list[nclk][1][1]
            t = np.array([t[n] for n in range(len(t)) 
                          if (t[n] >= MJDall[0]) & (t[n] <= MJDall[-1])])
            x = np.array([x[n] for n in range(len(t)) 
                          if (t[n] >= MJDall[0]) & (t[n] <= MJDall[-1])]) 
          # outlier => np.nan
            diffx = np.array(x0 - x)    # because x-x0 = (UTCk-x0) - (UTCk-x)
            for i in range(len(diffx)):                                            
                if (abs(diffx[i]) >= x_limit):                                     
                    diffx[i] = np.nan                                              
            diffdata.append(diffx)
            
          # data_diff_list[] is [["clk0-clk", [MJD], [clk0-clk]],,]
            data_diff_list.append([name, t, diffx])
            diffname = diffname + name + "     "
        
  # save diffdata ==> "all_clk-Clk0.dat"
    ncol = len(diffdata)
    np.savetxt("all_clk-" + name0 + ".dat", np.transpose(diffdata), 
               fmt=' '.join(['%i'] + ['%+.6e']*(ncol-1)), header=diffname)

    return data_diff_list

# ==============================================================
# DEF : prepare the initial data for starting TA calculation
def TA_make_init_x(infile_raw, MJDall, TAgen_dir, TArate_span, Days_for_Adev,
                   x_limit):                                                       

  # TA_MJD0 is the day of starting the TA calculation                             
    Prepared_span = max(TArate_span, Days_for_Adev)                               
    TA_MJD0 = MJDall[0] + Prepared_span + 1                                       

  # save the output path and filename
  # os.makedirs(TAgen_dir)                ########## if directry exists, error

  # data_list[] is the data of "all_raw_clk.dat"
    data_list = []
    data_list = read_jointdata(infile_raw, MJDall, x_limit)                        

  # prepare & save initial TAx data
    TAx_init_list = []
    t = []
    x = []
    for nclk in range(len(data_list)):
        name = data_list[nclk][0]
        name = name.split('.')[0]
        name = name.split('-')[1]                                                  
        t = data_list[nclk][1][0]
        x = data_list[nclk][1][1]
        t = np.array([t[n] for n in range(len(t)) 
                      if (t[n] >= MJDall[0]) & (t[n] <= TA_MJD0 - 1)])
        x = np.array([x[n] for n in range(len(t)) 
                      if (t[n] >= MJDall[0]) & (t[n] <= TA_MJD0 - 1)])
        x = (-1)*x    # (UTCk - CS) ==> (CS - UTCk)
        wgt = data_list[nclk][2]
        flg = data_list[nclk][3]
        cid = data_list[nclk][4]                                                   
        
      # data_TAx[] is np.array([MJD],[CS-UTCk])
        data_TAx = np.array([t, x])  
        TAx_init_list.append([name, data_TAx, wgt, flg, cid])                      

      # save each [[MJD],[CS-UTCk]] ==> "./TA_GEN/CSxx_TAx.dat"
        outfile = TAgen_dir + 'TAx_' + name +'.dat'                                
        title = "t0      x(t0)"
        np.savetxt(outfile, np.transpose([t, x]), header=title, 
                   fmt=' '.join(['%i'] + ['%+.6e']))

    return TAx_init_list, TA_MJD0

# ==============================================================
# DEF :  calculate a TAx prediction (= xhat) for next day
def TA_predict_x(TAx_list, MJD, TAgen_dir, TArate_span):                 
    
    # read from TAx_list[] (until one day before of the calculating date)
    TAxhat_list = []
    for nclk in range(len(TAx_list)):
        name = TAx_list[nclk][0]
        name = name.split('.')[0]
        t = TAx_list[nclk][1][0]
        x = TAx_list[nclk][1][1]
        wgt = TAx_list[nclk][2]
        flg = TAx_list[nclk][3]
        cid = TAx_list[nclk][4]                                                    
        
      # calculate rate & predicted-TAx on each MJD
        for i in range(len(t)):
            if (t[i] == (MJD-1)):  # rate calculation goes back from yesterday 

              # pick-up the the latest TAx for rate calculation
                t_block = np.array([t[i] for i in range(len(t)) 
                                    if (t[i] >= MJD-(TArate_span+1)) 
                                    & (t[i] <= MJD-1)])
                x_block = np.array([x[i] for i in range(len(t)) 
                                    if (t[i] >= MJD-(TArate_span+1)) 
                                    & (t[i] <= MJD-1)])
              # calculate xhat as the predicted_TAx for next day
                slope = (x_block[-1] - x_block[0]) / (t_block[-1] - t_block[0])
                xhat = x[i] + slope * (t[i] - t[i-1])
                
              # Caution to the time-tag of saved data! 
              # TAxhat_list[] is [[name,t1,x(t0),rate,xhat(t1),wgt,flg,cid],,,]
                xhat_data = [name, t[i]+1, x[i], slope/86400, xhat, 
                             wgt, flg, cid]                                                               
                TAxhat_list.append(xhat_data)

    return TAxhat_list

# ==============================================================
# DEF :  get UTCk-clk 
def TA_ref_TA(infile_raw, TAx_list, comment, TA_MJD0, x_limit):                    

  # get the "clk0-TA" data from TAx_list[]
    for nclk in range(len(TAx_list)):
        flg = TAx_list[nclk][3]
        if (flg == 1):
            MJD = TAx_list[nclk][1][0]
            clk0_TA = TAx_list[nclk][1][1]
            MJD0 = int(np.amin(MJD))
            MJD1 = int(np.amax(MJD))

  # get the "UTCk-clk0" data from the "all_raw_clk.dat"
    data_list = []
    data_list = read_jointdata(infile_raw, MJDall, x_limit)                        

    for nclk in range(len(data_list)):
        flg = data_list[nclk][3]
        if (flg == 1):   # pick up the clk0
            name0 = data_list[nclk][0]
            name0 = name0.split('.')[0]
            name0 = name0.split('-')[1]                                            
            t = data_list[nclk][1][0]
            x = data_list[nclk][1][1]
            nt = [t[i] for i in range(len(t)) 
                  if (t[i] >= MJD0) & (t[i] <= MJD1)]        
            nx = [x[i] for i in range(len(t)) 
                  if (t[i] >= MJD0) & (t[i] <= MJD1)]
            UTCk_clk0 = nx

  # make the "UTCk-TA" from "clk0-TA"+"UTCk-clk0"
    UTCk_TA = []
    UTCk_TA_list = data_list
    for i in range(len(MJD)):
        if (MJD[i] == nt[i]):
            sum = clk0_TA[i] + UTCk_clk0[i]
            if (MJD[i] < TA_MJD0):                                                 
                sum = np.nan                                                       
            UTCk_TA.append(sum)
    data = np.array([MJD, UTCk_TA])
    UTCk_TA_list.append(["UTCk-TA", data, "0.0", "3", "0"])                        

    return UTCk_TA_list, comment                                                   

# ==============================================================
# DEF : define the span for rate estimation 
def TA_set_span_for_RateWgt(MJDall):
    
  # TArate_span = interval for estimating the clock rate   
    print('\n**2** Set the span for clock rate calculation. ')
    chk1 = False
    while chk1 == False:
        try:                                                                        
            TArate_span = int(input('??2?? Span in days (recommended 10 or 30 or 90) = '))
            if (TArate_span >= len(MJDall)):
                print('!!2!! Data is not enough, select smaller span.')
            else:
                print('**2** The rate span = ', TArate_span)
                chk1 = True
        except ValueError:                                                         
            print("\n!!2!! Input is wrong. Try again !!")                          

  # TAwgt_span = ADEV_tau for calcluating clock weight                            
    print('\n**2** Set the ADEV_tau for setting the clock weight. ')               
    chk1 = False                                                                  
    while chk1 == False:                                                          
        try:                                                                       
            TAwgt_span = int(input('??2?? Span in days (16 or 32 or 64 or 128) = '))              
            if (TAwgt_span >= len(MJDall)):                                       
                print('!!2!! Data is not enough, select smaller span.')                
            else:                                                                 
                print('**2** The ADEV_tau = ', TAwgt_span)                        
                chk1 = True                                                       
        except ValueError:                                                         
            print("\n!!2!! Input is wrong. Try again !!")                          
            
    return TArate_span, TAwgt_span                                                
    
# ==============================================================   
# DEF : set the weight limit or not 
def TA_set_weight_limit(wgt_lmt0):                                                    

    while (wgtlmt_yn := input ('\n??2?? Set the weight_limit? (y/n) = ')) not in ('y','n'):
        print("!!!! Input is not correct. Please try again!")
    if wgtlmt_yn == 'y' : wgt_lmt = wgt_lmt0
    if wgtlmt_yn == 'n' : wgt_lmt = 'OFF'
    print('**2** weight limit = ', wgt_lmt)  

    return wgt_lmt
    
# ==============================================================
# main program
# ==============================================================
# select the process menu

menu_n = 0
data_list = []
TAx_list = []
diff_list = []
clk_initparam = []

# read parameter file                                                             
t_intvl, x_limit, Twgt_N, wgt_lmt0 = read_param_file()                            

# read initial clock data
print("\n" + "..... Reading initial data of 'UTC(k) - Clock'")                 
data_list, log, comment = read_rawdata(t_intvl, x_limit)                                

# Search maximum MJD span (instead of common span)
MJDall = MJD_max(data_list)

# joint clock data and output ==> "all_raw_clk.dat"
print("..... formatting the initial data for processing")
outfile_raw = "all_raw_clk.dat"

# outfile_raw, log = joint_rawdata(outfile_raw, data_list, log)                    
outfile_raw, log = joint_rawdata(outfile_raw, MJDall, data_list, log)
outfile = outfile_raw

# Show the initial status of all clocks
print("\n************ Initial status of all clocks ************")
menu_n = 0
adev_plot(outfile, "", [], MJDall, log, comment, menu_n, x_limit)                      

# select data processing while int(menu_n) < 10:
chk0 = False                                                                        
while (chk0 == False):                                                              

    print('')
    print('=======================================')
    print('      Menu of the data processing      ')
    print('---------------------------------------')
    print('  0.  Plot the "UTC(k)-Clock" data'     )
    print('  1.  Assign the MJD range'            )
    print('  2.  Make TA with daily prediction'    )
    print('              and dynamical weighting'  )
    print(' 10.  Exit')
    print('=======================================')
    
    menu_n = input('????? Input menu number ? >> ')
    if (menu_n not in ['0','1','2','10']):                                          
        print ("\n!!!!! Input is wrong, Try again! ")                              

    if menu_n == '0':
        data_list = []
        data_list, log, comment = read_rawdata(t_intvl, x_limit)
        MJDall = MJD_max(data_list)
        outfile, log = joint_rawdata(outfile_raw, MJDall, data_list, log)
        adev_plot(outfile, "", [], MJDall, log, comment, menu_n, x_limit)                 

    if menu_n == '1':
        MJDall, log = MJD_limit(MJDall, log)
        outfile, log = joint_rawdata(outfile, MJDall, data_list, log)                 
        adev_plot(outfile, "", [], MJDall, log, comment, menu_n, x_limit)                 
        
    if menu_n == '2':
      # get the initial clock weight: clk_initparam=[clkname,wgt]
        infile_clklist = "clock_file_list.txt"                                     
        clk_initparam = clk_init_wgt(infile_clklist)
        
      # prepare the TAxij (= clk_i - clk_j) from the measurement raw data
        infile_raw = "all_raw_clk.dat"
        diff_list = TA_clk_clk(infile_raw, MJDall, x_limit)                        

      # prepare the directory for data saving                                      
        TAgen_dir = "./TA_GEN/"                                                     
        dt_str, dt_com, TAgen_path = set_dirpath(TAgen_dir)                        
        try:                                                                       
            os.makedirs(TAgen_path)                                                  
        except FileExistsError:                                                      
            print("Directry already exists.")                                      
            
      # prepare the parameters and initial TAx data for rate calculation
        TArate_span, TAwgt_span = TA_set_span_for_RateWgt(MJDall)                 
        
      # set the weight limit?                                       
        wgt_lmt = TA_set_weight_limit(wgt_lmt0)                                       

      # prepare the initial bunch of data for starting TA calculation
      # initdate_span is necessary range for Allan-Variance(tau=TAwgt_span)
        initdate_span = TAwgt_span * Twgt_N                                       
        TAx_list, TA_MJD0 = TA_make_init_x(infile_raw, MJDall, TAgen_path, TArate_span, 
                                           initdate_span, x_limit)                 
        
      # calculate the TAx for the assigned MJD   
        TAwgt_list = []
        all_wgt_list = []
        for MJD in range(TA_MJD0, MJDall[-1]+1):
            print('\n' + "***** Calculating a timescale for ", MJD)
          
          # calculate the clock weight from the Allan-variance of past TAx
            TAx_list, wgt_addinfo = TA_calc_wgt(TAx_list, MJD, TAwgt_span,  
                                                clk_initparam, x_limit, Twgt_N, 
                                                wgt_lmt, t_intvl)                     
            
          # calculate the TAx predictaion (= xhat) for the next day
          # Here, clock data interval is considered as one day (=86400sec)
            TAxhat_list = TA_predict_x(TAx_list, MJD, TAgen_path, 
                                       TArate_span)
            
          # calculate and save the TAx
            TAx_list, wgt_list, comment, log_TA = TA_calc_x(TAx_list, TAxhat_list, 
                                   diff_list, MJD, TAgen_path, TArate_span, TAwgt_span, 
                                   wgt_addinfo, x_limit, Twgt_N, wgt_lmt)                   

          # merge the weights of all clocks in all MJD
            all_wgt_list.append([MJD,wgt_list])                                    

      # save the clock weights data
        wgtfile = save_weight(all_wgt_list)                                        
          
      # get the (TA - UTCk) from (x + (UTCk-CS))
        UTCk_TA_list, comment = TA_ref_TA(infile_raw,TAx_list,comment,
                                               TA_MJD0, x_limit)                   
        
      # joint the clock data & TA ==> outfile (= "all_raw_clk_TA.dat")
        log = log + log_TA                                                           
        outfile, log = joint_rawdata("all_raw_clk_TA.dat", MJDall, 
                                     UTCk_TA_list, log)                            
       
      # plot & save the final results 
        adev_plot(outfile,wgtfile,all_wgt_list,MJDall,log, comment,menu_n,x_limit)                
        
    if menu_n == '10':
        chk0 == 0                                                                  
        break
