import os
import sys
import openpyxl
from math import exp, sqrt, log, gamma
import numpy as np
#from numpy import arange
from scipy import optimize          # newton関数はscipy.optimizeモジュールに入っている
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
from scipy import optimize          # newton関数はscipy.optimizeモジュールに入っている
from scipy.optimize import minimize


from tklib.tkutils import terminate, getarg, getintarg, getfloatarg, safe_getelement, validate_error
from tklib.tkutils import pint, pfloat, sort_lists
from tklib.tkinifile import tkIniFile
from tklib.tkvariousdata import tkVariousData
from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2
from tklib.tkapplication import tkApplication

from tklib.tktransport.tkTransport import read_datafile, FermiIntegral_x2
from tklib.tktransport.tkTransport import FermiIntegral_fast as FermiIntegral
from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, meff2NC_FEA, meff2DC0_FEA, BMShift_FEA
from tklib.tktransport.tkWeightedMobility import weighted_mobility, weighted_mobility_new, weighted_mobility_exact
from tklib.tktransport.tkTransport import LorentzNumber_FEA
from tklib.tktransport.tkDOS_FEA import integrate_Simpson, integrate_Simpson_list, tkDOS
from tklib.tktransport.tkmobility_tau import tkMobility, split_optstr


"""
Estimate m* from S and T
# G.J. Snyder, A. Pereyra and R. Gurunathan, Adv. Funct. Mater. 32, 2112772 (2022)
# Effective mass from Seebeck coefficient
# https://onlinelibrary.wiley.com/doi/10.1002/adfm.202112772
"""


#================================
# Global variables
#================================
Debug = 0

L_FEA = LorentzNumber_FEA

#mode: 'basic', 'prop', 'Hall', 'init', 'sim', 'fit', 'T', 'calL', 'calF'
mode = 'sim'

infile      = 'SnSeTe-S-Hall.xlsx'
T_label     = r'^T[\s|\(|\[|$]'
sigma_label = r'^sigma[$|\S]?.*$'
n_label     = r'^N[$|\S]?.*$'
mu_label    = r'^mu[$|\s|\(\[]?.*$'
S_label     = r'^S[$|\s|\(\[]?.*$'

outxlsFjfile   = "Fj.xlsx"
outxlsfFDfile  = "fFD.xlsx"
outxlsPropfile = "properties.xlsx"
outxlsSLfile   = "S-L.xlsx"
outcsvFitfile  = ""
outcsvLfile    = ""

parameterfile   = None
parameterbkfile = None
outxlsfile      = None

# Calculation ranges for mode = 'basic'
# x = E / kB / T
xmin_basic  = -20.0
xmax_basic  =  20.0
nx_basic    = 401
xstep_basic = (xmax_basic - xmin_basic) / (nx_basic - 1)

# Calculation ranges for mode = 'prop'
# x = E / kB / T
xmin  = -30.0
xmax  = 250.0
nx    = 281
xstep = (xmax - xmin) / (nx - 1)

# Seebeck coefficient
Smin  = 0.0
Smax  = 1.0e-3  # V/K
nS    = 101
Sstep = (Smax - Smin) / (nS - 1)

# N range
Nmin = 1.0e15  # cm-3
Nmax = 1.0e22  # cm-3
nN   = 101

# T range
T0 = 300.0   # K
Tmin = 300.0
Tmax = 800.0
nT   = 6

# read from DOSCAR
dos = tkDOS()
dos.meeff = 0.3
#dos.mheff = 1.0
#dos.EV = 0.0
dos.EC = 0.0
#dos.EA = 0.05   # Acceptor level, eV
#dos.NA = 0.0e17
#dos.ED = 1.05   # Donor level
#dos.ND = 0.0e17
#dos.EF0 = 0.0   # initial EF to find EF


#移動度パラメータ
mobility = tkMobility()
# For debug purpose. 1 will use 3-terms polynomial for the inverse of mobility
mobility.debug      = 0
mobility.use_simple = 0
#mobility.use_simple = 1
# charge, effective mass, scattering factor
mobility.charge = 1.0            # in e, 1.0 for hole, -1.0 for electron
mobility.meff   = dos.meeff      # in me
mobility.rfac   = 0.5            # tau = (meff/2)^0.5 * l0(T) * E^(r-0.5)
mobility.l0     = 1.0e-8         # m

# Lattice thermal conductivity
klatt  = 5.0            # W/m/K

# EF range for mode = 'EF'
#dEFmin    = -1.0  # measured from EV, eV
#dEFmax    =  1.0  # measured from EC, eV
#dos.nEF   = 50
#dos.Estep = 0.01    # Integration step, eV

epsEF = 1.0e-5

label_sample = None
xsample      = None
label_S      = None
yS           = None
label_sigma  = None
ysigma       = None
label_N      = None
yN           = None
label_mu     = None
ymu          = None

ysigmaini = None
ymuini    = None
ySini     = None

# 最適化パラメータの初期値
varname = ["meff", "l0"]
varunit = [  "me",  "m"]
ai0     = []
optid   = [     1,    1]


app = None


#=============================================
# scipy.optimize.minimizeで使うアルゴリズム
#=============================================
#nelder-mead Downhill simplex
#powell Modified Powell
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
#newton-cg Newton-CG
#trust-ncg 信頼領域 Newton-CG 法
#dogleg 信頼領域 dog-leg 法
#L-BFGS-B’ (see here)
#TNC’ (see here)
#COBYLA’ (see here)
#SLSQP’ (see here)
#trust-constr’(see here)
#dogleg’ (see here)
#trust-exact’ (see here)
#trust-krylov’ (see here)
method = "nelder-mead"
#method = 'cg'
#method = 'powell'
#method = 'bfgs'

maxiter = 1000
tol    = 1.0e-4
h_diff = 1.0e-3
outputinterval = 1

Nmin_fit     = '*'
Nmax_fit     = '*'
sigmamin_fit = '*'
sigmamax_fit = '*'

#=============================
# Graph configuration
#=============================
fig = None
figsize       = (12, 8)
figsize_sim   = (10, 6)
figsize_small = (8, 6)
figsize_FD = (8, 6)
fontsize = 18
legend_fontsize = 12
graphupdateinterval = 10


#=============================
# Treat argments
#=============================
def parameter_list():
     return mobility.meff, mobility.l0

def set_parameters(ai):
    global dos, mobility

    dos.meeff     = ai[0]
    dos.NC        = meff2NC_FEA(dos.meeff, T0)
    dos.DC0       = meff2DC0_FEA(dos.meeff, T0)
    mobility.meff = ai[0]
    if len(ai) > 1:
#        print("ai1=", ai[1])
        mobility.l0   = ai[1]
        mobility.set_scattering_parameters()
#        print("l0=", mobility.l0)

def read_parameters(path):
    global ai0, optid, T0
    
    ini = tkIniFile()
    inf = ini.ReadAll(path, AddSection = 0)
    if inf is None:
        return

    ai0 = list(ai0)
    keylist = ["meff", "l0"]
    for i in range(2):
        key = keylist[i]
        str = inf.get(key, None)
        val, id = split_optstr(str)
        if val is not None:
            ai0[i] = val
            optid[i] = id

    mobility.meff = pfloat(ai0[0])
    dos.meeff     = pfloat(ai0[0])
    mobility.l0   = pfloat(ai0[1])
    
    mobility.charge = pfloat(safe_getelement(inf, "charge", mobility.charge))
    mobility.rfac   = pfloat(safe_getelement(inf, "rfac", mobility.rfac))
    dos.EC          = pfloat(safe_getelement(inf, "EC", dos.EC))
    
    T0              = pfloat(safe_getelement(inf, "T0", T0))

def save_parameters(path, ai, args):
    ini = tkIniFile(path)
    for key in args.keys():
        ini.WriteString('Preferences', key, args[key])

    for i in range(len(ai)):
        ini.WriteString('Parameters', varname[i], "{}:{}".format(ai[i], optid[i]))

def save_parameterfile(ai = ai0, S2 = ''):
    save_parameters(parameterfile, ai, 
            {"infile": infile, "T0": T0, 
             "charge": mobility.charge, "rfac": mobility.rfac, 
             "EC": dos.EC, "NC": dos.NC, "DC0": dos.DC0, 
             "S2": S2})

def print_parameters(ai = None):
    if ai is None:
        ai = ai0
    for i in range(len(varname)):
        print("  {:10s}: {:14.8g} {:6} optid={}".format(varname[i], ai[i], varunit[i], optid[i]))

def usage(app):
    print("")
    print("Usage: Variables in () are optional")
    print(" (i) python {} basic".format(sys.argv[0]))
    print("     Plot basic functions")
    print(" (ii) python {} prop T meff r l0 kappa_latt".format(sys.argv[0]))
    print("          meff in me0, l0 in m, kappa_latt in W/m/K")
    print("     Plot Seebeck related properties")
    print("     ex: python {} basic {} {} {} {} {}".format(sys.argv[0], T0, dos.meeff, mobility.rfac, mobility.l0, klatt))
    print(" (ii') python {} Hall T meff r l0".format(sys.argv[0]))
    print("     Plot Hall related properties")
    print("     ex: python {} Hall {} {} {} {}"
                    .format(sys.argv[0], T0, dos.meeff, mobility.rfac, mobility.l0))
    print(" (iii) python {} init infile".format(sys.argv[0]))
    print("     Create .in file")
    print(" (iv) python {} sim infile T meff r l0".format(sys.argv[0]))
    print("     Calculate weighted mobility etc from input file")
    print("     Plot Jonker / Pisarenko plots")
    print("     ex: python {} sim {} {} {} {} {}".format(sys.argv[0], infile, T0, dos.meeff, mobility.rfac, mobility.l0))
    print(" (v) python {} fit T meff r l0".format(sys.argv[0]))
    print("     Fit to Jonker / Pisarenko plots")
    print("     ex: python {} fit {} {} {} {} {}"
                    .format(sys.argv[0], infile, T0, mobility.meff, mobility.rfac, mobility.l0))
    print(" (vi) python {} T infile Tmin Tmax nT".format(sys.argv[0]))
    print("     Simulate T dependences")
    print("     ex: python {} fit {} {} {} {}"
                    .format(sys.argv[0], infile, Tmin, Tmax, nT))
    print(" (vii) python {} calL infile meff r l0 ".format(sys.argv[0]))
    print("     Calculate L and kappa,e from T, Ne and sigma")
    print("     ex: python {} calL {} {} {} {}"
                    .format(sys.argv[0], infile, dos.meeff, mobility.rfac, mobility.l0))
    print(" (viii) python {} calF infile meff r l0 T_label n_label mu_label sigma_label".format(sys.argv[0]))
    print("     Calculate FHall from T, Ne and mu")
    print("     ex: python {} calF {} {} {} {} {} {} {}"
                    .format(sys.argv[0], infile, dos.meeff, mobility.rfac, mobility.l0, 
                            T_label, n_label, mu_label, sigma_label))

def updatevars():
    global mode, infile, outxlsxfile, parameterfile, parameterbkfile
    global T_label, n_label, mu_label, sigma_label, S_label
    global mobility, dos, klatt
    global T0, Tmin, Tmax, nT
    global ai0, optid
    global method, tol, maxiter, h_diff
    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit

    argv = sys.argv
