###
# parameters.py (C) of thesis research "Modeling and Simulating Three Phases of Steel: Austenite, Feritte and Cementite" 2018 by Thijs Verbeek
# is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
#
# You should have received a copy of the license along with this
# work. If not, see <http://creativecommons.org/licenses/by-nc-sa/4.0/>.
###

import numpy as np

def t_end():
    return np.inf

def T0():
    return 995

def T_end():
    return 800

def tT(cooling_rate):
    return (T0()-T_end())/cooling_rate

def T_func(t,cooling_rate):
    if cooling_rate > 0:
        #temperature decrease ends at t = tT
        if type(t) == np.float64 or type(t) == np.float or type(t) == np.int:
            if t < tT(cooling_rate):
                T = T0() - t * cooling_rate
            else:
                T = T_end()
        else:
            T = T0() - t * cooling_rate
            T[T>T_end()] = T_end()

    elif cooling_rate == 0:
        #T is constant
        T = T0()
    return T

### Values for equilibrium concentrations might have to be approximated better. 
### Didn't have enough data to approximate all with good exponential fit.
### It works, but more data might give better results for lower temperatures 
### (which at some point makes the gamma/alpha non-linear iteration not converge)

def wt_alphagamma(T):
#    wt_ag = .022/100
#    wt_ag = (-1.0e-4*T + 1.19e-1)/100
    #1185.15 K - 1000 K lin fit, constraint on T = 1000 w_at = 0.019/100
    wt_ag = (-(.019/185.15) * T + 1185.15 * (.019/185.15))/100
    return wt_ag

def wt_gammaalpha(T):
#    wt_ga = .76/100
#    wt_ga = (-4.54e-3*T + 5.258)/100
    #1185.15 K - 1000 K lin fit
    wt_ga = (-(.76/185.15) * T + 1185.15 * (.76/185.15))/100
    return wt_ga

def wt_alphatheta(T):
#    wt_at = .017/100
#    wt_at = (4.03428571e-5*T - 2.92495238e-2)/100
#    #4 datapoints
#    wt_at = (1.65982194*np.exp(-8.99149745e3*(1/T-4.92633220e-4))+6.87671081e-5)/100
    #200 datapoints arrhenius fit, constraint on T = 1000 w_at = 0.019/100
    c = -9613.28161688
    a = np.exp(-c/1000)*0.019
    wt_at = a*np.exp(c/T)/100
    return wt_at

def wt_gammatheta(T):
    wt_gt = (2.8e-3*T - 2.05)/100
    #1185.15 K - 1000 K lin fit, constraint on T = 1000 w_at = 0.76/100
    wt_gt = (0.00274545*T - 1.98545455)/100
    #6 datapoints arrhenius fit, constraint on T = 1000 w_at = 0.76/100
    c = -3261.06836816
    a = np.exp(-c/1000)*0.76
    wt_gt = a*np.exp(c/T)/100
    return wt_gt

def wt_theta():
    return 6.6666666667/100

def c0_theta():
    return 7730.14*10**3*(wt_theta()/(1-wt_theta())) #g/(micro m)^3 #7730.14*10**3*(wt_theta()/(1-wt_theta())) #g/m^3



def c_alphagamma(T):
    #m
    a_alpha = .29e-9 # m
#    #micro m
#    a_alpha = .29e-3 #micro m
    NA = 6.022140857e23 #1/mol
    M_Fe = 55.845#g/mol
    c_ag = 2*M_Fe/(a_alpha**3*NA)*(wt_alphagamma(T)/(1-wt_alphagamma(T)))#/(1e-18) #g/m^3
    return c_ag

def c_gammaalpha(T):
    #m
    a_gamma = .36e-9 # m
#    #micro m
#    a_gamma = .36e-3 #micro m
    NA = 6.022140857e23 #1/mol
    M_Fe = 55.845#g/mol
    c_ga = 4*M_Fe/(a_gamma**3*NA)*(wt_gammaalpha(T)/(1-wt_gammaalpha(T)))#/(1e-18) #g/m^3
    return c_ga

def c_alphatheta(T):
    #m
    a_alpha = .29e-9 # m
#    #micro m
#    a_alpha = .29e-3 #micro m
    NA = 6.022140857e23 #1/mol
    M_Fe = 55.845#g/mol
    c_at = 2*M_Fe/(a_alpha**3*NA)*(wt_alphatheta(T)/(1-wt_alphatheta(T)))#/(1e-18) #g/m^3
    return c_at

def c_gammatheta(T):
    #m
    a_gamma = .36e-9 # m
#    #micro m
#    a_gamma = .36e-3 #micro m
    NA = 6.022140857e23 #1/mol
    M_Fe = 55.845#g/mol
    c_gt = 4*M_Fe/(a_gamma**3*NA)*(wt_gammatheta(T)/(1-wt_gammatheta(T)))#/(1e-18) #g/(micro m)^3  4*M_Fe/(a_gamma**3*NA)*(wt_gammatheta(T)/(1-wt_gammatheta(T))) #g/m^3
    return c_gt



