#---------------------Description--------------------------------------
#
#      Script to compute pCO2 in silicate/carbonate equilibrium
#       The Urey reactions considered are:
#          (1) MgSiO3 + CO2 = MgCO3 + SiO2
#          (2) CaSiO3 + CO2 = CaCO3 + SiO2
#          (3) FeSiO3 + CO2 = FeCO3 + SiO2
#       and the thermodynamic coefficients are from:
#           Adamcik and Draper (1963) Planet Space Sci. 11, pp1303-1307.
#
#
#-----------------------------------------------------------------------
#ToDo:

#

#Data on section of text which this script is associated with
Chapter = '8'
Section = ''
Figure = 'fig:UreyEq'
Table = ''
#

import phys
import math
from ClimateUtilities import *

#Free energy (298K), entropy (298K) and temperature coeffs
MgCoeffs = [79496.e3,174.5264e3,22.09e3,-63.43,-6.19e8]
CaCoeffs = [88700.8e3,161.08e3,4.2248e3,-32.384,-1.339e8]
FeCoeffs = [64852.e3,174.05e3,48.49e3,-100.91,-16.99e8]

#Answer is given in bars
def pCO2(T,coeffs):
    Rs = phys.Rstar
    dH,dS,A,B,C = coeffs
    logK = -dH/(Rs*T) + dS/Rs
    logK += A*( (298.-T)/(Rs*T) + math.log(T/298.)/Rs)
    logK += B*( (298.**2-T**2)/(2.*Rs*T) + (T-298.)/Rs)
    logK += C*((298.-T)/(298.*Rs*T**2) + (T**2-298.**2)/(2.*298.**2*Rs*T**2))
    return math.exp(logK)

#Make plot

c= Curve()
c.XlogAxis = c.YlogAxis = True
c.Ylabel = "pCO2 (bars)"
c.Xlabel = "T (K)"
T = [280.+2.*i for i in range(220)]
pCO2_Mg = [pCO2(TT,MgCoeffs) for TT in T]
pCO2_Ca = [pCO2(TT,CaCoeffs) for TT in T]
pCO2_Fe = [pCO2(TT,FeCoeffs) for TT in T]
c.addCurve(T,'T')
c.addCurve(pCO2_Mg,'Mg')
c.addCurve(pCO2_Ca,'Ca')
c.addCurve(pCO2_Fe,'Fe')
plot(c)