import os
import sys
from math import exp, sqrt, log
import numpy as np
from scipy import integrate             # 数値積分関数 integrateを読み込む
from scipy.interpolate import interp1d
from scipy.optimize import minimize
from matplotlib import pyplot as plt


from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg, analyze_varstr
from tklib.tkutils import is_exist, is_file, is_dir
from tklib.tkutils import terminate, getarg, getintarg, getfloatarg, safe_getelement
from tklib.tkapplication import tkApplication

from tklib.tkparams import tkParams
from tklib.tkvariousdata import tkVariousData
from tklib.tkinifile import tkIniFile

from tklib.tkexcel import tkExcel

from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me, NA, R
from tklib.tkcrystal.tkcif import tkCIF
from tklib.tksci.tkFit import tkFit
from tklib.tkgraphic.tkplotevent import tkPlotEvent




"""
Simulate and fit Cp-T data to determine Debye temperature
"""


#=======================================
# プログラム開始
#=======================================
def initialize():
    app = tkApplication()
    cparams = tkParams()
    app.cparams = cparams

    cparams.Debug = 0

#mode: 'sim', 'fit'
    cparams.mode = 'sim'

    cparams.infile  = None #'Ba3SiO-heat capacity.xlsx'
    cparams.ciffile = None #'Ba3SiO-11.cif'

    cparams.Tlabel  = 0
    cparams.Cplabel = 1

    cparams.outxlsfile      = None
    cparams.outxlsTDfile    = None
    cparams.parameterfile   = None
    cparams.parameterbkfile = None

    cparams.xT     = None
    cparams.yCp    = None
    cparams.yCpini = None
    cparams.Cpmax  = None

# Debye Temperature, Cp at high T
    cparams.TD  = 300.0 # K
    cparams.Cp0 = 0.0

# Simulation T range
    cparams.Tmin  =    0.0
    cparams.Tmax  = 2000.0
    cparams.nT    = 2001

# Fitting T range
    cparams.Tfitmin = 0.0
    cparams.Tfitmax = 0.0
    cparams.Tfit0   = None
    cparams.Tfit1   = None
    cparams.Tcalmin = '*'
    cparams.Tcalmax = '*'

# 最適化パラメータの初期値
    cparams.varname = [      "TD",       "Cp0"]
    cparams.varunit = [       "K",     "J/K/g"]
    cparams.ai0     = [cparams.TD, cparams.Cp0]
    cparams.optid   = [         1,           1]

# Debye function
    cparams.interp_type = 'linear'
#    cparams.interp_type = 'cubic'
    cparams.TD0     = 300.0
    cparams.TfDmin  = 0.0
    cparams.TfDmax  = 2000.0
    cparams.TfDstep = 1.0

# TD determination by secant method
    cparams.h_TD       = 1.0e-5
    cparams.tol_TD     = 1.0e-3
    cparams.maxiter_TD = 100
    cparams.dump_TD    = 0.1

# scipy.optimize.minimizeで使うアルゴリズム
#nelder-mead Downhill simplex
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
#newton-cg Newton-CG
    cparams.method = "nelder-mead"
    cparams.maxiter = 1000
    cparams.tol    = 1.0e-4
    cparams.h_diff = 1.0e-3
    cparams.outputinterval = 1

#=============================
# Graph configuration
#=============================
    cparams.fig = None
    cparams.figsize = (10, 6)
    cparams.fontsize = 14
    cparams.legend_fontsize = 8
    cparams.graphupdateinterval = 5

    return app, cparams

def usage(app):
    cparams = app.cparams
    print("")
    print(" python {} mode infile ciffile args".format(sys.argv[0]))
    print(" (i)  python {} infile(.xlsx) ciffile TD optid(TD) Cp0 optid(Cp0) Tlabel Cplabel".format(sys.argv[0], 'sim'))
    print("  ex: python {} {} {} {} {} {}".format(sys.argv[0], 'sim', cparams.infile, cparams.ciffile, cparams.TD, cparams.optid[0], cparams.Cp0, cparams.optid[1]))
    print(" (ii) python {} infile(.xlsx) ciffile TD optid(TD) Cp0 optid(Cp0) Tlabel Cplabel Tfitmin Tfitmax Tcalmin Tcalmax ".format(sys.argv[0], 'fit'))
    print("  ex: python {} {} {} {} {} {}".format(sys.argv[0], 'fit', cparams.infile, cparams.ciffile, cparams.TD, cparams.optid[0], cparams.Cp0, cparams.optid[1]))

