#Input
T = 923.15
R = 8.314
C =0.005
phi = 0.5
# 
from fipy.tools import numerix
#--------------------------------------------------------
# calculation
#--------------------------------------------------------
RT = R*T 
term1 = 0.02*1.e-4*numerix.exp(-10115.0/T)*numerix.exp(0.5898 \
                                               *(1.0+2.0/numerix.pi*numerix.arctan(14.985-15309.0/T)))
term2 = 4.529e-7*numerix.exp(-(1.0/T-2.221e-4)*17767.0)
term3 = term2/term1
term4 = (1.0/T-2.221e-4)*26436.0
M  = term3*numerix.exp(C*term4)
p = phi**3*(10.0 - 15.0*phi + 6.0*phi**2)
Mp = M**p
print 'M**p = ', Mp