def c_alphagamma_GT(T, curv):
    Rg = 8.31 #J/(K * mol)
    Eag = 0.39 #J/m^2
    NA = 6.022140857e23 #1/mol
    a_gamma = .36e-9 # m
    Vm = a_gamma**3 / 4 * NA #m^3/mol
    return c_alphagamma(T) * np.exp(Eag * Vm / Rg / T * curv*1e6)

def c_gammaalpha_GT(T, curv):
    Rg = 8.31 #J/(K * mol)
    Ega = 0.39 #Jm/^2
    NA = 6.022140857e23 #1/mol
    a_alpha = .29e-9 # m
    Vm = a_alpha**3 / 2 * NA  #(micro m)**3
    return c_gammaalpha(T) * np.exp(Ega * Vm / Rg / T * curv*1e6)

def c_alphatheta_GT(T, curv):
    Rg = 8.31 #J/(K * mol)
    Eat = 1.05 #J/m^2
    Vm = 1/(7730.14*10**3) #m^3/mol
    return c_alphatheta(T) * np.exp(Eat * Vm / Rg / T * curv*1e6)

def c_gammatheta_GT(T, curv):
    Rg = 8.31 #J/(K * mol)
    Egt = 0.94 #J/m^2
    Vm = 1/(7730.14*10**3) #m^3/mol
    return c_gammatheta(T) * np.exp(Egt * Vm / Rg / T * curv*1e6)



def D_gamma(T):
#    D0_gamma = 1.5e-5 #(micro m)^2/s 1.5e7 #m^2/s
    D0_gamma = 1.5e7 #(micro m)^2/s 1.5e7 #m^2/s
    Q_gamma = 142.1e3/12.011 #J/g
    Rg = 8.31/12.011 #J/(K * g)
    D_g = D0_gamma*np.exp(-Q_gamma/(Rg*T))
#    D_g = 2.98*(1e-6)**2
    return D_g

def D_alpha(T):
#    D0_alpha = 2.2e-4 #(micro m)^2/s 2.2e8 #m^2/s
    D0_alpha = 2.2e8 #(micro m)^2/s 2.2e8 #m^2/s
    Q_alpha = 122.5e3/12.011 #J/g
    Rg = 8.31/12.011 #J/(K * g)
    D_a = D0_alpha*np.exp(-Q_alpha/(Rg*T))
#    D_a = 78.7*(1e-6)**2
    return D_a

def k0_gammaalpha():
    return 1

def k0_alphatheta(T):
    c_gamma_wanted = 62800#90000
    T = T0()
    K_frac = (c_gammatheta(T) - c_gamma_wanted)  / (c_alphatheta(T) - c_alphagamma(T))
    #m
    a_gamma = .36e-9 # m
    a_alpha = .29e-9 # m
#    #micro m
#    a_gamma = .36e-3 #micro m
#    a_alpha = .29e-3 #micro m
    k_frac = K_frac * (a_alpha / a_gamma) * (D_gamma(T) / D_alpha(T))
    return  k0_gammatheta() * k_frac #5e-3#

def k0_alphatheta_K_frac(K_frac):
    #m
    a_gamma = .36e-9 # m
    a_alpha = .29e-9 # m
#    #micro m
#    a_gamma = .36e-3 #micro m
#    a_alpha = .29e-3 #micro m
    return  k0_gammatheta() * K_frac #5e-3#

def k0_gammatheta():
    return 1

def K_gammaalpha(T):
    #m
    a_gamma = .36e-9 # m
    #micro m
    a_gamma = .36e-3 #micro m
    K_ga = k0_gammaalpha()*D_gamma(T)/a_gamma
    return K_ga

def K_alphatheta(T):
    #m
    a_alpha = .29e-9 # m
    #micro m
    a_alpha = .29e-3 #micro m
    K_at = k0_alphatheta(T)*D_alpha(T)/a_alpha
    return K_at

def K_alphatheta_K_frac(T, K_frac):
    #m
    a_alpha = .29e-9 # m
    #micro m
    a_alpha = .29e-3 #micro m
    K_at = k0_alphatheta_K_frac(K_frac)*D_alpha(T)/a_alpha
    return K_at

def K_gammatheta(T):
    #m
    a_gamma = .36e-9 # m
    #micro m
    a_gamma = .36e-3 #micro m
    K_gt = k0_gammatheta()*D_gamma(T)/a_gamma
    return K_gt

def Mchi(T):
    ##        #Sietsma
#    M = 5.4e-8*12.011 #g m/J s
#    chi = 110/12.011 #J/(at%)g
#    print('Sietsma Mchi ', M*chi, 'm/(s at%)')
##    vn_mobility = M*chi*(c_gammaalpha_wt - c_alphagamma_wt)
###    #Mecozzi
    M0 = 2.4e-6 #or e-7 m^4/Js 2D
    chi = 7.4e7 #Jm^-3/wt%
    Q = 140*1e3/12.011
    Rg = 8.31/12.011 #J/(K * g)
    Mchi = M0 * np.exp(-Q/(Rg*T)) * chi
    return Mchi