#=============================
# Treat argments
#=============================
def split_optstr(str):
    try:
        val, id = str.split(':')
        return pfloat(val), pint(id)
    except:
        return str, None

def parameter_list(app):
    return [app.cparams.TD, app.cparams.Cp0]

def set_parameters(app, ai):
    cparams = app.cparams
    cparams.TD, cparams.Cp0 = ai

def read_parameters(app, path, varname, ai, optid):
    cparams = app.cparams

    ini = tkIniFile()
    inf = ini.ReadAll(path, AddSection = 0)
    if inf is None:
        return

    for i in range(len(varname)):
        key = varname[i]
        str = inf.get(key, None)
        val, id = split_optstr(str)
        print("i=", i, key, val, optid[i])
        if val is not None:
            ai[i] = pfloat(val)
        if id is not None:
            optid[i] = id

    cparams.Tfitmin = pfloat(inf.get('Tfitmin', cparams.Tfitmin))
    cparams.Tfitmax = pfloat(inf.get('Tfitmax', cparams.Tfitmax))

    return inf

def save_parameters(path, varname, ai = None, optid = None, section = 'Parameters'):
    ini = tkIniFile(path)
    if type(varname) is list:
        for i in range(len(varname)):
            ini.WriteString(section, varname[i], "{}:{}".format(ai[i], optid[i]))
    else:
        for key in varname:
            ini.WriteString(section, key, "{}".format(varname[key]))
    
def print_parameters(varname, varunit, ai, optid):
    for i in range(len(varname)):
        print("  {:10s}: {:14.8g} optid={}".format(varname[i], ai[i], optid[i]))


def update_vars(app):
    cparams = app.cparams

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

    cparams.mode    = getarg( 1, cparams.mode)
    cparams.infile  = getarg( 2, cparams.infile)
    cparams.ciffile = getarg( 3, cparams.ciffile)
    cparams.TD       = getfloatarg( 4, cparams.TD)
    cparams.optid[0] = getintarg  ( 5, cparams.optid[0])
    cparams.Cp0      = getfloatarg( 6, cparams.Cp0)
    cparams.optid[1] = getintarg  ( 7, cparams.optid[1])
    cparams.Tlabel   = getarg     ( 8, cparams.Tlabel)
    cparams.Cplabel  = getarg     ( 9, cparams.Cplabel)
    cparams.Tfitmin  = getfloatarg(10, cparams.Tfitmin)
    cparams.Tfitmax  = getfloatarg(11, cparams.Tfitmax)
    cparams.Tcalmin  = getarg     (12, cparams.Tcalmin)
    cparams.Tcalmax  = getarg     (13, cparams.Tcalmax)
    cparams.method   = getarg     (14, cparams.method)
    cparams.maxiter  = getintarg  (15, cparams.maxiter)
    cparams.tol      = getfloatarg(16, cparams.tol)


# 被積分関数を定義
def func(x):
    if x > 100.0:
        return 0.0
    if x < 1.0e-20:
        return x * x * exp(x)

    expx = exp(x)
    expx1 = expx - 1.0
    return pow(x, 4) * expx / expx1 / expx1

def debye_function(T, TD):
    if T <= 0.0:
        return 0.0
    else:
        y = TD / T
        ret = integrate.quad(func, 0.0, y)
        integ = ret[0]
#        integ = integrate_function_by_simpson(func, 0.0, y, nIntegMesh)
        f = 3.0 / pow(y, 3) * integ
        return f

def Cp(T, TD, Cp0):
    return Cp0 * debye_function(T, TD)

