import os
import sys
import glob
import openpyxl
import datetime
import numpy as np
from numpy import sin, cos, tan, cosh, pi, exp, log, sqrt
import pandas as pd
from matplotlib import pyplot as plt


from tklib.tkutils import mprint, IsDir, IsFile, SplitFilePath, modify_path, delete_file
from tklib.tkutils import terminate, safe_getelement, pint, pfloat, getarg, getintarg, getfloatarg, safe_getelement
from tklib.tkutils import save_csv
from tklib.tkfile import tkFile
from tklib.tkapplication import tkApplication
from tklib.tkexcel import tkExcel
import tklib.tkre as tkre
from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me, Ea_Arrhenius


#=======================================
# Read Resitest Excel files and show van der Pauw graphs
#=======================================

usage_str = '''
""
"(i)   usage: python {} f".format(sys.argv[0])
"        Plot van der Pauw correction factor f as a function of resistance ratio"
"   ex: python {} {}".format(sys.argv[0], 'f')
"(ii)   usage: python {} plot filemask".format(sys.argv[0])
"        Plot measured Hall-related values as a function of time / temperature"
"        Save summary CSV files"
"   ex: python {} {} {}".format(f'{sys.argv[0]} plot {cparams.filemask}')
'''[1:-1]


def initialize():
#================================
# Global variables
#================================
    app          = tkApplication(usage_str  = usage_str)
    argv, narg   = app.get_argv()

    cparams             = app.get_params()
    cparams.debug       = 0
    cparams.print_level = 0

# mode: 'f', 'plot'    
    cparams.mode = 'plot'

    cparams.figsize         = (8, 4)
    cparams.figsize_plot    = (12, 8)
    cparams.fontsize        = 16
    cparams.legend_fontsize = 8

    cparams.filemask = '*.xlsx'

    return app, cparams

#=============================
# Treat argments
#=============================
def usage(app):
    cparams = app.get_params()
#    app.usage(infile = cparams.infile)
    if app.usage_str is None:
        return

    for s in app.usage_str.split('\n'):
        cmd = 'print({})'.format(s.rstrip())
        eval(cmd)

def updatevars(app, cparams):
    argv = app.argv
#    if len(argv) <= 1:
#        app.terminate(usage = usage)

    cparams.mode = getarg(1, cparams.mode)
    if cparams.mode == 'f':
        pass
    elif cparams.mode == 'plot':
        cparams.filemask = getarg(2, cparams.filemask)
    else:
        app.terminate(f"Error in updatevars(): Invalide mode [{cparams.mode}]", usage = usage)

def cal_f(r, f0 = 1.0):
    h       = 1.0e-7
    dump    = 0.0
    eps     = 0.5e-6
    maxiter = 100

    f = f0
    for i in range(maxiter):
        d0 = delta(f, r)
        if abs(d0) < eps:
            return f

        dh = delta(f + h, r)
        d1 = (dh - d0) / h
        if d1 < 0.0:
            d1d = d1 - dump
        else:
            d1d = d1 + dump
        df = -d0 / d1d
#        print("r=", r, "f=", f, "  d(delta)/df=", d1, "  df=", df, "  d=", d0)
        f += df

    return None
    
def rho_vanderPaw(d, dV, dI):
    k = pi * d / log(2.0) # d in cm
    return k * dV / dI

def read_value_range2(irow, key, inf, sheet):
    inf[f"{key}_1"] = pfloat(sheet[irow][2].value)
    val = sheet[irow][3].value.split()
    inf[f"{key}_1_range"] = [pfloat(val[0]), pfloat(val[2])]

    inf[f"{key}_2"] = pfloat(sheet[irow][4].value)
    val = sheet[irow][5].value.split()
    inf[f"{key}_2_range"] = [pfloat(val[0]), pfloat(val[2])]

def normalize_str(str):
    str = str.replace('～', '-')
    str = str.replace('年', '/')
    str = str.replace('月', '/')
    str = str.replace('日', '')

    return str
    