#    if len(argv) == 1:
#        usage()
#        exit()

    mode   = getarg( 1, mode)
    if mode == 'basic':
        pass
    elif mode == 'prop':
        T0              = getfloatarg(2, T0)
        dos.meeff       = getfloatarg(3, dos.meeff)
        mobility.rfac   = getfloatarg(4, mobility.rfac)
        mobility.l0     = getfloatarg(5, mobility.l0)
        klatt           = getfloatarg(6, klatt)
    elif mode == 'Hall':
        T0              = getfloatarg(2, T0)
        dos.meeff       = getfloatarg(3, dos.meeff)
        mobility.rfac   = getfloatarg(4, mobility.rfac)
        mobility.l0     = getfloatarg(5, mobility.l0)
    elif mode == 'init':
        infile          = getarg     (2, infile)
    elif mode == 'sim':
        infile          = getarg     ( 2, infile)
        T_label         = getarg     ( 3, T_label)
        n_label         = getarg     ( 4, n_label)
        mu_label        = getarg     ( 5, mu_label)
        sigma_label     = getarg     ( 6, sigma_label)
        S_label         = getarg     ( 7, S_label)
        T0              = getfloatarg( 8, T0)
        dos.meeff       = getfloatarg( 9, dos.meeff)
        mobility.rfac   = getfloatarg(10, mobility.rfac)
        mobility.l0     = getfloatarg(11, mobility.l0)
    elif mode == 'fit':
        infile          = getarg     (2, infile)
        T_label         = getarg     (3, T_label)
        n_label         = getarg     (4, n_label)
        mu_label        = getarg     (5, mu_label)
        sigma_label     = getarg     (6, sigma_label)
        S_label         = getarg     (7, S_label)
        T0              = getfloatarg(8, T0)
        dos.meeff       = getfloatarg(9, dos.meeff)
        mobility.rfac   = getfloatarg(10, mobility.rfac)
        mobility.l0     = getfloatarg(11, mobility.l0)
        method          = getarg     (12, method)
        tol             = getfloatarg(13, tol)
        maxiter         = getintarg  (14, maxiter)
        Nmin_fit        = getarg     (15, Nmin_fit)
        Nmax_fit        = getarg     (16, Nmax_fit)
        sigmamin_fit    = getarg     (17, sigmamin_fit)
        sigmamax_fit    = getarg     (18, sigmamax_fit)
    elif mode == 'T':
        infile          = getarg     (2, infile)
        Tmin            = getfloatarg(3, Tmin)
        Tmax            = getfloatarg(4, Tmax)
        nT              = getintarg  (5, nT)
    elif mode == 'calL' or mode == 'calF':
        infile          = getarg     (2, infile)
        T_label         = getarg     (3, T_label)
        n_label         = getarg     (4, n_label)
        mu_label        = getarg     (5, mu_label)
        sigma_label     = getarg     (6, sigma_label)
        dos.meeff       = getfloatarg(7, dos.meeff)
        mobility.rfac   = getfloatarg(8, mobility.rfac)
        mobility.l0     = getfloatarg(9, mobility.l0)
    elif mode == 'help' or mode == 'usage':
        app.terminate("", usage = usage, pause = True)
    else:
        app.terminate("Error in updatevars(): Invalid mode {}".format(mode), usage = usage, pause = True)

    mobility.meff = dos.meeff
    mobility.set_scattering_parameters(mobility.l0, mobility.rfac)
    
    ai0 = parameter_list()


def basic():
    global mobility, dos
    
    print("")
    print("mode:", mode)

    print("")
    print("Fermi-Dirac distribution")
    xx         = []
    yfFD       = []
    yfFDh      = []
    ymdfdx     = []
    yfFDapprox = []
    print("{:8}\t{:12}\t{:12}\t{:12}\t{:12}".format("x=E/kBT", "fFDe(exact)", "fFDe(approx)", "fFDh(exact)", "-df/dx"))
    for i in range(nx_basic):
        x  = xmin_basic + i * xstep_basic
        xx.append(x)
        fFD   = 1.0 / (exp(x) + 1.0)
        fFDh  = 1.0 / (exp(-x) + 1.0)
        mdfdx = fFD * fFDh
        fFDa  = 0.5 - 0.25 * x + 1.0 / 48.0 * x**3 - 1.0 / 64.0 * x**5
        
        yfFD.append(fFD)
        yfFDh.append(fFDh)
        ymdfdx.append(mdfdx)
        yfFDapprox.append(fFDa)

        print("{:8.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}".format(x, fFD, fFDa, fFDh, mdfdx))

    print("")
    print("Fermi integrals")
    yF00 = []
    yF05 = []
    yF10 = []
    yF15 = []
    yF20 = []
    yG00 = []
    yG05 = []
    yG10 = []
    yG15 = []
    yG20 = []
    print("{:12}\t{:12}\t{:12}\t{:12}\t{:12}\t{:12}\t{:12}".format("x=EF/kBT", "F0", "F1/2", "F1", "F3/2", "F2" "exp(x)", "Gamma(1)"))
    for i in range(nx_basic):
        x  = xx[i]
        for j in range(3):
            eta = j / 2.0
            y0 = FermiIntegral(x, j)
            y1 = FermiIntegral_x2(x, j)
            validate_error(y0, y1, 2.0e-6, "Error in basic() for F{}: ".format(eta))

        yF00.append(FermiIntegral(x, 0.0))
        yF05.append(FermiIntegral(x, 0.5))
        yF10.append(FermiIntegral(x, 1.0))
        yF15.append(FermiIntegral(x, 1.5))
        yF20.append(FermiIntegral(x, 2.0))
        yG00.append(exp(x) * gamma(1.0))
        yG05.append(exp(x) * gamma(1.5))
        yG10.append(exp(x) * gamma(2.0))
        yG15.append(exp(x) * gamma(2.5))
        yG20.append(exp(x) * gamma(3.0))
        if i % 10 == 0:
            print("{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}"
                    .format(x, yF00[i], yF05[i], yF10[i], yF15[i], yF20[i], yG00[i]))

    print("")
    print("Save data to [{}]".format(outxlsfFDfile))
    tkVariousData().to_excel(outxlsfFDfile, ["x=E/kBT", "fFD(exact)", "fFD(approx)", "fFDh(exact)", "-df/dx"], 
                            [xx, yfFD, yfFDapprox, yfFDh, ymdfdx])

    print("Save data to [{}]".format(outxlsFjfile))
    tkVariousData().to_excel(outxlsFjfile, ["x=EF/kBT", "F0", "F1/2", "F1", "F3/2", "F2", "exp(x)*G(1)", 
                            "exp(x)*G(1.5)", "exp(x)*G(2)", "exp(x)*G(2.5)", "exp(x)*G(3)"], 
                            [xx, yF00, yF05, yF10, yF15, yF20, yG00, yG05, yG10, yG15, yG20])