def fit(app):
    cparams = app.cparams

    print("")
    print("Fitting to Cp-T data by nonlinear least-squares method")

    print("")
    print("mode    : {}".format(cparams.mode))
    print("infile  : {}".format(cparams.infile))
    print("cif file: {}".format(cparams.ciffile))
    print(f"Tlabel : {cparams.Tlabel}")
    print(f"Cplabel: {cparams.Cplabel}")
    print(f"T fit range: {cparams.Tfitmin} - {cparams.Tfitmax} K")
    print(f"method     : {cparams.method}")
    print(f"maxiter    : {cparams.maxiter}")
    print(f"tol        : {cparams.tol}")
    print(f"output interval      : {cparams.outputinterval}")
    print(f"graph update interval: {cparams.graphupdateinterval}")

    if '***' in cparams.method:
        app.terminate("\nError: Choose method\n", pause = True)

    if cparams.Tlabel == '' or '***' in cparams.Tlabel:
        app.terminate("\nError: Null Tlabel is not allowed\n", pause = True)
    if cparams.Cplabel == '' or '***' in cparams.Cplabel:
        app.terminate("\nError: Null Cplabel is not allowed\n", pause = True)

    fit = tkFit(tol = cparams.tol, nmaxiter = cparams.maxiter,
                print_interval = cparams.outputinterval, plot_interval = cparams.graphupdateinterval)

    print("")
    print(f"Read data from [{cparams.infile}]")
    fit.read_data(cparams.infile, xlabel = cparams.Tlabel, ylabel = cparams.Cplabel, 
                    xmin = cparams.Tfitmin, xmax = cparams.Tfitmax, usage = lambda: usage(app))

    cparams.Tcalmin = pfloat(cparams.Tcalmin, defval = min(fit.x))
    cparams.Tcalmax = pfloat(cparams.Tcalmax, defval = max(fit.x))

    print(f"T cal range: {cparams.Tcalmin} - {cparams.Tcalmax} K")
    nTcal = 201
    Tcalstep = (cparams.Tcalmax - cparams.Tcalmin) / (nTcal - 1)
    
    print(f"TD : {cparams.TD} K")
    print(f"Cp0: {cparams.Cp0}")
    
    print("")
    print(f"Read cif file from {cparams.ciffile}")
    cif = tkCIF()
    cifdata = cif.ReadCIF(cparams.ciffile, print_level = 0)
    cif.Close()
    if not cifdata:
        app.terminate(f"Error: Could not get cif data from infile [{cifdata}]", pause = True)
#    cifdata.Print("sample|cell|properties|atomtypes|asymsites")

    cry  = cifdata.GetCrystal()
    cry.PrintInf("sample|cell|atomtypes|asymsites")

    fu   = cry.FormulaUnit()
    UMW  = cry.UnitcellWeight()
    d    = cry.Density()
    ad   = cry.AtomDensity()
    molg = ad / NA / d          # Atomic density in mol/g
    CvDP     = 3.0 * R          # J/K/mol
    CvDPmolg = CvDP * molg      # J/K/g
    print("")
    print("Formura unit      : {:12.6g}".format(fu))
    print("Unit cell weight  : {:12.6g}".format(UMW))
    print("Density           : {:12.6g} g/cm^3".format(d))
    print("Atom density      : {:12.6g} cm^-3".format(ad))
    print("                    {:12.6g} mol/g".format(molg))
    print("Dulong-Petit limit: {:12.6g} J/K/mol".format(CvDP))
    print("                    {:12.6g} J/K/g".format(CvDPmolg))

    print("")
    print("Fitting configuration:")
    cparams.Tfit0 = cparams.Tfitmin
    if cparams.Tfitmax == 0.0:
        cparams.Tfit1 = max(xT)
    else:
        cparams.Tfit1 = cparams.Tfitmax
    print("Fitting T range: {:10.4g} - {:10.4g} K".format(cparams.Tfit0, cparams.Tfit1))
    print("Algorism:", cparams.method)
    print("# of max interation:", cparams.maxiter)
    print("eps:", cparams.tol)
    print("output interval:", cparams.outputinterval)

    if cparams.Cp0 == 0.0:
        cparams.Cp0 = CvDPmolg