def get_sec_passed(daystr, timestr):
#    date_obj = datetime.datetime.strptime(f"{daystr} {timestr}", '%Y年%m月%d日 %H:%M:%S')
    date_obj = datetime.datetime.strptime(f"{daystr} {timestr}", '%Y/%m/%d %H:%M:%S')
    delta = date_obj - datetime.datetime(2022, 1, 1, 0, 0, 0)
    dsec = delta.days * 24 * 60 * 60 + delta.seconds
    return dsec

def read_resitest_excel(f):
    print(f"Read [{f}]")

    inf = {}

    fp = tkExcel(f)
    if fp is None:
        print("")
        print(f"Error in read_resitest_excel(): Can not read [{f}]")
        return None

    inf['file_path'] = f

    s_config = fp.Read_sheet(sheet_name = 0)
    try:
        inf["measurement_type"] = s_config[9][1].value
        inf["comment"]          = s_config[16][1].value
        inf["source_path"]      = s_config[2][0].value
        inf["thickness"]        = pfloat(s_config[10][1].value)
    except:
        fp.Close()
        return None

    val = s_config[10][0].value
    if 'nm' in val:
        inf["thickness"] *= 1.0e-7    # nm => cm
        inf["thickness_unit"] = 'nm'
    elif 'mm' in val:
        inf["thickness"] *= 0.01      # mm => cm
        inf["thickness_unit"] = 'mm'
    elif 'um' in val:
        inf["thickness"] *= 1.0e-4    # um => cm
        inf["thickness_unit"] = 'um'
    else:
        print("")
        print("Error in read_resitest_excel(): Invalid unit [{}]".format(inf["thickness_unit"]))
        print("")
        exit()

    s_summary = fp.Read_sheet(sheet_name = 2)

    inf["carrier_type"] = s_summary[23][2].value
    val = s_summary[23][3].value.split('%')
    inf["carrier_type_accuracy"] = pfloat(val[0]) * 0.01  # % => float

    def read_value_range(irow, key, inf, sheet):
        inf[key] = pfloat(sheet[irow][2].value)
        val = sheet[irow][3].value.split()
        inf[f"{key}_range"] = [pfloat(val[0]), pfloat(val[2])]
    
    read_value_range(22, 'mu_Hall',  inf, s_summary)
    read_value_range(24, 'N',        inf, s_summary)
    read_value_range(25, 'Nsheet',   inf, s_summary)
    read_value_range(26, 'RH',       inf, s_summary)
    read_value_range(27, 'RHsheet',  inf, s_summary)
    read_value_range(28, 'rho',      inf, s_summary)
    read_value_range(29, 'rhosheet', inf, s_summary)
    read_value_range(29, 'VH',       inf, s_summary)
    inf['VH_phase'] = pfloat(s_summary[31][2].value)

    s_rho = fp.Read_sheet(sheet_name = 3)
    inf['day_start']  = normalize_str(s_rho[14][2].value)
    inf['day_end']    = normalize_str(s_rho[14][3].value)
    inf['time_start'] = s_rho[15][2].value
    inf['time_end']   = s_rho[15][3].value
    inf['T_start']    = pfloat(s_rho[16][2].value)
    inf['T_end']      = pfloat(s_rho[16][3].value)
    inf['T_average']  = pfloat(s_rho[16][4].value)

    inf["dsec_start"] = get_sec_passed(inf['day_start'], inf['time_start'])
    inf["dsec_end"]   = get_sec_passed(inf['day_end'],   inf['time_end'])

    read_value_range2(24, 'rho_final',      inf, s_rho)
    read_value_range2(25, 'rhosheet_final', inf, s_rho)
    inf["F_1"] = pfloat(s_rho[26][2].value)
    inf["F_2"] = pfloat(s_rho[26][4].value)

    inf["rho2134"] = pfloat(s_rho[29][2].value)
    inf["rho3241"] = pfloat(s_rho[29][3].value)
    inf["rho4312"] = pfloat(s_rho[29][4].value)
    inf["rho1423"] = pfloat(s_rho[29][5].value)
    inf["rho1423"] = pfloat(s_rho[29][5].value)

    def read_vanderPauw(icolumn, key, inf, sheet):
        i0 = 38
        ic = icolumn
        inf[f"{key}_+_averageV"] = pfloat(sheet[i0  ][ic].value)
        inf[f"{key}_+_errorV"]   = pfloat(sheet[i0+1][ic].value)
        inf[f"{key}_+_averageI"] = pfloat(sheet[i0+2][ic].value)
        inf[f"{key}_+_errorI"]   = pfloat(sheet[i0+3][ic].value)
        inf[f"{key}_+_currentV"] = pfloat(sheet[i0+4][ic].value)

        i0 = 44
        inf[f"{key}_-_averageV"] = pfloat(sheet[i0  ][ic].value)
        inf[f"{key}_-_errorV"]   = pfloat(sheet[i0+1][ic].value)
        inf[f"{key}_-_averageI"] = pfloat(sheet[i0+2][ic].value)
        inf[f"{key}_-_errorI"]   = pfloat(sheet[i0+3][ic].value)
        inf[f"{key}_-_currentV"] = pfloat(sheet[i0+4][ic].value)

    read_vanderPauw(2, 'R2134', inf, s_rho)
    read_vanderPauw(3, 'R3241', inf, s_rho)
    read_vanderPauw(4, 'R4312', inf, s_rho)
    read_vanderPauw(5, 'R1423', inf, s_rho)

    s_Hall = fp.Read_sheet(sheet_name = 4)
    inf['day_start_Hall']  = normalize_str(s_Hall[19][2].value)
    inf['day_end_Hall']    = normalize_str(s_Hall[19][3].value)
    inf['time_start_Hall'] = s_Hall[20][2].value
    inf['time_end_Hall']   = s_Hall[20][3].value
    inf['T_start_Hall']    = pfloat(s_Hall[21][2].value)
    inf['T_end_Hall']      = pfloat(s_Hall[21][3].value)
    inf['T_average_Hall']  = pfloat(s_Hall[21][4].value)
    inf['B_start_Hall']    = pfloat(s_Hall[22][2].value)
    inf['B_end_Hall']      = pfloat(s_Hall[22][3].value)

    inf["dsec_start_Hall"] = get_sec_passed(inf['day_start_Hall'], inf['time_start_Hall'])
    inf["dsec_end_Hall"]   = get_sec_passed(inf['day_end_Hall'],   inf['time_end_Hall'])

    inf["C_averageVH"] = pfloat(s_Hall[29][2].value)
    val = s_Hall[29][3].value.split()
    inf["C_VH_range"]  = [pfloat(val[0]), pfloat(val[2])]
    inf["D_averageVH"] = pfloat(s_Hall[29][5].value)
    val = s_Hall[29][6].value.split()
    inf["D_VH_range"]  = [pfloat(val[0]), pfloat(val[2])]

    
    """
    R2134 = rho_vanderPaw(inf['thickness'], 
                    inf['R2134_+_averageV'] - inf['R2134_-_averageV'],
                    inf['R2134_+_averageI'] - inf['R2134_-_averageI'])
    print("R2134=", R2134)
    """

    fp.Close()
    
    return inf