#=============================
# グラフの表示
#=============================
    print("")

    fig = plt.figure(figsize = figsize_FD)

    ax2   = fig.add_subplot(2, 1, 1)
    ax3   = fig.add_subplot(2, 1, 2)

    ax2.plot(xx, yfFD,       label = '$f_{FD,e}(exact)$', linestyle = '-',      color = 'black', linewidth = 0.5)
    ax2.plot(xx, yfFDh,      label = '$f_{FD,h}(exact)$', linestyle = '-',      color = 'blue', linewidth = 0.5)
    ax2.plot(xx, yfFDapprox, label = '$f_{FD}(approx)$',  linestyle = '-',      color = 'green', linewidth = 0.5)
    ax2.plot(xx, ymdfdx,     label = '$-df/dx$',          linestyle = 'dashed', color = 'red', linewidth = 0.5)
    ax2.set_xlabel("$x=(E-E_F)/k_BT$ (eV)", fontsize = fontsize)
    ax2.set_ylabel("$f_{FD}$, $-df/dx$", fontsize = fontsize)
    ax2.set_xlim([-10.0, 10.0])
    ax2.set_ylim([-0.1, 1.1])
    ax2.legend(fontsize = legend_fontsize)
    ax2.tick_params(labelsize = fontsize)

    ax3.plot(xx, yF00, label = '$F_0$',     linestyle = '-', color = 'black', linewidth = 0.5)
    ax3.plot(xx, yF05, label = '$F_{1/2}$', linestyle = '-', color = 'red',   linewidth = 0.5)
    ax3.plot(xx, yF10, label = '$F_1$',     linestyle = '-', color = 'blue',  linewidth = 0.5)
    ax3.plot(xx, yF15, label = '$F_{3/2}$', linestyle = '-', color = 'purple', linewidth = 0.5)
    ax3.plot(xx, yF20, label = '$F_2$',     linestyle = '-', color = 'green', linewidth = 0.5)
    ax3.plot(xx, yG00, label = r'$e^x$$\Gamma(1)$',   linestyle = '', marker = 'o', markerfacecolor = 'black',  markersize = 1)
    ax3.plot(xx, yG05, label = r'$e^x$$\Gamma(3/2)$', linestyle = '', marker = 'o', markerfacecolor = 'red',    markersize = 1)
    ax3.plot(xx, yG10, label = r'$e^x$$\Gamma(2)$',   linestyle = '', marker = 'o', markerfacecolor = 'blue',   markersize = 1)
    ax3.plot(xx, yG15, label = r'$e^x$$\Gamma(5/2)$', linestyle = '', marker = 'o', markerfacecolor = 'purple', markersize = 1)
    ax3.plot(xx, yG20, label = r'$e^x$$\Gamma(3)$',   linestyle = '', marker = 'o', markerfacecolor = 'green',  markersize = 1)
    ax3.set_xlabel("$x=E_F/k_BT$ (eV)", fontsize = fontsize)
    ax3.set_ylabel("$F_r$", fontsize = fontsize)
    ax3.set_xlim([-10.0, 10.0])
    ax3.set_ylim([1.0e-5, 0.1e4])
    ax3.set_yscale('log')
    ax3.legend(fontsize = legend_fontsize)
    ax3.tick_params(labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("", pause = True)

def Hall():
    print("")
    print("Carrier:")
    print("  q={} e".format(mobility.charge))
    print("  meff={}me".format(dos.meeff))
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    if mobility.charge < 0.0:
        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    else:
        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
    print("  mean free path prefactor: l0={} m".format(mobility.l0))

    print("")
    print(f"Calculat Hall coefficient, Hall factor, and Hall mobility at {T0} K")
    xx         = []
    yEF        = []
    yNe        = []
    ysigma     = []
    ymu        = []
    ytau       = []
    yRH0       = []
    yRH        = []
    yFH        = []
    yNe_Hall    = []
    ymu_Hall   = []
    print("{:>8} {:>8} {:>12} {:>12} {:>12} {:>12} {:>12}"
            .format("x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "FH", "RH(m^3/C"))
    for i in range(nx):
        x  = xmin + i * xstep
        EF = x * kB * T0 / e
        
        inf = dos.cal_Hall_properteis(T0, EF, mobility, validate_error_str = 'Error in properteis()')

        sigma   = inf["sigma"]
        Ne      = inf["n"]
        mu      = inf["mu"]
        tau     = inf["<tau>"]
        FH      = inf["FH"]
        RH      = inf["RH"]
        RH0     = inf["RH0"]
        n_Hall  = inf["nHall"]
        mu_Hall = inf["muHall"]

        xx.append(x)
        yEF.append(EF)
        yNe.append(Ne)
        ysigma.append(sigma)
        ymu.append(mu)
        ytau.append(FH)
        yRH0.append(RH0)
        yRH.append(RH)
        yFH.append(FH)
        yNe_Hall.append(n_Hall)
        ymu_Hall.append(mu_Hall)

        if i % 10 == 0:
            print("{:8.3g} {:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
                    .format(x, EF, yNe[i], ysigma[i], ymu[i], yFH[i], yRH[i]))

    print("")
    print("Save data to [{}]".format(outxlsPropfile))
    tkVariousData().to_excel(outxlsPropfile, ["x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "FH", "RH(m^3/C"], 
                             [xx, yEF, yNe, ysigma, ymu, yFH, yRH])


#=============================
# グラフの表示
#=============================
    print("")
    fig = plt.figure(figsize = figsize)

    axtau   = fig.add_subplot(3, 3, 1)
    axNe    = fig.add_subplot(3, 3, 2)
    axsigma = fig.add_subplot(3, 3, 3)
    axmu    = fig.add_subplot(3, 3, 4)
    aRH     = fig.add_subplot(3, 3, 5)
    aFH     = fig.add_subplot(3, 3, 6)
    aNeRH   = fig.add_subplot(3, 3, 7)
    aNeFH   = fig.add_subplot(3, 3, 8)

    axsigma.set_title("$T_0$={} K $m^*$={}$m_e$ r={} $k_l$$_a$$_t$$_t$={} W/m/K".format(T0, dos.meeff, mobility.rfac, klatt))

    axtau.plot(xx, ytau,    label = '$\\tau(EF)$',  linestyle = '-', color = 'red',  linewidth = 0.5)
#    axtau.plot(xx,     ytauavg, label = '<$\\tau$>',    linestyle = '-', color = 'blue', linewidth = 0.5)
    axtau.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axtau.set_ylabel("$\\tau$ (fs)", fontsize = fontsize)
    axtau.set_ylim([0.0, max(ytau) * 1.1])
    axtau.legend(fontsize = legend_fontsize, loc = 'best')
    axtau.tick_params(labelsize = fontsize)

    axNe.plot(xx, yNe,     label = '$N_e$',           linestyle = '-', color = 'red',  linewidth = 1.0)
    axNe.plot(xx, yNe_Hall, label = '$N_e$$_{,Hall}$', linestyle = '-', color = 'blue',  linewidth = 1.0)
    axNe.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axNe.set_yscale('log')
    axNe.legend(fontsize = legend_fontsize, loc = 'best')
    axNe.tick_params(labelsize = fontsize)

#    axsigma.plot(xx, ysigma, label = '$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
    axsigma.plot(yNe_Hall, ysigma, label = r'$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
#    axsigma.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axsigma.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
    axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
    axsigma.set_xscale('log')
    axsigma.set_yscale('log')
    axsigma.set_xlim([1.0e10, 1.0e23])
    axsigma.legend(fontsize = legend_fontsize, loc = 'best')
    axsigma.tick_params(labelsize = fontsize)

#    axmu.plot(xx, ymu, label = r'$\\mu$',                linestyle = '-', color = 'red',  linewidth = 1.0)
    axmu.plot(yNe_Hall, ymu, label = r'$\\mu$',                linestyle = '-', color = 'red',  linewidth = 1.0)
#    axmu.plot(xx, ymu_Hall, label = r'$\\mu$$_{,Hall}$', linestyle = '-', color = 'blue',  linewidth = 1.0)
    axmu.plot(yNe_Hall, ymu_Hall, label = r'$\\mu$$_{,Hall}$', linestyle = '-', color = 'blue',  linewidth = 1.0)
#    axmu.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axmu.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
    axmu.set_ylabel("$\\mu$ (cm$^2$/V/s)", fontsize = fontsize)
    axmu.set_xscale('log')
    axmu.set_xlim([1.0e10, 1.0e23])
    axmu.legend(fontsize = legend_fontsize, loc = 'best')
    axmu.tick_params(labelsize = fontsize)

#    aRH.plot(xx, yRH0, label = r'$R_H$$_0$=$1/e/N_e$', linestyle = '-', color = 'red',  linewidth = 1.0)
    aRH.plot(yNe_Hall, yRH0, label = r'$R_H$$_0$=$1/e/N_e$', linestyle = '-', color = 'red',  linewidth = 1.0)
#    aRH.plot(xx, yRH,  label = r'$R_H$$_{,Hall}$',     linestyle = '-', color = 'blue',  linewidth = 1.0)
    aRH.plot(yNe_Hall, yRH,  label = r'$R_H$$_{,Hall}$',     linestyle = '-', color = 'blue',  linewidth = 1.0)
#    aRH.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    aRH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
    aRH.set_ylabel(r"$R_H$ (cm$^3$/C)", fontsize = fontsize)
    aRH.set_xscale('log')
    aRH.set_yscale('log')
    aRH.set_xlim([1.0e10, 1.0e23])
    aRH.legend(fontsize = legend_fontsize, loc = 'best')
    aRH.tick_params(labelsize = fontsize)

#    aFH.plot(xx, yFH, label = r'$F_{Hall}$', linestyle = '-', color = 'red',  linewidth = 1.0)
    aFH.plot(yNe_Hall, yFH, label = r'$F_{Hall}$', linestyle = '-', color = 'red',  linewidth = 1.0)
#    aFH.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    aFH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
    aFH.set_ylabel("$F_H$", fontsize = fontsize)
    aFH.set_xscale('log')
    aFH.set_xlim([1.0e10, 1.0e23])
    aFH.legend(fontsize = legend_fontsize, loc = 'best')
    aFH.tick_params(labelsize = fontsize)

    aNeRH.plot(yNe_Hall, yRH, label = r'$R_H$$_{,Hall}$', linestyle = '-', color = 'red',  linewidth = 1.0)
    aNeRH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
    aNeRH.set_ylabel(r"$R_H$ (cm$^3$/C)", fontsize = fontsize)
    aNeRH.set_xscale('log')
    aNeRH.set_yscale('log')
    aNeRH.set_xlim([1.0e10, 1.0e23])
#    aNeRH.set_ylim([min(yRH) * 0.5, min(yRH) * 1.0e3])
#    aNeRH.legend(fontsize = legend_fontsize, loc = 'best')
    aNeRH.tick_params(labelsize = fontsize)

    aNeFH.plot(yNe_Hall, yFH, label = r'$N_e$', linestyle = '-', color = 'red',  linewidth = 1.0)
    aNeFH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
    aNeFH.set_ylabel(r"$F_H$ (cm$^{-3}$)", fontsize = fontsize)
    aNeFH.set_xscale('log')
    aNeFH.set_xlim([1.0e10, 1.0e23])
#    aNeFH.legend(fontsize = legend_fontsize, loc = 'best')
    aNeFH.tick_params(labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("", pause = True)

def properties():
    global mobility, dos
    
    print("")
    print("mode:", mode)
    print("T={} K".format(T0))
    print("meff:", dos.meeff)
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    if mobility.charge < 0.0:
        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    else:
        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor in tau(E) prop to E^(r-0.5) (non-degenerated) r={}".format(mobility.rfac))
    print("Lorentz number")
    print("  Free electron mode: {:12.8g} Wohm/K^2".format(L_FEA))
    print("Mean free path prefactor: l0={} m".format(mobility.l0))
    print("Lattice thermal conductivity: klatt={} W/K".format(klatt))

    print("")
    print("Properties")
    xx         = []
    xx_tau     = []
    yEF        = []
    yNe        = []
    ysigma     = []
    ymu        = []
    ytau       = []
    ytauavg    = []
    yS         = []
    ySndeg     = []
    ySdeg      = []
    ykappa     = []
    ykappatot  = []
    yL         = []
    yLapprox   = []
    yPF        = []
    yZT        = []
    print("{:>8} {:>8} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}"
            .format("x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "<tau>(fs)",
                    "S(uV/K)", "Snon-deg", "Sdeg", "kappa,e(W/m/K)", "kappa,tot", "L(Wohm/K^2)", "Lapprox", "PF(W/m/K^2)", "ZT"))
    for i in range(nx):
        x  = xmin + i * xstep
        EF = x * kB * T0 / e
        
#        print("T=", T0)
        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
            = dos.cal_transport_properteis(T0, EF, mobility, klatt, validate_error_str = 'Error in properteis()')

# Relaxation time tau(EF)
        if EF <= 0.0:
            tau = None
        else:
            tau = mobility.tau(T0, EF) * 1.0e15      # fs, E in eV
        if EF > 0.0:
            xx_tau.append(x)
            ytau.append(tau)

        Sndeg = dos.cal_S_nondegenerated_from_Ne(n, mobility.rfac, charge = mobility.charge) * 1.0e6    # uV/K
        Sdeg  = dos.cal_S_degenerated_from_Ne(T0, n, mobility.rfac, charge = mobility.charge) * 1.0e6  # uV/K

# Approximated Lorentz number
        Lapprox = dos.cal_LorentzNumber_from_S_approximate(S)

        xx.append(x)
        yEF.append(EF)
        yNe.append(n)
        ysigma.append(sigma)
        ymu.append(mu)
        ytauavg.append(tau_avg)
        yS.append(S)
        ySndeg.append(Sndeg)
        ySdeg.append(Sdeg)
        ykappa.append(kappa)
        ykappatot.append(kappa + klatt)
        yL.append(L)
        yLapprox.append(Lapprox)
        yPF.append(PF)
        yZT.append(ZT)
        if i % 10 == 0:
            print("{:8.3g} {:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
                    .format(x, EF, n, sigma, mu, tau_avg, S, Sndeg, Sdeg, kappa, kappa_tot, L, Lapprox, PF, ZT))

    print("")
    print("Save data to [{}]".format(outxlsPropfile))
    tkVariousData().to_excel(outxlsPropfile, ["x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "<tau>(fs)",
                    "S(uV/K)", "Snon-deg", "Sdeg", "kappa,e(W/m/K)", "kappa,tot(W/m/K)", "L(Wohm/K^2)", "Lapprox", "PF(uW/cm/K^2)", "ZT"], 
                             [xx, yEF, yNe, ysigma, ymu, ytauavg, yS, ySndeg, ySdeg, ykappa, ykappatot, yL, yLapprox, yPF, yZT])


#=============================
# グラフの表示
#=============================
    print("")
    maxS = max(yS)
    minS = min(yS)
    if maxS > 0.0:
        maxS *= 1.2
    else:
        maxS = 0.0
    if minS > 0.0:
        minS = 0.0
    else:
        minS *= 2.0

    fig = plt.figure(figsize = figsize)

    axtau   = fig.add_subplot(3, 4, 1)
    axNe    = fig.add_subplot(3, 4, 2)
    axsigma = fig.add_subplot(3, 4, 3)
    axmu    = fig.add_subplot(3, 4, 4)
    axS     = fig.add_subplot(3, 4, 5)
    axke    = fig.add_subplot(3, 4, 6)
    axPF    = fig.add_subplot(3, 4, 7)
    axZT    = fig.add_subplot(3, 4, 8)
    axL1    = fig.add_subplot(3, 4, 9)
    axL2    = fig.add_subplot(3, 4, 10)
    axLS    = fig.add_subplot(3, 4, 11)

    axsigma.set_title(r"$T_0$={} K $m^*$={}$m_e$ r={} $k_l$$_a$$_t$$_t$={} W/m/K".format(T0, dos.meeff, mobility.rfac, klatt))

    axtau.plot(xx_tau, ytau,    label = '$\\tau(EF)$',  linestyle = '-', color = 'red',  linewidth = 0.5)
    axtau.plot(xx,     ytauavg, label = '<$\\tau$>',    linestyle = '-', color = 'blue', linewidth = 0.5)
    axtau.set_xlabel("$E$ (eV)", fontsize = fontsize)
    axtau.set_ylabel("$\\tau$ (fs)", fontsize = fontsize)
    axtau.set_ylim([0.0, max(ytau) * 1.1])
    axtau.legend(fontsize = legend_fontsize, loc = 'best')
    axtau.tick_params(labelsize = fontsize)

    axNe.plot(xx, yNe, label = '$N_e$',    linestyle = '-', color = 'red',  linewidth = 1.0)
    axNe.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axNe.set_yscale('log')
    axNe.legend(fontsize = legend_fontsize, loc = 'best')
    axNe.tick_params(labelsize = fontsize)

#    axsigma.plot(xx, ysigma, label = r'$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
    axsigma.plot(yNe, ysigma, label = r'$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
#    axsigma.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axsigma.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
    axsigma.set_xscale('log')
    axsigma.set_yscale('log')
    axsigma.legend(fontsize = legend_fontsize, loc = 'best')
    axsigma.tick_params(labelsize = fontsize)

#    axmu.plot(xx, ymu, label = '$\\mu$',    linestyle = '-', color = 'red',  linewidth = 1.0)
    axmu.plot(yNe, ymu, label = '$\\mu$',    linestyle = '-', color = 'red',  linewidth = 1.0)
#    axmu.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axmu.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axmu.set_ylabel("$\\mu$ (cm$^2$/V/s)", fontsize = fontsize)
    axmu.set_xscale('log')
    axmu.legend(fontsize = legend_fontsize, loc = 'best')
    axmu.tick_params(labelsize = fontsize)

    axke.plot(yNe, ykappa, label = '$\\kappa$$_e$',    linestyle = '-', color = 'red',  linewidth = 1.0)
    axke.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
#    axke.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
    axke.set_ylabel("$\\kappa$$_e$ (W/m/K)", fontsize = fontsize)
    axke.set_xscale('log')
    axke.set_yscale('log')
    axke.legend(fontsize = legend_fontsize, loc = 'best')
    axke.tick_params(labelsize = fontsize)

    axS.plot(yNe, yS,     label = '$S$',           linestyle = '-',      color = 'black', linewidth = 1.0)
    axS.plot(yNe, ySndeg, label = '$S_{non-deg}$', linestyle = 'dashed', color = 'red',   linewidth = 0.5)
    axS.plot(yNe, ySdeg,  label = '$S_{deg}$',     linestyle = 'dashed', color = 'blue',  linewidth = 0.5)
    axS.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axS.set_ylabel(r"$S$ ($\mu$V/K)", fontsize = fontsize)
    axS.set_xscale('log')
    axS.set_ylim([minS, maxS])
    axS.legend(fontsize = legend_fontsize, loc = 'best')
    axS.tick_params(labelsize = fontsize)

    axPF.plot(yNe, yPF, label = '$PF$',           linestyle = '-',      color = 'black', linewidth = 1.0)
    axPF.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axPF.set_ylabel(r"$PF$ ($\mu$W/m/K$^2$)", fontsize = fontsize)
    axPF.set_xscale('log')
    axPF.legend(fontsize = legend_fontsize, loc = 'best')
    axPF.tick_params(labelsize = fontsize)

    axZT.plot(yNe, yZT, label = '$ZT$',           linestyle = '-',      color = 'black', linewidth = 1.0)
    axZT.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axZT.set_ylabel("$ZT$", fontsize = fontsize)
    axZT.set_xscale('log')
    axZT.legend(fontsize = legend_fontsize, loc = 'best')
    axZT.tick_params(labelsize = fontsize)

#    axL1.plot(xx, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
    axL1.plot(yNe, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
    axL1.plot([min(xx), max(xx)], [L_FEA, L_FEA], label = 'L(FEA)', linestyle = '-', color = 'red', linewidth = 0.5)
#    axL1.set_xlabel("$x=(E_F-E_C)/k_B/T$", fontsize = fontsize)
    axL1.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axL1.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
    axL1.set_ylim([0.0, 3.2e-8])
    axL1.legend(fontsize = legend_fontsize, loc = 'best')
    axL1.tick_params(labelsize = fontsize)

    axL2.plot(yNe, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
    axL2.plot([min(yNe), max(yNe)], [L_FEA, L_FEA], label = 'L(FEA)', linestyle = '-', color = 'red', linewidth = 0.5)
    axL2.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axL2.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
    axL2.set_xscale('log')
    axL2.set_xlim([1.0e10, 1.0e23])
    axL2.set_ylim([0.0, 3.2e-8])
    axL2.legend(fontsize = legend_fontsize, loc = 'best')
    axL2.tick_params(labelsize = fontsize)

    axLS.plot(yS, yL,       label = 'L(cal)',         linestyle = '-', color = 'red',  linewidth = 1.0)
    axLS.plot(yS, yLapprox, label = 'L(approx, r=0)', linestyle = '-', color = 'blue', linewidth = 0.5)
    axLS.plot([min(yS), max(yS)], [L_FEA, L_FEA], label = 'L(FEA)', linestyle = '-', color = 'red', linewidth = 0.5)
    axLS.set_xlabel(r"$S$ ($\mu$V/K)", fontsize = fontsize)
    axLS.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
    axLS.set_ylim([0.0, 3.2e-8])
    axLS.legend(fontsize = legend_fontsize, loc = 'best')
    axLS.tick_params(labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("", pause = True)

def recover_parameters(ais, optid, aidef0):
    aidef = list(aidef0).copy()
    
    ai = []
    c = 0
    for i in range(len(optid)):
        if optid[i] == 1:
            ai.append(ais[c])
            c += 1
        else:
            ai.append(aidef[i])

    return ai

def choose_parameters(ais0, optid):
    ais = list(ais0).copy()

    ai = []
    c = 0
    for i in range(len(optid)):
        if optid[i] == 1:
            ai.append(ais[i])

    return ai

def calS2_sigma(ai):
    global mu, dos, ysigma

    aiorg = parameter_list()

    ainew = [ai[0]]
    if ainew[0] <= 0.0:
        ainew[0] = 1.0e-4

    set_parameters([mobility.meff, ainew[0]])

    S2 = 0.0
    ndata = len(ysigma)
    eps = 1.0e-300
    for i in range(ndata):
        S      = yS[i]
        sigma  = ysigma[i]
        n      = yN[i]
        mu     = ymu[i]

        if Nmin_fit is not None and n < Nmin_fit:
            continue
        if Nmax_fit is not None and Nmax_fit < n:
            continue
        if sigmamin_fit is not None and sigma < sigmamin_fit:
            continue
        if sigmamax_fit is not None and sigmamax_fit < sigma:
            continue

        EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
        if not flag:
            terminate("Error in callback(): EF calculation did not converge. diffEF={}  eps={} in {} iteration"
                    .format(diffEF, epsEF, maxiter))

        sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
            = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
#        print("m,l0=", mobility.meff, mobility.l0)
#        print("  sigma=", i, ysigma[i], sigmafin)

        dlogsigma = ysigma[i] - sigmafin
#        dlogsigma = log(ysigma[i]+eps) - log(sigmafin+eps)
        S2 += dlogsigma * dlogsigma

    S2 /= (ndata - 1)

    set_parameters(aiorg)
#    print("S2=", S2)
    
    return S2

def calS2_S(ai):
    global mu, dos, ysigma
    global method, tol, maxiter, h_diff
    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit

    aiorg = parameter_list()

    ainew = [ai[0]]
    if ainew[0] <= 0.0:
        ainew[0] = 1.0e-4

    set_parameters(ainew)
#    print_parameters()
#    print("m=", dos.meeff, dos.NC, dos.DC0)

    S2 = 0.0
    ndata = len(ysigma)
    eps = 1.0e-300
    for i in range(ndata):
        S      = yS[i]
        sigma  = ysigma[i]
        n      = yN[i]
        mu     = ymu[i]

        if Nmin_fit is not None and n < Nmin_fit:
            continue
        if Nmax_fit is not None and Nmax_fit < n:
            continue
        if sigmamin_fit is not None and sigma < sigmamin_fit:
            continue
        if sigmamax_fit is not None and sigmamax_fit < sigma:
            continue

        EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
        if not flag:
            app.terminate("Error in callback(): EF calculation did not converge. diffEF={}  eps={} in {} iteration"
                    .format(diffEF, epsEF, maxiter), usage = usage, pause = True)

        sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
            = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
        
        dlogS = yS[i] - Sfin
        """
        if yS[i] <= 0.0:
            dlogS = log(-yS[i]+eps) - log(-Sfin+eps)
        else:
            dlogS = log(yS[i]+eps) - log(Sfin+eps)
        """
        S2 += dlogS * dlogS

    S2 /= (ndata - 1)

    set_parameters(aiorg)
    
    return S2

# First derivatives to be used e.g. for cg method
# Approximate by forward difference method with the delta h = h_diff
def diff1_sigma(ai):
    n = len(ai)
    f0 = calS2_sigma(ai)
    df = np.empty(n)
    for i in range(n):
        aii = ai
        aii[i] = ai[i] + h_diff
        df[i] = (calS2_sigma(aii) - f0) / h_diff
    return df

def diff1_S(ai):
    n = len(ai)
    f0 = calS2_S(ai)
    df = np.empty(n)
    for i in range(n):
        aii = ai
        aii[i] = ai[i] + h_diff
        df[i] = (calS2_S(aii) - f0) / h_diff
    return df

def plot(fig, yn, ysigma, ymu, yS, ysigmaini, ymuini, ySini, 
                  ysigmafin = None, ymufin = None, ySfin = None):
    global plt
    global graphupdateinterval

    plt.clf()

    ax1  = fig.add_subplot(1, 2, 1)
    ax2  = fig.add_subplot(1, 2, 2)
    plt.title(infile, fontsize = fontsize)

    Srange = [min(yS) * 0.9, max(yS) * 1.1]
    yn2, ySi2 = sort_lists([yn, ySini])
    ins1 = ax1.plot(yn, yS,       label = 'Pisarenko(obs)', linestyle = 'none', marker = 'o', markerfacecolor = 'blue', markeredgecolor = 'blue')
    ins2 = ax1.plot(yn2, ySi2,    label = 'Pisarenko(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
#    ins2 = ax1.plot(yn, ySini,    label = 'Pisarenko(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
    if ysigmafin:
        yn2, ySf2 = sort_lists([yn, ySini])
        ins3 = ax1.plot(yn2, ySf2,    label = 'Pisarenko(fin)', linestyle = '-',      color = 'red')
#        ins3 = ax1.plot(yn, ySfin,    label = 'Pisarenko(fin)', linestyle = '-',      color = 'red')
    ax1.set_xscale('log')
    ax1.set_ylim([0.0, Srange[1]])
    ax1.set_xlabel(r'$N_e$ (cm$^{-3}$)', fontsize = fontsize)
    ax1.set_ylabel(r'S ($\mu$V/K)', fontsize = fontsize)
    ax1.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
    ax1.tick_params(labelsize = fontsize)

    ysi, ySi = sort_lists([ysigmaini, ySini])
    ins1 = ax2.plot(ysigma,    yS,       label = 'Jonker(obs)', linestyle = 'none', marker = 'o', markerfacecolor = 'blue', markeredgecolor = 'blue')
    ins2 = ax2.plot(ysi, ySi,    label = 'Jonker(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
#    ins2 = ax2.plot(ysigmaini, ySini,    label = 'Jonker(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
    if ysigmafin:
        ysf, ySf = sort_lists([ysigmafin, ySfin])
        ins3 = ax2.plot(ysf, ySf,    label = 'Jonker(fin)', linestyle = '-',      color = 'red')
#        ins3 = ax2.plot(ysigmafin, ySfin,    label = 'Jonker(fin)', linestyle = '-',      color = 'red')
    ax2.set_xscale('log')
    ax2.set_ylim([0.0, Srange[1]])
    ax2.set_xlabel(r'$\sigma$ (S/cm)', fontsize = fontsize)
    ax2.set_ylabel(r'S ($\mu$V/K)', fontsize = fontsize)
    ax2.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
    ax2.tick_params(labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.01)


# Callback function for scipy.optimize.minimize()
# Print variables every iteration, and update graph for every graphupdateinterval iterationsおｔ
iter = 0
def callback(fig, S2):
    global iter, graphupdateinterval
    global outputinterval
    global dos, mu
    global ysigma
    global method, tol, maxiter, h_diff
    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit

    if iter % outputinterval == 0:
        print("callback {:04d}: {:10.4g} {:10.4g}  S2={:10.6g}".format(iter, mobility.meff, mobility.l0, S2))

    if iter % graphupdateinterval == 0:
        ndata = len(ysigma)
        ysigmafin = []
        ymufin    = []
        ySfin     = []
        for i in range(ndata):
            S      = yS[i]
            sigma  = ysigma[i]
            n      = yN[i]
            mu     = ymu[i]

            if Nmin_fit is not None and n < Nmin_fit:
                continue
            if Nmax_fit is not None and Nmax_fit < n:
                continue
            if sigmamin_fit is not None and sigma < sigmamin_fit:
                continue
            if sigmamax_fit is not None and sigmamax_fit < sigma:
                continue

            EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
            if not flag:
                app.terminate("Error in callback(): EF calculation did not converge. diffEF={}  eps={} in {} iteration"
                        .format(diffEF, epsEF, maxiter), usage = usage, pause = True)

            sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
                = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
        
            ysigmafin.append(sigmafin)
            ymufin.append(mufin)
            ySfin.append(Sfin)

        plot(fig, yN, ysigma, ymu, yS, ysigmaini, ymuini, ySini, ysigmafin, ymufin, ySfin)

    iter += 1

def callback_S(fig, xk):
    S2 = calS2_S(xk)
    ainew = [xk[0]]
    set_parameters(ainew)

    callback(fig, S2)

def callback_sigma(fig, xk):
    global iter, graphupdateinterval
    global outputinterval
    global dos, mu
    global ysigma
    
    S2 = calS2_sigma(xk)
    ainew = [mobility.meff, xk[0]]
    set_parameters(ainew)

    callback(fig, S2)

def construct_lists(label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu):
    ndata = len(yS)
    if label_sample is None:
        label_sample = 'sample'
        xsample = []
        for i in range(ndata):
            xsample.append(i + 1)
    if label_sigma is None:
        ysigma = []
        for i in range(ndata):
            ysigma.append(None)
    if label_N is None:
        yN = []
        for i in range(ndata):
            if ymu is not None:
                yN.append(None)
            else:
                yN.append(ysigma[i] / e / ymu[i])
    if label_mu is None:
        ymu = []
        for i in range(ndata):
            if yN[i] is not None:
                ymu.append(None)
            else:
                ymu.append(ysigma[i] / e / yN[i])

    return label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu 

def fit():
    global fig
    global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
    global ysigmaini, ymuini, ySini
    global method, tol, maxiter, h_diff
    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
    
    print("infile     :", infile)
    print("  T_label    : ", T_label)
    print("  S_label    : ", S_label)
    print("  n_label    : ", n_label)
    print("  mu_label   : ", mu_label)
    print("  sigma_label: ", sigma_label)
    print("")
    print(f"T0: {T0} K")
    print("Carrier:")
    print("  q={} e".format(mobility.charge))
    print("  meff={}me".format(dos.meeff))
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
    print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
    print("  mean free path prefactor: l0={} m".format(mobility.l0))
    print("")
    print( "Fitting condition")
    print(f"  method : {method}")
    print(f"  tol    : {tol}")
    print(f"  maxiter: {maxiter}")
    print(f"  h_diff: {h_diff}")
    print(f"  outputinterval: {outputinterval}")
    print(f"  N range    : {Nmin_fit} - {Nmax_fit} cm^-3")
    print(f"  sigma range: {sigmamin_fit} - {sigmamax_fit} S/cm")
    def set_default_range(var):
        if var == '' or var == '*':
            return None
        return pfloat(var, defval = None)

    Nmin_fit = set_default_range(Nmin_fit)
    Nmax_fit = set_default_range(Nmax_fit)
    sigmamin_fit = set_default_range(sigmamin_fit)
    sigmamax_fit = set_default_range(sigmamax_fit)

    if '***' in method:
        app.terminate(f"***Error: Choose method", pause = True)

    if '***' in S_label:
        app.terminate(f"***Error: Choose S_label", pause = True)
    if '***' in n_label:
        app.terminate(f"***Error: Choose n_label", pause = True)
    if '***' in mu_label:
        app.terminate(f"***Error: Choose mu_label", pause = True)

    if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
        app.terminate(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]", pause = True)

    print("")
    print("Read S data from {}".format(infile))

    datafile = tkVariousData(infile)
    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
    label_S, yS           = datafile.FindDataArray(S_label, flag = 'i')
# Convert V/K to uV/K
    for i in range(len(yS)):
        yS[i] *= 1.0e6 
    label_sigma, ysigma   = datafile.FindDataArray(sigma_label, flag = 'i')
    label_N, yN           = datafile.FindDataArray(n_label, flag = 'i')
    label_mu, ymu         = datafile.FindDataArray(mu_label, flag = 'i')
    label_sample, xsample = labels[0], datalist[0]

    ndata = len(yS)
    print("ndata(all): ", ndata)
    ndata_fit = 0
    print("data to be fitted")
    print(f"{'sample':15}  {'n(cm-3)':12}  {'mu(cm2/Vs)':10}  {'sigma(S/cm)':10}  {'S(uV/K)':10}")
    for i in range(ndata):
        sigma = ysigma[i]
        n = yN[i]
        if Nmin_fit is not None and n < Nmin_fit:
            continue
        if Nmax_fit is not None and Nmax_fit < n:
            continue
        if sigmamin_fit is not None and sigma < sigmamin_fit:
            continue
        if sigmamax_fit is not None and sigmamax_fit < sigma:
            continue

        print(f"{xsample[i]:15}  {n:12.4g}  {ymu[i]:10.4g}  {sigma:10.4g}  {yS[i]:10.4g}")

    print("")
    print("Calculate initial values")
    print_parameters()
    yEFini    = []
    ysigmaini = []
    ymuini    = []
    ySini     = []
    print("{:>8} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}"
            .format("EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "Ne,ini", "sigma,ini", "mu,ini", "S,ini"))
    for i in range(ndata):
        sample = xsample[i]
        S      = yS[i]
        sigma  = ysigma[i]
        n      = yN[i]
        mu     = ymu[i]

        EFini, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
        if not flag:
            app.terminate("Error in fit(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
                        usage = usage, pause = True)

#        x = (EF - doc.EC) * e / kB / T0
        sigmaini, nini, muini, tau_avg, Sini, kappa, kappa_tot, L, PF, ZT, inf = \
                dos.cal_transport_properteis(T0, EFini, mobility, klatt, validate_error_str = 'Error in fit()')
        
        yEFini.append(EFini)
        ysigmaini.append(sigmaini)
        ymuini.append(muini)
        ySini.append(Sini)

        print("{:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
                    .format(EFini, n, sigma, mu, S, nini, sigmaini, muini, Sini))

#=============================
# グラフの表示
#=============================
    print("")
    print("plot")
    fig = plt.figure(figsize = figsize_sim)

    plot(fig, yN, ysigma, ymu, yS, ysigmaini, ymuini, ySini)

#=============================
# Optimization
#=============================
    print("")
    print("Nonlinear least-squares fitting for effective mass by ", method)
    print("  tol=", tol)
    ai = [mobility.meff]
    ret = minimize(calS2_S, ai, method = method, jac = diff1_S, tol = tol, callback = lambda xk: callback_S(fig, xk),
                options = {'maxiter':maxiter, "disp":True})
    if method == 'nelder-mead':
        simplex = ret['final_simplex']
        ai = simplex[0][0]    
        fmin = ret['fun']
    aifin = [ai[0], mobility.l0] #ai0[1]]
    set_parameters(aifin)
    print("Optimized at S2={:12.6g}:".format(fmin))
    print_parameters(parameter_list())

    print("")
    print("Nonlinear least-squares fitting for l0 by ", method)
    print("  tol=", tol)
    ai_sigma = [mobility.l0]
    print("meff(ini)=", mobility.meff)
#    ret = minimize(calS2_S, ai_sigma, method = method, jac = diff1_sigma, tol = tol, callback = lambda xk: callback_sigma(fig, xk),
    ret = minimize(calS2_sigma, ai_sigma, method = method, jac = diff1_sigma, tol = tol, callback = lambda xk: callback_sigma(fig, xk),
                options = {'maxiter':maxiter, "disp":True})
    if method == 'nelder-mead':
        simplex = ret['final_simplex']
        ai_sigma = simplex[0][0]    
        fmin = ret['fun']
    aifin = [mobility.meff, ai_sigma[0]]
    set_parameters(aifin)
    print("Optimized at S2={:12.6g}:".format(fmin))
#    print_parameters()
    print_parameters(parameter_list())
    
    print("")
    print("Calculate final values")
    print_parameters(parameter_list())
    yEFfin    = []
    ysigmafin = []
    ySfin     = []
    print("{:>8} {:>12} {:>12} {:>12} {:>12} {:>12}"
            .format("EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "S,fin"))
    for i in range(ndata):
        sample = xsample[i]
        S      = yS[i]
        sigma  = ysigma[i]
        n      = yN[i]
        mu     = ymu[i]

        EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
        if not flag:
            app.terminate("Error in fit(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
                        usage = usage, pause = True)

        sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
            = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
        
        mufin = sigmafin / sigma / mufin
        sigmafin = e * nfin * mufin
        
        yEFfin.append(EFfin)
        ysigmafin.append(sigmafin)
        ySfin.append(Sfin)
        print("{:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
                    .format(EFfin, n, sigma, mu, S, Sfin))

    print("")
    print("Save final parameters to [{}]".format(parameterfile))
    save_parameterfile(ai = aifin, S2 = fmin)

    print("Save fitting data to [{}]".format(outxlsfile))
    tkVariousData().to_excel(outxlsfile, ["EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "S,fin"], 
            [yEFfin, yN, ysigma, ymu, yS, ySfin])

    app.terminate("", pause = True)

def sim():
    global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
    global Nmin, Nmax
    
    print("")
    print("infile     :", infile)
    print("  T_label    : ", T_label)
    print("  S_label    : ", S_label)
    print("  n_label    : ", n_label)
    print("  mu_label   : ", mu_label)
    print("  sigma_label: ", sigma_label)
    print(f"T0        : {T0} K")
    print("Carrier:")
    print("  q={} e".format(mobility.charge))
    print("  meff={}me".format(dos.meeff))
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    if mobility.charge < 0.0:
        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    else:
        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
    print("  mean free path prefactor: l0={} m".format(mobility.l0))

    if '***' in S_label:
        app.terminate(f"***Error: Choose S_label", pause = True)
    if '***' in n_label:
        app.terminate(f"***Error: Choose n_label", pause = True)
    if '***' in mu_label:
        app.terminate(f"***Error: Choose mu_label", pause = True)

    if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
        print("")
        print(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]")
        app.terminate(pause = True)

    print("")
    print("Read S data from {}".format(infile))
    datafile = tkVariousData(infile)
    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
    label_S, yS           = datafile.FindDataArray(S_label, flag = 'i')
# Convert V/K to uV/K
    for i in range(len(yS)):
        yS[i] *= 1.0e6 
    label_sigma, ysigma   = datafile.FindDataArray(sigma_label, flag = 'i')
    label_N, yN           = datafile.FindDataArray(n_label, flag = 'i')
    label_mu, ymu         = datafile.FindDataArray(mu_label, flag = 'i')
    label_sample, xsample = labels[0], datalist[0]
    if T_label == '' or '***' in T_label:
        label_T = 'T(K)'
        yT = []
        for i in range(len(xsample)):
            yT.append(T0)
    else:
        label_T, yT           = datafile.FindDataArray(T_label, flag = 'i')

#    print("x=", xsample)
    ndata = len(yS)
    print("ndata(all): ", ndata)
#    print("xsample=", label_sample, xsample)
#    print("yS=", label_S, yS)
#    print("ysigma=", label_sigma, ysigma)
    print(f"{'sample':15}  {'n(cm-3)':12}  {'mu(cm2/Vs)':10}  {'sigma(S/cm)':10}  {'S(uV/K)':10}")
    for i in range(ndata):
        sigma = ysigma[i]
        n = yN[i]
        print(f"{xsample[i]:15}  {n:12.4g}  {ymu[i]:10.4g}  {sigma:10.4g}  {yS[i]:10.4g}")


    print("")
    print("Weighted mobility")
    ysigma_w = []
    yn_w     = []
    ymu_w    = []
    ymu_w0   = []
    yEF      = []
    if yN[0] is not None:
        print("{:>10}\t{:>8}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>15}\t{:>15}"
            .format("sample", "T(K)", "EF(eV)", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)", "sigma_w", "mu_w(cm2/Vs)", "mu_w*(me/m*)^(3/2)", "n_w(cm-3)"))
    else:
        print("{:>10}\t{:>10}\t{:>10}\t{:>10}".format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)"))
    for i in range(ndata):
        sample = xsample[i]
        S      = yS[i]
        sigma  = ysigma[i]
        n      = yN[i]
        mu     = ymu[i]
        T      = yT[i]
        print(f"sample={sample} S={S} uV/K sigma={sigma} S/cm")
        mu_w, n_w, sign, carriertype = weighted_mobility(sigma / 0.01, S * 1.0e-6, T0)
        mu_w   *= 1.0 * 1.0e4                       # cm2/Vs
        mu_w0  = mu_w * pow(mobility.meff, -1.5)
        n_w    *= 1.0e-6                            # cm^-3
        sigma_w = e * n_w * mu_w

        if n is not None:
            EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
        else:
            EF, diffEF, flag = dos.EF_from_electronSeebeck(S, T, mobility, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
        if not flag:
            app.terminate("Error in sim(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
                                usage = usage, pause = True)
#        print("n=", n, EF)       

        ysigma_w.append(sigma_w)
        yn_w.append(n_w)
        ymu_w.append(mu_w)
        ymu_w0.append(mu_w0)
        yEF.append(EF)
        if n is not None:
            print("{:10.4g}\t{:8.3g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:15.4g}\t{:15.4g}"
                .format(sample, T, EF, S, sigma, mu, n, sigma_w, mu_w, mu_w0, n_w))
        else:
            print("{:10.4g}\t{:8.3g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, T, EF, S, sigma))

    Nmin_data = min(yN)
    Nmax_data = max(yN)
    if Nmin_data < Nmin:
        Nmin = Nmin_data
    if Nmax < Nmax_data:
        Nmax = Nmax_data

    lnNmin = log(Nmin)
    lnNmax = log(Nmax)
    lnNstep = (lnNmax - lnNmin) / (nN - 1)
    print("")
    print("N range: {:10.4g} - {:10.4g} cm-3, {:10.4g} step in ln(N), nN={}".format(Nmin, Nmax, lnNstep, nN))
    xNsim        = []
    xsigmasim    = []
    yScal        = []
    ySsim_nondeg = []
    ySsim_deg    = []
    print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
            .format("N(cm-3)", "EF(eV)", "S,cal(uV/K)", "S,non-deg(uV/K)", "S,deg(uV/K)", "sigma,cal(S/cm)", "Ne,cal(cm-3)", "mu,cal(cm2/Vs)"))
    for iN in range(nN):
        lnN = lnNmin + lnNstep * iN
        N   = exp(lnN)
#        print("N=", N, dos.NC)
        if N < dos.NC:
            EF = kB * T0 * log(N / dos.NC) / e
        else:
            EF = BMShift_FEA(N, dos.DC0)

        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
            = dos.cal_transport_properteis(T0, EF, mobility, klatt, validate_error_str = 'Error in sim()')

        Sndeg = dos.cal_S_nondegenerated_from_Ne(n, mobility.rfac, charge = mobility.charge) * 1.0e6    # uV/K
        Sdeg  = dos.cal_S_degenerated_from_Ne(T0, n, mobility.rfac, charge = mobility.charge) * 1.0e6  # uV/K

        xNsim.append(n)
        xsigmasim.append(sigma)
        yScal.append(S)
        ySsim_nondeg.append(Sndeg)
        ySsim_deg.append(Sdeg)
        print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
                .format(N, EF, S, Sndeg, Sdeg, sigma, n, mu))

    print("")
    print("Save parameters to [{}]".format(parameterfile))
    save_parameterfile(S2 = '')

    print("")
    print("Save data to [{}]".format(outxlsfile))
    if yN[0] is not None:
        tkVariousData().to_excel(outxlsfile, ["sample", 'S(uV/K)', 'sigma(S/cm)', 'mu(cm2/Vs)', 'N(cm-3)', "sigma_w", "mu_w", "mu_w*(me/m*)^(3/2)", "n_w"],
                         [xsample, yS, ysigma, ymu, yN, ysigma_w, ymu_w, ymu_w0, yn_w])
    else:
        tkVariousData().to_excel(outxlsfile, ["sample", 'S(uV/K)', 'sigma(S/cm)'], [xsample, yS, ysigma])

#=============================
# グラフの表示
#=============================
    print("")
    fig = plt.figure(figsize = figsize_sim)

    ax1  = fig.add_subplot(2, 2, 1)
    ax1b = ax1.twinx()
    ax2  = fig.add_subplot(2, 2, 2)
    ax4  = fig.add_subplot(2, 2, 3)
    ax4b = ax4.twinx()
    ax3  = fig.add_subplot(2, 2, 4)

    maxS = max(yS)
    minS = min(yS)
    if maxS > 0.0:
        maxS *= 2.0
    else:
        maxS = 0.0
    if minS > 0.0:
        minS = 0.0
    else:
        minS *= 2.0

    ins1 = ax1.plot( xsample, yS,       label = 'S',          linestyle = '-',      linewidth = 1.0, color = 'red',   marker = 'o', markersize = 5.0)
    ins2 = ax1b.plot(xsample, ysigma,   label = r'$\sigma$',   linestyle = '-',      linewidth = 1.0, color = 'blue',  marker = '^', markersize = 10.0)
    ins3 = ax1b.plot(xsample, ysigma_w, label = r'$\sigma_w$', linestyle = 'dashed', linewidth = 0.5, color = 'green', marker = '+', markersize = 15.0, markerfacecolor = 'green')
    ax1.set_xlabel("sample", fontsize = fontsize)
    ax1.set_ylabel(r"S (K/$\mu$V)", fontsize = fontsize)
    ax1b.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
#    ax1.set_ylim([minS, maxS])
#    h1a, l1a = ax1.get_legend_handles_labels()
#    h1b, l1b = ax1b.get_legend_handles_labels()
#    ax1.legend(h1a + h1b, l1a + h1b, fontsize = legend_fontsize, loc = 0)
    ins = ins1 + ins2 + ins3
    ax1.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
#    ax1.legend(loc = 0)
#    ax1b.legend(loc = 0)
    ax1.tick_params(labelsize = fontsize)
    ax1b.tick_params(labelsize = fontsize)

    if yN[0] is not None:
        ins1 = ax4.plot( xsample, ymu,    label = r'$\mu$',     linestyle = '-',      color = 'red',  linewidth = 0.5, marker = 'o', markersize = 10.0)
        ins2 = ax4.plot( xsample, ymu_w,  label = r'$\mu_w$',   linestyle = 'dashed', color = 'red',  linewidth = 0.5, marker = '^', markersize = 5.0)
        ins3 = ax4.plot( xsample, ymu_w0, label = r'$\mu_w^0$', linestyle = 'dashed', color = 'red',  linewidth = 0.5, marker = '+', markersize = 5.0)
        ins4 = ax4b.plot(xsample, yN,     label = '$N_e$',     linestyle = '-',      color = 'blue', linewidth = 0.5, marker = 'o', markersize = 10.0)
        ins5 = ax4b.plot(xsample, yn_w,   label = '$N_{e,w}$', linestyle = 'dashed', color = 'blue', linewidth = 0.5, marker = '^', markersize = 5.0)
        ax4.set_xlabel("sample", fontsize = fontsize)
        ax4.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
        ax4b.set_ylabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
#    ax4.set_ylim([minS, maxS])
        ax4b.set_yscale('log')
        ins = ins1 + ins2 + ins3 + ins4 + ins5  
        ax4.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
        ax4.tick_params(labelsize = fontsize)
        ax4b.tick_params(labelsize = fontsize)

    ax2.plot(ysigma,    yS,           label = 'Jonker plot',  linestyle = '',       color = 'red',   marker = 'o', markersize = 5.0)
    ax2.plot(xsigmasim, yScal,        label = 'S,cal',        linestyle = '-',      color = 'red',   linewidth = 1.0)
    ax2.plot(xsigmasim, ySsim_nondeg, label = 'S,sim,nondeg', linestyle = 'dashed', color = 'blue',  linewidth = 0.5)
    ax2.plot(xsigmasim, ySsim_deg,    label = 'S,sim,deg',    linestyle = 'dashed', color = 'green', linewidth = 0.5)
    ax2.set_xlabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
    ax2.set_ylabel(r"S ($\mu$V/K)", fontsize = fontsize)
    ax2.set_xscale('log')
    ax2.set_ylim([minS, maxS])
    ax2.legend(fontsize = legend_fontsize, loc = 'best')
    ax2.tick_params(labelsize = fontsize)

    if yN[0] is not None:
        ax3.plot(yN,    yS,           label = 'Pisarenko plot', linestyle = '', marker = 'o', markersize = 5.0)
    ax3.plot(xNsim, yScal,        label = 'S,cal',          linestyle = '-', color = 'red', linewidth = 1.0)
    ax3.plot(xNsim, ySsim_nondeg, label = 'S,sim,nondeg',   linestyle = 'dashed', color = 'blue', linewidth = 0.5)
    ax3.plot(xNsim, ySsim_deg,    label = 'S,sim,deg',      linestyle = 'dashed', color = 'green', linewidth = 0.5)
    ax3.set_xlabel(r"$N$ (cm$^{-3}$)", fontsize = fontsize)
    ax3.set_ylabel(r"S ($\mu$V/K)", fontsize = fontsize)
    ax3.set_xscale('log')
    ax3.set_ylim([minS, maxS])
    ax3.legend(fontsize = legend_fontsize, loc = 'best')
    ax3.tick_params(labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("", pause = True)

def init():
    print("")
    print("Carrier:")
    print("  q={} e".format(mobility.charge))
    print("  meff={}me".format(dos.meeff))
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
    print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
    print("  mean free path prefactor: l0={} m".format(mobility.l0))

    if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
        print("")
        print(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]")
        app.terminate(pause = True)

    print("")
    print("Read S data from {}".format(infile))
    label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
            = read_datafile(infile, usage = usage)
# Convert V/K to uV/K
    for i in range(len(yS)):
        yS[i] *= 1.0e6 

    ndata = len(xsample)
    print("ndata: ", ndata)

    if yN[0] is not None:
        print("{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}"
            .format("sample", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)"))
    else:
        print("{:>10}\t{:>10}\t{:>10}".format("sample", "S(uV/K)", "sigma(S/cm)"))
    for i in range(ndata):
        sample = xsample[i]
        S      = yS[i]
        sigma  = ysigma[i]
        n      = yN[i]
        mu     = ymu[i]
        if n is not None:
            print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma, mu, n))
        else:
            print("{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma))


    if yS[0] > 0.0:
        mobility.charge = 1.0
    elif yS[0] < 0.0:
        mobility.charge = -1.0
    for i in range(1, ndata):
        if yS[i] * mobility.charge < 0.0:
            app.terminate("Error in init(): S data include both positive and negative values: S[{}}={}, S[{}]={}".format(0, yS[0], i, yS[i]),
                                usage = usage, pause = True)

    print("")
    print("Save parameters to [{}]".format(parameterfile))
    save_parameterfile(S2 = '')

def calL():
    print("")
    print("Carrier:")
    print("  q={} e".format(mobility.charge))
    print("  meff={}me".format(dos.meeff))
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    if mobility.charge < 0.0:
        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    else:
        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
    print("  mean free path prefactor: l0={} m".format(mobility.l0))

    print("")
    print("Read T, Ne, and sigma data from {}".format(infile))
    datafile = tkVariousData(infile)
    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
#    print("labels: ", labels)
#    ncol  = len(datalist)
#    print("ncol: ", ncol)
    label_T, xT           = datafile.FindDataArray(r'^T[\s|\(|\[|$]', flag = 'i')
    label_sigma, ysigma   = datafile.FindDataArray(r'^sigma[$|\S]?.*$', flag = 'i')
    label_N, yN           = datafile.FindDataArray(r'^N[$|\S]?.*$', flag = 'i')
    ndata = len(xT)
    print("ndata: ", ndata)
#    print("xT=", label_T, xT)
#    print("ysigma=", label_sigma, ysigma)
#    print("yN=", label_N, yN)

    print("Calculat EF, L and kappa,e")
    yEF     = []
    yL      = []
    ykappa  = []
    print("kappa,e=sigma,obs*L,cal*T")
    print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
        .format("T(K)", "EF(eV)", "sigma(S/cm)", "Ne(cm^-3)", "L(Wohm/K^2)", "kappa,e"))
    for i in range(len(xT)):
        T = xT[i]
        n = yN[i]
        s = ysigma[i]

        if n is not None:
            EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
        if not flag:
            app.terminate("Error in calL(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
                            usage = usage, pause = True)
#        print("n=", n, EF)       

        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
            = dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in properteis()')

        kappa = L * s / 0.01 * T

        yEF.append(EF)
        yL.append(L)
        ykappa.append(kappa)

        print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
            .format(xT[i], EF, ysigma[i], yN[i], L, kappa))

    header, ext     = os.path.splitext(infile)
    filebody        = os.path.basename(header)
    outxlsLfile     = f'{filebody}-L-me{dos.meeff}-r{mobility.rfac}.xlsx'
    print("")
    print("Save data to [{}]".format(outxlsLfile))
    tkVariousData().to_excel(outxlsLfile, ["T(K)", "EF,cal(eV)", "sigma,obs(S/cm)", "Ne,obs(cm^-3)", "L,cal(Wohm/K^2)", "kappa,e,cal=sigma,obs*L,cal*T(W/m/K)"], 
                          [xT, yEF, ysigma, yN, yL, ykappa])

#=============================
# グラフの表示
#=============================
    print("")
    fig = plt.figure(figsize = figsize_sim)

    if min(xT) == max(xT):
        xX = yN
        xlabel = '$N$ (cm$^{-3}$)'
        xscale = 'log'
    else:
        xX = xT
        xlabel = '$T$ (K)'
        xscale = 'linear'

    axEF    = fig.add_subplot(2, 2, 1)
    axNe    = fig.add_subplot(2, 2, 2)
    axsigma = axNe.twinx()
    axL     = fig.add_subplot(2, 2, 3)
    axkappa = axL.twinx()
    axkappa = fig.add_subplot(2, 2, 4)

    axEF.plot(xX, yEF, label = 'EF(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
    axEF.set_xlabel(xlabel, fontsize = fontsize)
    axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize)
    axEF.set_xscale(xscale)
    axEF.legend(fontsize = legend_fontsize, loc = 'best')
    axEF.tick_params(labelsize = fontsize)

    ins1  = axNe.plot   (xX, yN,     label = '$N_e$(obs)',           linestyle = '-', linewidth = 1.0, color = 'red',    marker = 'o', markersize = 5.0)
    ins2  = axsigma.plot(xX, ysigma, label = r'$\sigma$(obs)',        linestyle = 'dashed', linewidth = 1.0, color = 'blue', marker = 's', markersize = 5.0)
    axNe.set_xlabel(xlabel, fontsize = fontsize)
    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axNe.set_xscale(xscale)
    axNe.set_yscale('log')
    axsigma.set_ylabel("$\sigma$ (S/cm)", fontsize = fontsize)
    axsigma.set_yscale('log')
    ins = ins1 + ins2
    axNe.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
    axNe.tick_params(labelsize = fontsize)
    axsigma.tick_params(labelsize = fontsize)

    axL.plot(xX, yL, label = 'L(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
    axL.set_xlabel(xlabel, fontsize = fontsize)
    axL.set_ylabel(r"$L$ (W$\Omega$/K$^2$)", fontsize = fontsize)
    axL.set_xscale(xscale)
    axL.legend(fontsize = legend_fontsize, loc = 'best')
    axL.tick_params(labelsize = fontsize)
    
    axkappa.plot(xX, ykappa, label = r'$\kappa_e$(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
    axkappa.set_xlabel(xlabel, fontsize = fontsize)
    axkappa.set_ylabel(r"$\kappa_e$=$\sigma$$L$$T$ (W/m/K)", fontsize = fontsize)
    axkappa.set_xscale(xscale)
    axkappa.set_yscale('log')
    axkappa.legend(fontsize = legend_fontsize, loc = 'best')
    axkappa.tick_params(labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("", pause = True)

def calF():
    print("")
    print("infile     :", infile)
    print("T_label    : ", T_label)
    print("n_label    : ", n_label)
    print("mu_label   : ", mu_label)
    print("sigma_label: ", sigma_label)
    print("Carrier:")
    print("  q={} e".format(mobility.charge))
    print("  meff={}me".format(dos.meeff))
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    if mobility.charge < 0.0:
        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    else:
        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
    print("  mean free path prefactor: l0={} m".format(mobility.l0))

    print("")
    print("Read T, Ne, and sigma data from {}".format(infile))
    datafile = tkVariousData(infile)
    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
    label_T, xT           = datafile.FindDataArray(T_label, flag = 'i')
    label_sigma, ysigma   = datafile.FindDataArray(sigma_label, flag = 'i')
    label_N, yN           = datafile.FindDataArray(n_label, flag = 'i')
    label_mu, ymu         = datafile.FindDataArray(mu_label, flag = 'i')

    ndata = len(xT)
    print("ndata: ", ndata)

    print("")
    print("Calculat EF, FH, mu_drift, and Ne")
    yEF      = []
    yFH      = []
    yRH      = []
    yndrift  = []
    ymudrift = []
    print("kappa,e=sigma,obs*L,cal*T")
    print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
        .format("T(K)", "EF(eV)", "sigma(S/cm)", "Ne,Hall(obs)", "mu,Hall(obs)", "Ne(cm-3)", "mu,drift(cm2/Vs)", "FH", "RH(cm3/C)"))
    for i in range(len(xT)):
        T = xT[i]
        n = yN[i]
        s = ysigma[i]

        if n is not None:
            EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
        if not flag:
            app.terminate("Error in calF(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter), 
                        usage = usage, pause = True)

#        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
#            = dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in properteis()')
        inf = dos.cal_Hall_properteis(T0, EF, mobility, validate_error_str = 'Error in properteis()')
        FH       = inf["FH"]
        RH       = inf["RH"]
        RH0      = inf["RH0"]
        n_drift  = inf["nHall"]
        mu_drift = inf["muHall"]

        yEF.append(EF)
        yFH.append(FH)
        yRH.append(RH)
        yndrift.append(yN[i] * FH)
        ymudrift.append(ymu[i] / FH)

        print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
            .format(xT[i], EF, ysigma[i], yN[i], ymu[i], yndrift[i], ymudrift[i], FH, RH))

    header, ext     = os.path.splitext(infile)
    filebody        = os.path.basename(header)
    outxlsLfile     = f'{filebody}-FH-me{dos.meeff}-r{mobility.rfac}.xlsx'
    print("")
    print("Save data to [{}]".format(outxlsLfile))
    tkVariousData().to_excel(outxlsLfile, ["T(K)", "EF,cal(eV)", "sigma,obs(S/cm)", "Ne,Hall(obs)(cm-3)", "mu,Hall(obs)(cm2/Vs)", 
                           "Ne=Ne,Hall*FH(cm-3)", "mu,drift=mu,Hall/FH(cm2/Vs)", "FH", "RH(cm3/C)"], 
                          [xT, yEF, ysigma, yN, ymu, yndrift, ymudrift, yFH, yRH])

#=============================
# グラフの表示
#=============================
    print("")
    fig = plt.figure(figsize = figsize_sim)

    if min(xT) == max(xT):
        xX = yN
        xlabel = '$N$ (cm$^{-3}$)'
        xscale = 'log'
    else:
        xX = xT
        xlabel = '$T$ (K)'
        xscale = 'linear'

    axEF    = fig.add_subplot(2, 3, 1)
    axsigma = fig.add_subplot(2, 3, 2)
    axNe    = fig.add_subplot(2, 3, 3)
    axmu    = fig.add_subplot(2, 3, 4)
    axFH    = fig.add_subplot(2, 3, 5)
    axRH    = fig.add_subplot(2, 3, 6)

    axEF.plot(xX, yEF, label = 'EF(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
    axEF.set_xlabel(xlabel, fontsize = fontsize)
    axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize)
    axEF.set_xscale(xscale)
    axEF.legend(fontsize = legend_fontsize, loc = 'best')
    axEF.tick_params(labelsize = fontsize)

    axsigma.plot(xX, ysigma, label = r'$\sigma$(obs)',        linestyle = 'dashed', linewidth = 1.0, color = 'blue', marker = 's', markersize = 5.0)
    axsigma.set_xlabel(xlabel, fontsize = fontsize)
    axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
    axsigma.set_xscale(xscale)
    axsigma.set_yscale('log')
    axsigma.legend(fontsize = legend_fontsize, loc = 'best')
    axsigma.tick_params(labelsize = fontsize)

    axNe.plot   (xX, yN,      label = '$N_e$$_{,Hall}(obs)',   linestyle = '-', linewidth = 1.0, color = 'red',   marker = 'o', markersize = 5.0)
    axNe.plot   (xX, yndrift, label = '$N_e$$_{,drift}$(obs)', linestyle = '-', linewidth = 0.5, color = 'green', marker = '^', markersize = 3.0)
    axNe.set_xlabel(xlabel, fontsize = fontsize)
    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
    axNe.set_xscale(xscale)
    axNe.set_yscale('log')
    axNe.legend(fontsize = legend_fontsize, loc = 'best')
    axNe.tick_params(labelsize = fontsize)

    axmu.plot(xX, ymu,      label = r'$\mu_{Hall}$(obs)',  linestyle = '-', linewidth = 1.0,  color = 'red',   marker = 'o', markersize = 5.0)
    axmu.plot(xX, ymudrift, label = r'$\mu_{drift}$(obs)', linestyle = '-', linewidth = 0.5,  color = 'green', marker = '^', markersize = 3.0)
    axmu.set_xlabel(xlabel, fontsize = fontsize)
    axmu.set_ylabel(r"$\mu$ (cm$^2$/Vs)", fontsize = fontsize)
    axmu.set_xscale(xscale)
    axmu.legend(fontsize = legend_fontsize, loc = 'best')
    axmu.tick_params(labelsize = fontsize)

    axFH.plot(xX, yFH, label = '$F_{Hall}$', linestyle = '-', linewidth = 1.0,  color = 'blue', marker = 'o', markerfacecolor = 'blue', markersize = 3.0)
    axFH.set_xlabel(xlabel, fontsize = fontsize)
    axFH.set_ylabel("$F_{Hall}$", fontsize = fontsize)
    axFH.set_xscale(xscale)
    axFH.legend(fontsize = legend_fontsize, loc = 'best')
    axFH.tick_params(labelsize = fontsize)

    axRH.plot(xX, yRH, label = '$R_H$',      linestyle = '-', linewidth = 1.0,  color = 'blue',  marker = 'o', markerfacecolor = 'blue', markersize = 3.0)
    axRH.set_xlabel(xlabel, fontsize = fontsize)
    axRH.set_ylabel("$R_H$ (cm$^3$/C)", fontsize = fontsize)
    axRH.set_xscale(xscale)
    axRH.set_yscale('log')
    axRH.legend(fontsize = legend_fontsize, loc = 'best')
    axRH.tick_params(labelsize = fontsize)

    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("", pause = True)

def T():
    global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu

    print("")
    print("Carrier:")
    print("  q={} e".format(mobility.charge))
    print("  meff={}me".format(dos.meeff))
    dos.NC  = meff2NC_FEA(dos.meeff, T0)
    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
    if mobility.charge < 0.0:
        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    else:
        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
    print("  mean free path prefactor: l0={} m".format(mobility.l0))

    print("")
    print("Read S data from {}".format(infile))
    label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
            = read_datafile(infile, usage = usage)
    label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
            = construct_lists(label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu)
    print("x=", xsample)
    ndata = len(xsample)
    print("ndata: ", ndata)

    if yN[0] is not None:
        print("{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}"
            .format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)"))
    else:
        print("{:>10}\t{:>10}\t{:>10}\t{:>10}".format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)"))
    for i in range(ndata):
        sample = xsample[i]
        S      = yS[i]
        sigma  = ysigma[i]
        n      = yN[i]
        mu     = ymu[i]
        if n is not None:
            print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma, mu, n))
        else:
            print("{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma))

    Nmin_data = min(yN)
    Nmax_data = max(yN)
    if Nmin_data < Nmin:
        Nmin = Nmin_data
    if Nmax < Nmax_data:
        Nmax = Nmax_data
    print("")
    print(f"Plot N range: {Nmin:12.4g} - {Nmax:12.4g}")

    lnNmin = log(Nmin)
    lnNmax = log(Nmax)
    lnNstep = (lnNmax - lnNmin) / (nN - 1)
    print("")
    print("N range: {:10.4g} - {:10.4g} cm-3, {:10.4g} step in ln(N), nN={}".format(Nmin, Nmax, lnNstep, nN))
    xEFcal       = make_matrix2(nT, nN, 0.0)
    xNsim        = make_matrix2(nT, nN, 0.0)
    ysigmasim    = make_matrix2(nT, nN, 0.0)
    ymusim       = make_matrix2(nT, nN, 0.0)
    yScal        = make_matrix2(nT, nN, 0.0)
    ykappacal    = make_matrix2(nT, nN, 0.0)
    yPFcal       = make_matrix2(nT, nN, 0.0)
    yZTcal       = make_matrix2(nT, nN, 0.0)
    Tstep = (Tmax - Tmin) / (nT - 1)
    for iT in range(nT):
        T = Tmin + iT * Tstep
        dos.NC  = meff2NC_FEA(dos.meeff, T)
        dos.DC0 = meff2DC0_FEA(dos.meeff, T)

        print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
                .format("T(K)", "N(cm-3)", "EF(eV)", "S,cal(uV/K)", "sigma,cal(S/cm)", "Ne,cal(cm-3)", "mu,cal(cm2/Vs)", "kappa,e,cal(W/m/K)", "PF(uW/cm/K^2)", "ZT"))
        for iN in range(nN):
            lnN = lnNmin + lnNstep * iN
            N   = exp(lnN)
#           print("N=", N, dos.NC)
            if N < dos.NC:
                EF = kB * T0 * log(N / dos.NC) / e
            else:
                EF = BMShift_FEA(N, dos.DC0)

            EF, diffEF, flag = dos.EF_from_electrondensity(N, T, EF0 = EF, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
            if not flag:
               app.terminate("Error in sim(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
                                usage = usage, pause = True)

            sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
                = dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in sim()')

            xEFcal[iT][iN]     = EF
            xNsim[iT][iN]      = n
            ysigmasim[iT][iN]  = sigma
            ymusim[iT][iN]     = mu
            yScal[iT][iN]      = S
            ykappacal[iT][iN]  = kappa
            yPFcal[iT][iN]     = PF
            yZTcal[iT][iN]     = ZT
            print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
                    .format(T, N, EF, S, sigma, n, mu, kappa, PF, ZT))

    print("")
    print("Save parameters to [{}]".format(parameterfile))
    save_parameterfile(S2 = '')


    print("")
    labels = ["N(cm-3)"] + ["{:.4g} K".format(Tmin + iT * Tstep) for iT in range(nT)]
    header, ext     = os.path.splitext(infile)
    filebody        = os.path.basename(header)
    xlsfile      = filebody + '-EF.xlsx'.format(T)
    print("Save data to [{}]".format(xlsfile))
#    print("labels=", labels)
#    print("datalist=", [xNsim[0]])# + [xEFcal[iT] for iT in range(nT)])
    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + xEFcal) #[xEFcal[iT] for iT in range(nT)])
    xlsfile   = filebody + '-sigma.xlsx'.format(T)
    print("Save data to [{}]".format(xlsfile))
    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ysigmasim)
    xlsfile   = filebody + '-mu.xlsx'.format(T)
    print("Save data to [{}]".format(xlsfile))
    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ymusim)
    xlsfile   = filebody + '-S.xlsx'.format(T)
    print("Save data to [{}]".format(xlsfile))
    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yScal)
    xlsfile   = filebody + '-kappa.xlsx'.format(T)
    print("Save data to [{}]".format(xlsfile))
    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ykappacal)
    xlsfile   = filebody + '-PF.xlsx'.format(T)
    print("Save data to [{}]".format(xlsfile))
    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yPFcal)
    xlsfile   = filebody + '-ZT.xlsx'.format(T)
    print("Save data to [{}]".format(xlsfile))
    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yZTcal)

#=============================
# グラフの表示
#=============================
    print("")
    fig = plt.figure(figsize = figsize)

    axSsobs    = fig.add_subplot(3, 3, 1)
    axSsobsb   = axSsobs.twinx()
    axmuNobs   = fig.add_subplot(3, 3, 2)
    axmuNobsb  = axmuNobs.twinx()
    axNs       = fig.add_subplot(3, 3, 3)
    axNmu      = fig.add_subplot(3, 3, 4)
    axNS       = fig.add_subplot(3, 3, 5)
    axNkappa   = fig.add_subplot(3, 3, 6)
    axNPF      = fig.add_subplot(3, 3, 7)
    axNZT      = fig.add_subplot(3, 3, 8)
    axNEF      = fig.add_subplot(3, 3, 9)

    maxS = max(yS)
    minS = min(yS)
    if maxS > 0.0:
        maxS *= 1.2
    else:
        maxS = 0.0
    if minS > 0.0:
        minS = 0.0
    else:
        minS *= 1.2

    ins1 = axSsobs.plot (xsample, yS,       label = 'S',          linestyle = '-',      linewidth = 1.0, color = 'red',   marker = 'o', markersize = 5.0)
    ins2 = axSsobsb.plot(xsample, ysigma,   label = r'$\sigma$',   linestyle = '-',      linewidth = 1.0, color = 'blue',  marker = '^', markersize = 10.0)
    axSsobs.set_xlabel("sample", fontsize = fontsize)
    axSsobs.set_ylabel(r"S (K/$\mu$V)", fontsize = fontsize)
    axSsobsb.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
    ins = ins1 + ins2 
    axSsobs.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
    axSsobs.tick_params(labelsize = fontsize)
    axSsobsb.tick_params(labelsize = fontsize)

    if yN[0] is not None:
        ins1 = axmuNobs.plot( xsample, ymu,    label = r'$\mu$',     linestyle = '-',      color = 'red',  linewidth = 0.5, marker = 'o', markersize = 10.0)
        ins2 = axmuNobsb.plot(xsample, yN,     label = '$N_e$',     linestyle = '-',      color = 'blue', linewidth = 0.5, marker = 'o', markersize = 10.0)
        axmuNobs.set_xlabel("sample", fontsize = fontsize)
        axmuNobs.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
        axmuNobsb.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
        axmuNobsb.set_yscale('log')
        ins = ins1 + ins2
        axmuNobs.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
        axmuNobs.tick_params(labelsize = fontsize)
        axmuNobsb.tick_params(labelsize = fontsize)

    if yN[0] is not None:
        axNs.plot( yN, ysigma, label = r'$\sigma_{,obs}$', linestyle = '', marker = 'o', markersize = 5.0)
        axNmu.plot(yN, ymu,    label = '$N_e$$_{,obs}$', linestyle = '', marker = 'o', markersize = 5.0)
        axNS.plot( yN, yS,     label = '$S_{obs}$', linestyle = '', marker = 'o', markersize = 5.0)
    for iT in range(nT):
        T= Tmin + iT * Tstep
        axNs.plot(    xNsim[iT], ysigmasim[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
        axNmu.plot(   xNsim[iT], ymusim[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
        axNS.plot(    xNsim[iT], yScal[iT],     label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
        axNkappa.plot(xNsim[iT], ykappacal[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
        axNPF.plot(   xNsim[iT], yPFcal[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
        axNZT.plot(   xNsim[iT], yZTcal[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
        axNEF.plot(   xNsim[iT], xEFcal[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
    axNs.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
    axNmu.set_xlabel(r"$N$ (cm$^{-3}$)", fontsize = fontsize)
    axNS.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
    axNkappa.set_xlabel( "$N$ (cm$^{-3}$)", fontsize = fontsize)
    axNPF.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
    axNZT.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
    axNEF.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
    axNs.set_ylabel( r"$\sigma$ (S/cm)", fontsize = fontsize)
    axNmu.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
    axNS.set_ylabel( "S (V/K)", fontsize = fontsize)
    axNkappa.set_ylabel( r"$\kappa$$_{,e}$ (W/m/K)", fontsize = fontsize)
    axNPF.set_ylabel( r"PF ($\mu$W/cm/K$^2$)", fontsize = fontsize)
    axNZT.set_ylabel( "ZT", fontsize = fontsize)
    axNEF.set_ylabel( "$E_F$", fontsize = fontsize)
    axNs.set_xscale('log')
    axNmu.set_xscale('log')
    axNS.set_xscale('log')
    axNkappa.set_xscale('log')
    axNPF.set_xscale('log')
    axNZT.set_xscale('log')
    axNEF.set_xscale('log')
    axNs.set_yscale('log')
    axNkappa.set_yscale('log')
    axNs.legend(fontsize = legend_fontsize, loc = 'best')
    axNmu.legend(fontsize = legend_fontsize, loc = 'best')
    axNS.legend(fontsize = legend_fontsize, loc = 'best')
    axNkappa.legend(fontsize = legend_fontsize, loc = 'best')
    axNPF.legend(fontsize = legend_fontsize, loc = 'best')
    axNZT.legend(fontsize = legend_fontsize, loc = 'best')
    axNEF.legend(fontsize = legend_fontsize, loc = 'best')
    axNs.tick_params(labelsize = fontsize)
    axNmu.tick_params(labelsize = fontsize)
    axNS.tick_params(labelsize = fontsize)
    axNkappa.tick_params(labelsize = fontsize)
    axNPF.tick_params(labelsize = fontsize)
    axNZT.tick_params(labelsize = fontsize)
    axNEF.tick_params(labelsize = fontsize)
    
    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("", pause = True)


def main():
    global app
    global parameterfile, parameterbkfile, outxlsfile, outxlsPropfile
    global outxlsfFDfile, outxlsFjfile
    
    app = tkApplication()
    logfile = app.replace_path(infile)
    print(f"Open logfile [{logfile}]")
    app.redirect(targets = ["stdout", logfile], mode = 'w')

    updatevars()

    parameterfile   = app.replace_path(infile, template = ["{dirname}", "{filebody}.in"])
    parameterbkfile = app.replace_path(infile, template = ["{dirname}", "{filebody}.in.bak"])
    outxlsfile      = app.replace_path(infile, template = ["{dirname}", "{filebody}-Seebeck-{mode}.xlsx"], ext_dict = {"mode": mode})
    outxlsPropfile  = app.replace_path(infile, template = ["{dirname}", "{filebody}-properties-{mode}.xlsx"], ext_dict = {"mode": mode})

    outxlsfFDfile   = app.replace_path(infile, template = ["{dirname}", "fFD.xlsx"])
    outxlsFjfile    = app.replace_path(infile, template = ["{dirname}", "Fj.xlsx"])

    print("")
    print("mode: ", mode)
    print("infile               : {}".format(infile))
    print("out xlsx file        : {}".format(outxlsfile))
    print("parameter file       : {}".format(parameterfile))
    print("parameter backup file: {}".format(parameterbkfile))
    print("")

    print("")
    print("Read [{}]".format(parameterfile))
    read_parameters(parameterfile)
# check args again to use given parameters
    updatevars()

    print("")
    print("Electronic structure")
#    print(f"Eg={dos.EC - dos.EV} eV")
#    print(f"EV={dos.EV} eV")
    print(f"EC={dos.EC} eV")
#    print(f"EA={dos.EA} eV")
#    print(f"NA={dos.NA} cm^-3")
#    print(f"ED={dos.ED} eV")
#    print(f"ND={dos.ND} cm^-3")
#    print(f"mh_eff={dos.mheff} me_0")
    print(f"me_eff={dos.meeff} me_0")
    print(f"Initial parameter")
    print(f"  EF0={dos.EF0} eV")
    print("Transport parameters")
    print(f"charge={mobility.charge} e")
    print(f"meff={mobility.meff}  me_0")
    print(f"rfac={mobility.rfac}")
    print(f"l0  ={mobility.l0} m")

    if mode == 'basic':
        basic()
    elif mode == 'prop':
        properties()
    elif mode == 'Hall':
        Hall()
    elif mode == 'init':
        init()
    elif mode == 'sim':
        sim()
    elif mode == 'fit':
        fit()
    elif mode == 'T':
        T()
    elif mode == 'calL':
        calL()
    elif mode == 'calF':
        calF()
    else:
        app.terminate("Error in main: Invalide mode [{}]".format(mode), usage = usage, pause = True)

if __name__ == "__main__":
    main()