#        cparams.Cp0 = 3.0 * R * molg

    fit.varname = [      "TD",       "Cp0"]
    fit.unit    = [       "K",          ""]
    fit.pk      = [cparams.TD,       cparams.Cp0]
    fit.optid   = [cparams.optid[0], cparams.optid[1]]
    fit.func    = lambda T, *pk: Cp(T, pk[0], pk[1]) #debye_function(T, pk[0])
    nvars = len(fit.pk)

    print("")
    print("Calculate initial Cp-T")
    yini = fit.cal_ylist(fit.pk)
    fini = fit.minimize_func(fit.pk)

#=============================
# グラフの表示
#=============================
    print("")
    print("plot")
    cparams.xT  = fit.x
    cparams.yCp = fit.y
    cparams.fig, axes = plt.subplots(1, 2, figsize = cparams.figsize)
#    axes = [axes]
#    axes = axes.flatten()

    fit.initial_plot(data_axis = axes[0], error_axis = axes[1], yini = yini, fmin = fini, 
                fig = cparams.fig, fontsize = cparams.fontsize)

#=============================
# Optimization
#=============================
    print("")
    print(f"Nonlinear least-squares fitting by [{cparams.method}]")
    print(f"  tol={cparams.tol}")

    ai = fit.extract_parameters(fit.pk, fit.optid)
    print("  ai0  =", fit.pk)
    print("  optid=", fit.optid)
    print("  ai   =", ai)

    print("")
    pfin, ffin, success, res = fit.minimize(cparams.method)
    if success:
        print("")
        print(f"Converged at iteration: {res.nit}")
    else:
         print(f"Function did not converge")

    print("Final parameters")
    for i in range(nvars):
         print(f"  {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}")
    print(f"    f={ffin:12.6g}")

    set_parameters(app, pfin)

    print("")
    print("Optimized at S2={:12.6g}:".format(ffin))
    print_parameters(cparams.varname, cparams.varunit, pfin, fit.optid)

    print("Dulong-Petit limit: {:12.6g} J/K/mol".format(CvDP))
    print("                    {:12.6g} J/K/g".format(CvDPmolg))

#========================================
# Final result
#========================================
    yfin = fit.cal_ylist(pfin)
    print("")
    ndata = len(fit.x)
    nskip = int(ndata / 100)
    print("Final result")
    print(f"{'i':4} {'T':10} {'Cp':12} {'Cp(ini)':12} {'Cp(fin)':12}")
    for i in range(0, ndata, nskip):
        print(f"{i:4} {fit.x[i]:10.4g} {fit.y[i]:12.4g} {yini[i]:12.4g} {yfin[i]:12.4g}")

    xTcal = [cparams.Tcalmin + i * Tcalstep for i in range(nTcal)]
    ycal = fit.cal_ylist(pfin, xTcal)

    print("")
    print(f"Save simulation data to [{cparams.outxlsfile}]")
    tkVariousData().to_excel(cparams.outxlsfile, 
                ["T(K)", "Cp,obs(J/K/g)", "Cp,ini(J/K/g)", "Cp,fin(J/K/g)", "T(K)", "Cp,cal(J/K/g)"], 
                [fit.x, fit.y, yini, yfin, xTcal, ycal])

    print("")
    print(f"Save final parameters to [{cparams.parameterfile}]")
    save_parameters(cparams.parameterfile, fit.varname, pfin, fit.optid, section = 'Parameters')
    save_parameters(cparams.parameterfile, {'Tfitmin': cparams.Tfitmin, 'Tfitmax': cparams.Tfitmax, 'fmin': ffin}, section = 'Configuration')

#=============================
# グラフの表示
#=============================
    print("")
    print("Plot optimized")
#    fit.finalize_plot(yfin, iter = res.nit, fmin = ffin)
    axes[0].plot(xTcal, ycal, label = 'fit', color = 'blue', linestyle = '-', linewidth = 0.5)
    axes[0].plot(axes[0].get_xlim(), [CvDPmolg, CvDPmolg], label = 'Dulong-Petit limit', color = 'green', linestyle = 'dashed', linewidth = 0.5)
    ylim = axes[0].get_ylim()
    axes[0].set_ylim([0.0, ylim[1] * 1.05])
    plt.pause(1.0e-5)

    app.terminate("", pause = True)