def plot_Hall(app, cparams):
    print("")
    print("################################################")
    print("# Plot Hall data measured by ResiTest 8400")
    print("################################################")
    print(f"Search [{cparams.filemask}]")

    files = glob.glob(cparams.filemask)
#    print("files=", files)
    list = []
    for f in files:
        if tkre.Match(r'~', f):
            continue

        print(f"** Read [{f}]")
        inf = read_resitest_excel(f)
        if inf is None:
            print(f"Warning from read_resitest_excel(): File [{f}] is not a ResiTest 8400 file.")
#            print(f"Warning: [{f}] does not exist. Skip.")
            continue

        list.append(inf)

#    print("list=", list)
    list_sorted = sorted(list, key = lambda x: x['dsec_start'])

    t0 = list_sorted[0]['dsec_start']
    xt       = []
    yT       = []
    yF       = []
    yRR      = []

    xtav     = []
    xTav     = []
    yrhoav   = []
    yrho_err = []
    ysigmaav = []
    ymuav    = []
    ymu_err  = []
    yNav     = []
    yN_err   = []
    yRHav    = []
    yR2134   = []
    yV_Ip2134 = []
    yV_Ip3241 = []
    yV_Ip4312 = []
    yV_Ip1423 = []
    yV_Im2134 = []
    yV_Im3241 = []
    yV_Im4312 = []
    yV_Im1423 = []
    yV_IpCur2134 = []
    yV_IpCur3241 = []
    yV_IpCur4312 = []
    yV_IpCur1423 = []
    yV_ImCur2134 = []
    yV_ImCur3241 = []
    yV_ImCur4312 = []
    yV_ImCur1423 = []
    
    xtHall  = []
    yVH     = []
    yVH_err = []
    for inf in list_sorted:
        f = inf['file_path']
        print(f"{f}: {inf['day_start']} {inf['time_start']}")
        xt.append((inf['dsec_start'] - t0) / 60.0 / 60.0)
        yT.append(inf['T_start'])
        yF.append(inf['F_1'])
        yRR.append(inf['rho2134'] / inf['rho3241'])
        yRR.append(inf['rho4312'] / inf['rho1423'])

        xt.append((inf['dsec_end'] - t0) / 60.0 / 60.0)
        yT.append(inf['T_end'])
        yF.append(inf['F_2'])

        xtav.append(((inf['dsec_start'] + inf['dsec_end']) / 2.0 - t0) / 60.0 / 60.0)
        xTav.append((inf['T_start'] + inf['T_end'] + inf['T_start_Hall'] + inf['T_end_Hall']) / 4.0)

        yrhoav.append(inf['rho'])
        ysigmaav.append(1.0 / inf['rho'])
        maxerr = max([abs(inf['rho'] - inf['rho_range'][0]), abs(inf['rho'] - inf['rho_range'][1])])
        yrho_err.append(maxerr)

        ymuav.append(inf['mu_Hall'])
        ymu_err.append(max([abs(inf['mu_Hall'] - inf['mu_Hall_range'][0]), abs(inf['mu_Hall'] - inf['mu_Hall_range'][1])]))

        yNav.append(inf['N'])
        if inf['carrier_type'] == 'N':
            yRHav.append(-1.0 / e / inf['N'])
        else:\
            yRHav.append(1.0 / e / inf['N'])
        yN_err.append(max([abs(inf['N'] - inf['N_range'][0]), abs(inf['N'] - inf['N_range'][1])]))

        yV_Ip2134.append(inf['R2134_+_averageV'] / inf['R2134_+_averageI'])
        yV_Ip3241.append(inf['R3241_+_averageV'] / inf['R3241_+_averageI'])
        yV_Ip4312.append(inf['R4312_+_averageV'] / inf['R4312_+_averageI'])
        yV_Ip1423.append(inf['R1423_+_averageV'] / inf['R1423_+_averageI'])
        yV_Im2134.append(inf['R2134_-_averageV'] / inf['R2134_-_averageI'])
        yV_Im3241.append(inf['R3241_-_averageV'] / inf['R3241_-_averageI'])
        yV_Im4312.append(inf['R4312_-_averageV'] / inf['R4312_-_averageI'])
        yV_Im1423.append(inf['R1423_-_averageV'] / inf['R1423_-_averageI'])

        yV_IpCur2134.append(inf['R2134_+_currentV'] / inf['R2134_+_averageI'])
        yV_IpCur3241.append(inf['R3241_+_currentV'] / inf['R3241_+_averageI'])
        yV_IpCur4312.append(inf['R4312_+_currentV'] / inf['R4312_+_averageI'])
        yV_IpCur1423.append(inf['R1423_+_currentV'] / inf['R1423_+_averageI'])
        yV_ImCur2134.append(inf['R2134_-_currentV'] / inf['R2134_-_averageI'])
        yV_ImCur3241.append(inf['R3241_-_currentV'] / inf['R3241_-_averageI'])
        yV_ImCur4312.append(inf['R4312_-_currentV'] / inf['R4312_-_averageI'])
        yV_ImCur1423.append(inf['R1423_-_currentV'] / inf['R1423_-_averageI'])

        xtHall.append((inf['dsec_start_Hall'] - t0) / 60.0 / 60.0)
        yVH.append(inf['C_averageVH'])
        yVH_err.append(max([abs(inf['C_averageVH'] - inf['C_VH_range'][0]), abs(inf['C_averageVH'] - inf['C_VH_range'][1])]))

        xtHall.append((inf['dsec_end_Hall'] - t0) / 60.0 / 60.0)
        yVH.append(inf['D_averageVH'])
        yVH_err.append(max([abs(inf['D_averageVH'] - inf['D_VH_range'][0]), abs(inf['D_averageVH'] - inf['D_VH_range'][1])]))

    outcsvfile = 'summary.csv'
    keys = [key for key in list_sorted[0].keys()]
    print("")
    print(f"Save summry CSV to [{outcsvfile}]")
    fp = tkFile(outcsvfile, 'w')
    for key in keys:
        fp.Write(f'"{key}",')
    fp.Write('\n')
    for inf in list_sorted:
        for key in keys:
            t = type(inf[key])
