import sys
import numpy as np
from math import log, exp
from matplotlib import pyplot as plt


"""
Regular solution model for a A - B two element system.
"""


#===================================
# physical constants
#===================================
pi = 3.14159265358979323846
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
e0 = 8.854418782e-12; # C²N^-1m^-2";
kB = 1.380658e-23 # JK^-1";
me = 9.1093897e-31 # kg";
R = 8.314462618 # J/K/mol

# eps to avoid log(0.0)
xeps = 1.0e-6

#===================================
# parameters
#===================================
HAA = 10.0 # kJ/mol
HBB = 5.0
Omg = 20.0

#===================================
# ranges
#===================================
xBmin = 0.0
xBmax = 1.0
xBstep = 0.05

Tmin = 600.0 # K
Tmax = 1500.0
Tstep = 100.0

#===================================
# figure configuration
#===================================
fontsize = 16
legend_fontsize = 12


#==============================================
# fundamental functions
#==============================================
# 実数値に変換できない文字列をfloat()で変換するとエラーになってプログラムが終了する
# この関数は、変換できなかったらNoneを返すが、プログラムは終了させない
def pfloat(str):
    try:
        return float(str)
    except:
        return None

# pfloat()のint版
def pint(str):
    try:
        return int(str)
    except:
        return None

# 起動時引数を取得するsys.argリスト変数は、範囲外のindexを渡すとエラーになってプログラムが終了する
# egtarg()では、範囲外のindexを渡したときは、defvalを返す
def getarg(position, defval = None):
    try:
        return sys.argv[position]
    except:
        return defval

# 起動時引数を実数に変換して返す
def getfloatarg(position, defval = None):
    return pfloat(getarg(position, defval))

# 起動時引数を整数値に変換して返す
def getintarg(position, defval = None):
    return pint(getarg(position, defval))

def usage():
    print("")
    print("Usage: Variables in () are optional")
    print(" python {} (HAA HBB Omega xBmin xBmax xBstep Tmin Tmax Tstep)".format(sys.argv[0]))
    print(" ex: python {} {} {} {} {} {} {} {} {} {}"
            .format(sys.argv[0], HAA, HBB, Omg, xBmin, xBmax, xBstep, Tmin, Tmax, Tstep))

def terminate():
    print("")
    usage()
    exit()

#==============================================
# update default values by startup arguments
#==============================================
argv = sys.argv
#if len(argv) == 1:
# terminate()

HAA = getfloatarg(1, HAA)
HBB = getfloatarg(2, HBB)
Omg = getfloatarg(3, Omg)
xBmin = getfloatarg(4, xBmin)
xBmax = getfloatarg(5, xBmax)
xBstep = getfloatarg(6, xBstep)
Tmin = getfloatarg(7, Tmin)
Tmax = getfloatarg(8, Tmax)
Tstep = getfloatarg(9, Tstep)


def main():
    global HAA, HBB, Omg
    global xBmin, xBmax, xBstep
    global Tmin, Tmax, Tstep

    nxB = int((xBmax - xBmin) / xBstep + 1.00001)
    nT = int((Tmax - Tmin) / Tstep + 1.00001)
    xBarray = [xBmin + ixB * xBstep for ixB in range(nxB)]

    print("")
    print("xB range: {} - {}, {} step, {} data".format(xBmin, xBmax, xBstep, nxB))
    print("T range : {} - {}, {} step, {} data".format(Tmin, Tmax, Tstep, nT))

    fig = plt.figure(figsize = (8, 8))
    ax1 = fig.add_subplot(1, 1, 1)
# ax2 = fig.add_subplot(2, 3, 2)

    dGmarray = np.empty([nT, nxB])
    print("")
    print("{:7} {:5} {:8} {:8} {:8}".format('T', 'xB', 'dHm', 'dSm*T', 'dGm'))
    for iT in range(nT):
        T = Tmin + iT * Tstep
        for ixB in range(nxB):
            xB = xBmin + ixB * xBstep
            xA = 1.0 - xB
            dHm = Omg * xA * xB

            if xA < xeps or xB < xeps:
                dSm = 0.0
            else:
                dSm = -1.0e-3 * R * (xA * log(xA) + xB * log(xB))

            dGm = dHm - dSm * T
            dGmarray[iT][ixB] = dGm
            print("{:7.1f} {:5.3f} {:8.4f} {:8.4f} {:8.4f}".format(T, xB, dHm, dSm * T, dGm))

        ax1.plot(xBarray, dGmarray[iT], label = '{} K'.format(T))

    print("")
    print("plot")
# ax1.set_ylim([-10, 10])
    ax1.set_xlabel("xB")
    ax1.set_ylabel("dG, dH, dS (kJ/mol)")
    ax1.legend(fontsize = legend_fontsize)
    ax1.tick_params(labelsize = fontsize)
# plt.tight_layout()

    plt.pause(0.1)
    print("Press ENTER to exit>>", end = '')
    input()

    terminate()


if __name__ == '__main__':
    main()