def sim(app):
    cparams = app.cparams
    
    print("")
    print("mode    : {}".format(cparams.mode))
    print("infile  : {}".format(cparams.infile))
    print("cif file: {}".format(cparams.ciffile))
    print("TD      : {} K".format(cparams.TD))
    print("Cp0     : {} J/K/g".format(cparams.Cp0))
    print(f"h_TD   : {cparams.h_TD} K")

    if cparams.Tlabel == '' or '***' in cparams.Tlabel:
        app.terminate("\nError: Null Tlabel is not allowed\n", pause = True)
    if cparams.Cplabel == '' or '***' in cparams.Cplabel:
        app.terminate("\nError: Null Cplabel is not allowed\n", pause = True)

    print("")
    print(f"Read cif file from {cparams.ciffile}")
    cif = tkCIF()
    cifdata = cif.ReadCIF(cparams.ciffile, print_level = 0)
    cif.Close()
    if not cifdata:
        app.terminate(f"Error: Could not get cif data from infile [{cparams.ciffile}]")
#    cifdata.Print()

    cry  = cifdata.GetCrystal()
    cry.PrintInf("sample|cell|atomtypes|asymsites")
    fu   = cry.FormulaUnit()
    UMW  = cry.UnitcellWeight()
    d    = cry.Density()
    ad   = cry.AtomDensity()
    molg = ad / NA / d          # Atomic density in mol/g
    CvDP     = 3.0 * R          # J/K/mol
    CvDPmolg = CvDP * molg      # J/K/g
    if cparams.Cp0 == 0.0:
        cparams.Cp0 = 3.0 * R * molg        # J/g/K
    print("")
    print("Formura unit      : {:12.6g}".format(fu))
    print("Unit cell weight  : {:12.6g}".format(UMW))
    print("Density           : {:12.6g} g/cm^3".format(d))
    print("Atom density      : {:12.6g} cm^-3".format(ad))
    print("                    {:12.6g} mol/g".format(molg))

    print("")
    print(f"Read Cp-T data from {cparams.infile}")
    datafile = tkVariousData(cparams.infile)
    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, force_numeric = False, usage = lambda: usage(app))
    _Tlabel, xT   = datafile.FindDataArray(cparams.Tlabel, flag = 'i')
    _Cplabel, yCp = datafile.FindDataArray(cparams.Cplabel, flag = 'i')

    Cpmax = CvDPmolg #max(yCp)
    Cpapproxlimit = 0.036 * Cpmax
    print("")
    print("Dulong-Petit limit: {:12.6g} J/K/mol".format(CvDP))
    print("                    {:12.6g} J/K/g".format(CvDPmolg))


    def cal_TD(T, _Cp, Cpmax, TD0, dump, h, maxiter, tol):
        TD = TD0
        d = _Cp - Cp(T, TD, Cpmax)
        for i in range(maxiter):
            TDp = TD + h
            dp = _Cp - Cp(T, TDp, Cpmax)
            dfdx = (dp - d) / h
            dx = -d / dfdx * dump
            if abs(dx) < tol:
                break

            TD += dx
            d = dp

        return TD, d, dx

    cparams.Tcalmin = pfloat(cparams.Tcalmin, defval = min(xT))
    cparams.Tcalmax = pfloat(cparams.Tcalmax, defval = max(xT))
    print("")
    print( "T dependence of Debye tempearature assuming the Debye model")
    print(f"Calculation T range: {cparams.Tcalmin} - {cparams.Tcalmax} K")
    print(f"  Cp(max): {Cpmax:10.4g} J/K/g")
    print( "  Cp lower limit for TD approximation: {:10.4g} J/K/g".format(Cpapproxlimit))
    print(f"  dump : {cparams.dump_TD}")
    print(f"  h    : {cparams.h_TD}")
    print(f"  maxiter: {cparams.maxiter_TD}")
    print(f"  tol    : {cparams.tol_TD}")
    
    yTdapprox = []
    yTdexact = []