#            print("key:", key, t, t is list, t is float, t is int)
            if inf[key] is None:
                fp.Write(',')
            elif type(inf[key]) is float or type(inf[key]) is int or type(inf[key]) is str:
                fp.Write(f"{inf[key]},")
            else:
                fp.Write(f'"[{inf[key][0]} - {inf[key][1]}]",')
        fp.Write('\n')
    fp.Close()

    HallTcsvfile = 'Hall-T.csv'
    print("")
    print("Save Hall - T data to [{}]".format(HallTcsvfile))
    save_csv(HallTcsvfile, ["T(K)", "sigma(S/cm)", "mu(cm2/Vs)", "N(cm-3)", "RH(cm3/C)"],
            [xTav, ysigmaav, ymuav, yNav, yRHav])

    HallTxlsfile = 'Hall-T.xlsx'
    df = pd.DataFrame(np.array([xTav, ysigmaav, ymuav, yNav, yRHav]).T, columns = ["T(K)", "sigma(S/cm)", "mu(cm2/Vs)", "N(cm-3)", "RH(cm3/C)"])
    print("")
    print("Save Hall - T data to [{}]".format(HallTxlsfile))
    df.to_excel(HallTxlsfile, index = False, header = True)
    
#=============================
# グラフの表示
#=============================
    app.fig = plt.figure(figsize = cparams.figsize_plot)
    fig             = app.fig
    fontsize        = cparams.fontsize
    legend_fontsize = cparams.legend_fontsize
    
    axTt   = fig.add_subplot(2, 3, 1)
    axFt   = fig.add_subplot(2, 3, 2)
    axRRt  = axFt.twinx()
    axVHt  = fig.add_subplot(2, 3, 3)
    axRt   = fig.add_subplot(2, 3, 4)
    axRCt  = axRt.twinx()
    axrhot = fig.add_subplot(2, 3, 5)
    axNt   = fig.add_subplot(2, 3, 6)
    axmut  = axNt.twinx()

    axTt.set_title(list[0]['file_path'], fontsize = fontsize)

    ins1 = axTt.plot(xt, yT, label = 'T (K)', color = 'black',  linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    axTt.set_xlabel("t (h)", fontsize = fontsize)
    axTt.set_ylabel("T (K)", fontsize = fontsize)
    ins = ins1 #+ ins2
    axTt.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
    axTt.tick_params(labelsize = fontsize)

    ins1 = axFt.plot (xt, yF,  label = '$f$',                                                color = 'black', linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins2 = axRRt.plot(xt, yRR, label = '$R_{2134}$ $/ R_{3241}$, $R_{4312}$ $/ R_{1423}$',   color = 'red',   linestyle = '-', linewidth = 0.5, marker = 'x', markersize = 5.0)
    axFt.set_xlabel  ("t (h)", fontsize = fontsize)
    axFt.set_ylabel  ("$f$ ", fontsize = fontsize)
    axRRt.set_ylabel ("$R_{2134}$ $/ R_{3241}$, $R_{4312}$ $/ R_{1423}$ ", fontsize = fontsize)
    axRRt.yaxis.label.set_color('red')   
    ins = ins1 + ins2
    axFt.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
    axFt.tick_params (labelsize = fontsize)
    axRRt.tick_params(labelsize = fontsize)
    axRRt.tick_params(axis = 'y', colors = 'red', labelsize = fontsize)

    axVHt.plot    (xtHall, yVH, label = r'$V_H$', color = 'black',  linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0) # capsize = 2
    axVHt.errorbar(xtHall, yVH, yerr = yVH_err, capsize = 2, color = 'black',  linestyle = '',  linewidth = 0.5) # capsize = 2
    axVHt.set_xlabel("$t$ (h)", fontsize = fontsize)
    axVHt.set_ylabel(r"$V_H$ (V)", fontsize = fontsize)
    axVHt.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
    axVHt.tick_params(labelsize = fontsize)

    ins1 = axRt.plot (xtav, yV_Ip2134,  label = '$V_{2134}/I_{2134}+$', color = 'green', linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins2 = axRt.plot (xtav, yV_Ip3241,  label = '$V_{3241}/I_{3241}+$', color = 'green', linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins3 = axRt.plot (xtav, yV_Ip4312,  label = '$V_{4312}/I_{4312}+$', color = 'green', linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins4 = axRt.plot (xtav, yV_Ip1423,  label = '$V_{1423}/I_{1423}+$', color = 'green', linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins5 = axRt.plot (xtav, yV_Im2134,  label = '$V_{2134}/I_{2134}-$', color = 'blue',  linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins6 = axRt.plot (xtav, yV_Im3241,  label = '$V_{3241}/I_{3241}-$', color = 'blue',  linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins7 = axRt.plot (xtav, yV_Im4312,  label = '$V_{4312}/I_{4312}-$', color = 'blue',  linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins8 = axRt.plot (xtav, yV_Im1423,  label = '$V_{1423}/I_{1423}-$', color = 'blue',  linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0)
    ins9  = axRCt.plot(xtav, yV_IpCur2134, label = '$V_{I,2134}/I_{2134}+$', color = 'purple',   linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    ins10 = axRCt.plot(xtav, yV_IpCur3241, label = '$V_{I,3241}/I_{3241}+$', color = 'purple',   linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    ins11 = axRCt.plot(xtav, yV_IpCur4312, label = '$V_{I,4312}/I_{4312}+$', color = 'purple',   linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    ins12 = axRCt.plot(xtav, yV_IpCur1423, label = '$V_{I,1423}/I_{1423}+$', color = 'purple',   linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    ins13 = axRCt.plot(xtav, yV_ImCur2134, label = '$V_{I,2134}/I_{2134}-$', color = 'red', linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    ins14 = axRCt.plot(xtav, yV_ImCur3241, label = '$V_{I,3241}/I_{3241}-$', color = 'red', linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    ins15 = axRCt.plot(xtav, yV_ImCur4312, label = '$V_{I,4312}/I_{4312}-$', color = 'red', linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    ins16 = axRCt.plot(xtav, yV_ImCur1423, label = '$V_{I,1423}/I_{1423}-$', color = 'red', linestyle = '-', linewidth = 0.5, marker = '^', markersize = 5.0)
    axRt.set_xlabel ("t (h)", fontsize = fontsize)
    axRt.set_ylabel (r"$V/I$ ($\Omega$)", fontsize = fontsize)
    axRCt.set_ylabel(r"$V_I/I$ ($\Omega$)", fontsize = fontsize)
    axRt.yaxis.label.set_color('blue')   
    axRCt.yaxis.label.set_color('red')   
    ins = ins1 + ins2 + ins3 + ins4 + ins5 + ins9 + ins13
    axRt.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
    axRt.tick_params(labelsize = fontsize)
    axRt.tick_params(axis = 'y', colors = 'blue', labelsize = fontsize)
    axRCt.tick_params(axis = 'y', colors = 'red', labelsize = fontsize)

    ins1 = axrhot.plot(xTav, yrhoav, label = r'$\rho$', color = 'black',  linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 5.0) # capsize = 2
    axrhot.errorbar   (xTav, yrhoav, yerr = yrho_err, capsize = 2, color = 'black',  linestyle = '',  linewidth = 0.5) # capsize = 2
    axrhot.set_xlabel("$T$ (K)", fontsize = fontsize)
    axrhot.set_ylabel(r"$\rho$ ($\Omega$cm)", fontsize = fontsize)
#    ins = ins1 + in
#    axrhot.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
    axrhot.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
    axrhot.tick_params(labelsize = fontsize)

    ins1 = axNt.plot (xTav, yNav,  label = '$N$',   color = 'black', linestyle = '-', linewidth = 0.5)#, marker = 'o', markersize = 2.0)
    axNt.errorbar    (xTav, yNav,  yerr = yN_err, capsize = 2, color = 'black', linestyle = '',  linewidth = 0.5)
    ins2 = axmut.plot(xTav, ymuav, label = '$\mu$', color = 'red',   linestyle = '-', linewidth = 0.5)#, marker = '^', markersize = 2.0)
    axmut.errorbar   (xTav, ymuav, yerr = ymu_err, capsize = 2, color = 'red',   linestyle = '',  linewidth = 0.5)
    axNt.set_xlabel ("$T$ (K)", fontsize = fontsize)
    axNt.set_ylabel ("$N_{Hall}$ (cm$^{-3}$)",   fontsize = fontsize)
    axNt.set_yscale('log')
    axmut.set_ylabel("$\mu_{Hall}$ (cm$^2$/Vs)", fontsize = fontsize)
    axmut.yaxis.label.set_color('red')   
    ins = ins1 + ins2
    axNt.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
    axNt.tick_params(labelsize = fontsize)
    axNt.tick_params(axis = 'y', labelsize = fontsize)
    axmut.tick_params(labelsize = fontsize)
    axmut.tick_params(axis = 'y', colors = 'red', labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.1)

    print("Press ENTER to exit>>", end = '')
    input()


def delta(f, r):
    x = log(2.0) / f
    kr = (r - 1.0) / (r + 1.0)
    return exp(x) / 2.0 - cosh(x * kr)

def plot_f(app, cparams):
    print("")
    print("Calculate f(R = RAB,CD/RBC,DA)")

    outputinterval = 50    
    outcsvfile = 'f-R.csv'

    r1 = 1000.0
    r0 = 1.0 / r1
    nr = 1001
    lnrstep = (log(r1) - log(r0)) / (nr - 1)
    
    xr = []
    yf = []
    f  = 0.3
    print("{:>10}\t{:>10}".format('R', 'f'))
    for i in range(nr):
        lnr = log(r0) + i * lnrstep
        r   = exp(lnr)
        f = cal_f(r, f)
        if f is None:
            app.terminate(f"Error in cal_f(): Could not converge for R={r}", usage = usage)                    

        if i % outputinterval == 0:
            print("{:10.4g}\t{:10.6g}".format(r, f))
        xr.append(r)
        yf.append(f)

    print("")
    print("Save fitting data to [{}]".format(outcsvfile))
    yfstr = ["{:10.6f}".format(yf[i]) for i in range(len(yf))]
    save_csv(outcsvfile, ["R", "f"], [xr, yfstr])

#=============================
# Plot graphs
#=============================
    fig = plt.figure(figsize = cparams.figsize)

    ax = fig.add_subplot(1, 1, 1)
    ax.plot(xr, yf, label = '$f$', color = 'black', linewidth = 0.5)
    ax.set_xscale('log')
    ax.set_xlabel("$R = R_{AB,CD} / R_{BC,DA}$",fontsize = cparams.fontsize)
    ax.set_ylabel("$f$",fontsize = cparams.fontsize)
    ax.set_ylim([0.0, 1.1])
#    ax.legend(fontsize = cparams.legend_fontsize, loc = 'best')
    ax.tick_params(labelsize = cparams.fontsize)

# Rearange the graph axes so that they are not overlapped
    plt.tight_layout()
    plt.pause(0.1)

    print("Press ENTER to exit>>", end = '')
    input()

    
def main(app, cparams):
    updatevars(app, cparams)

    if cparams.mode == 'f':
        plot_f(app, cparams)
    elif cparams.mode == 'plot':
        plot_Hall(app, cparams)
    else:
        app.terminate(f"Error in main(): Invalide mode [{cparams.mode}]", usage = usage)

    app.terminate("", usage = usage)


if __name__ == "__main__":
    app, cparams = initialize()
    main(app, cparams)