#    print("{:10}\t{:10}\t{:10}".format("T(KT)", "Cp(obs)(J/K/g)", "TD(approx)(K)"))
#    print(f"{'T':10}\t{'Cp(obs)':10}\t{'TD(approx)':10}\t{'TD(exact)':10}\t{'d_Cp':8}\t{'dT':8}")
    for i in range(len(xT)):
        T  = xT[i]
        _Cp = yCp[i]
        Td = 49.0 + 56.0 * (Cpmax / _Cp - 1.0)
        Td = -35.0 + 5.0 * sqrt(Td)
        Td = sqrt(Td) * T
        yTdapprox.append(Td)
        if _Cp >= Cpmax:
            Cp_approx = None
            TD_exact = None
            d = None
            dx = None
            yTdexact.append(None)
        else:
            Cp_approx = Cp(T, Td, Cpmax)
            TD_exact, d, dx = cal_TD(T, _Cp, Cpmax, Td, cparams.dump_TD, cparams.h_TD, cparams.maxiter_TD, cparams.tol_TD)
            yTdexact.append(TD_exact)

#        if dx is not None and abs(dx) > cparams.tol_TD:
#            print(f"***Warning: Cp calculation from the approxmated TD={Td:10.4g} did not converge")
#        print(f"{T:10.4g}\t{yCp[i]:10.4g}\t{yTdapprox[i]:10.4g}\t{TD_exact:10.4g}\t{d:8.3g}\t{dx:8.3g}")

    print("")
    Tstep = (cparams.Tcalmax - cparams.Tcalmin) / (cparams.nT - 1)
    print(f"T range: {cparams.Tmin} - {cparams.Tmax} K, {Tstep} K step, nT={cparams.nT}")
    print(f"Calculate Cp-T data by Debye model with Cp0={cparams.Cp0:10.4g} J/K/g and TD={cparams.TD:10.4g} K")
    xTsim  = []
    yfD    = []
    yCpsim = []
#    print("  T(K)   T/TD        fD(T/TD)     Cp,sim(J/g/K)")
    for i in range(cparams.nT):
        Ti = cparams.Tcalmin + i * Tstep
# 温度が 0 K のときは y=>∞ となって計算できないので、極限値 0 を与える
        xTsim.append(Ti)
        fD = debye_function(Ti, cparams.TD)
        yfD.append(fD)
        yCpsim.append(yfD[i] * cparams.Cp0)         # J/g/K
#        print("%6.1f\t%6.4f\t%12.6g\t%8.4f" % (Ti, y, yfD[i], yCpsim[i]))

    print("Dulong-Petit limit: {:12.6g} J/K/mol".format(CvDP))
    print("                    {:12.6g} J/K/g".format(CvDPmolg))

    print("")
    print(f"Save observd Cp and TD data to [{cparams.outxlsTDfile}]")
#    save_xls(cparams.outxlsTDfile, ["T(K)", "Cp(J/K/g)", "TD,approx(K)", "TD,exact(K)"], [xT, yCp, yTdapprox, yTdexact])
    tkVariousData().to_excel(cparams.outxlsTDfile, 
                ["T(K)", "Cp(J/K/g)", "TD,approx(K)", "TD,exact(K)"], [xT, yCp, yTdapprox, yTdexact])

    print("")
    print(f"Save simulation data to [{cparams.outxlsfile}]")
    tkVariousData().to_excel(cparams.outxlsTDfile, 
                ["T(K)", "Debye function", "Cp,sim(J/K/g)"], [xTsim, yfD, yCpsim])

    print("")
    print(f"Save parameters to [{cparams.parameterfile}]")
    save_parameters(cparams.parameterfile, cparams.varname, [cparams.TD, cparams.Cp0], cparams.optid, section = 'Parameters')
    save_parameters(cparams.parameterfile, {'Tfitmin': cparams.Tfitmin, 'Tfitmax': cparams.Tfitmax, 
                                            'Cpmax': cparams.Cpmax}, section = 'Configuration')

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

    plot_event = tkPlotEvent(plt)

    ax1 = fig.add_subplot(1, 2, 1)
    ax2 = fig.add_subplot(1, 2, 2)
    plt.title(f"{cparams.infile} (TD={cparams.TD:8.3f} K)", fontsize = cparams.fontsize)

    input_data = ax1.plot(xT, yCp, label = '$C_p$(obs)', linestyle = '', 
                    marker = 'o', markeredgewidth = 0.0, markerfacecolor = 'blue', markersize = 2.0)
    simulaton_data = ax1.plot(xTsim, yCpsim, label = '$C_p$(sim)', linestyle = 'dashed', color = 'red', linewidth = 0.5)
    xlim = ax1.get_xlim()
    ax1.plot(xlim, [0.0, 0.0], color = 'red', linestyle = '-', linewidth = 1.0)
    DP_data = ax1.plot(xlim, [CvDPmolg, CvDPmolg], label = 'Dulong-Petit limit', color = 'green', linestyle = 'dashed', linewidth = 0.5)
    Cpappr_data = ax1.plot(xlim, [Cpapproxlimit, Cpapproxlimit], label = '$T_D$ approx lower limit', color = 'purple', linestyle = 'dashed', linewidth = 0.5)
    ax1.set_xlabel("T (K)")
    ax1.set_ylabel("$C_p$ (J/K/g)")
    ax1.set_xlim(xlim)
    ax1.legend(fontsize = cparams.legend_fontsize)
    ax1.tick_params(labelsize = cparams.fontsize)

    TDappr_data  = ax2.plot(xT, yTdapprox, label = '$T_D$(approx)', linewidth = 0.5)
    TDexact_data = ax2.plot(xT, yTdexact,  label = '$T_D$(exact)', linewidth = 0.5)
    ax2.set_xlabel("T (K)")
    ax2.set_ylabel("$T_D$ (K)")
    ax2.set_xlim(xlim)
    ax2.legend(fontsize = cparams.legend_fontsize)
    ax2.tick_params(labelsize = cparams.fontsize)

    plot_event.add_data({"label": "input",                  "plot_type": "2D", "axis": ax1, "data": input_data})
    plot_event.add_data({"label": "simulation",             "plot_type": "2D", "axis": ax1, "data": simulaton_data})
    plot_event.add_data({"label": "Dulong-Petit limit",     "plot_type": "2D", "axis": ax1, "data": DP_data})
    plot_event.add_data({"label": "Cp approximation limit", "plot_type": "2D", "axis": ax1, "data": Cpappr_data})
    plot_event.add_data({"label": "TD(approximation)",      "plot_type": "2D", "axis": ax2, "data": TDappr_data})
    plot_event.add_data({"label": "TD(exact)",              "plot_type": "2D", "axis": ax2, "data": TDexact_data})

    plot_event.register_event(fig, event = "button_press_event", 
                    callback = lambda event: plot_event.onclick(event))

    plt.tight_layout()
    plt.pause(0.1)

    app.terminate("\n", pause = True)

def main():
    print("")
    print("########################################################")
    print("# Simulate / fitting to Cp - T data by Debye model")
    print("########################################################")

    app, cparams = initialize()
    update_vars(app)

    cparams.logfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-out.txt"])

    print(f"Open logfile [{cparams.logfile}]")
    app.redirect(targets = ["stdout", cparams.logfile], mode = 'w')

    cparams.ai0 = parameter_list(app)
    cparams.parameterfile   = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in"])
    cparams.parameterbkfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in.bak"])
    cparams.outxlsfile      = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-{mode}.xlsx"], ext_dict = {"mode": cparams.mode})
    cparams.outxlsTDfile    = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-TD.xlsx"])

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

    print("")
    print(f"Read [{cparams.parameterfile}]")
    read_parameters(app, cparams.parameterfile, cparams.varname, cparams.ai0, cparams.optid)
    set_parameters(app, cparams.ai0)

# check args again to use given parameters
    update_vars(app)
    ai0 = parameter_list(app)

    print("")
    print("Initial parameters")
    print_parameters(cparams.varname, cparams.varunit, ai0, cparams.optid)
    
    if cparams.mode == 'sim':
        sim(app)
    elif cparams.mode == 'fit':
        fit(app)
    else:
        app.terminate("Error in main: Invalide mode [{}]".format(mode), usage = lambda: usage(app), pause = True)


if __name__ == "__main__":
    main()