"""
eafEnergy.py
Copyright 2016, Nikolas Martelaro
http://large.stanford.edu/courses/2016/ph240/martelaro1
Released for public domain use with no restrictions.

Computes the energy to heat and melt steel in an Electric arc furnace

Usage: python eafEnergy.py

Assumptions:
Steel assumed to be pure iron
"""

### CONSTANTS ###
M_ATM_IRON = 55.845             # atomic mass of iron [ g/mol ]
M_IRON = 1e6                    # mass of iron = 1 ton


### Iron goes through three phases, alpha, gamma, and delta before melting. To
### calculate the total energy to melt steel, we compute the heating for each
### phase as well as the heat required for each phase change.

### Temperatures [ K ] ###
T_alpha = [298.15, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, \
           900, 925, 950, 975, 1000, 1010, 1020, 1030, 1035, 1040, 1042, 1043, \
           1044, 1046, 1048, 1050, 1055, 1060, 1070, 1080, 1090, 1100, 1125, \
           1150, 1175, 1185]
T_gamma = [1185, 1200, 1250, 1300, 1350, 1400, 1450, 1500, 1550, 1600, 1650, \
           1667]
T_delta = [1667, 1700, 1750, 1800, 1811]

### Specific Heats [ J/(mol*K) ] ###
C_alpha = [25.084, 25.131, 26.321, 27.472, 28.532, 29.639, 30.810, 32.008, \
           33.260, 34.618, 36.153, 37.924, 40.170, 43.157, 45.152, 47.642, \
           50.691, 54.458, 56.815, 60.140, 65.490, 69.420, 74.900, 78.690, \
           83.770, 74.540, 68.270, 63.980, 60.900, 56.348, 53.850, 50.138, \
           47.462, 45.675, 44.350, 42.223, 41.063, 40.262, 40.000]
C_gamma = [33.775, 33.905, 34.353, 34.809, 35.280, 35.750, 36.220, 36.690, \
           37.161, 37.630, 38.100, 38.260]
C_delta = [40.400, 41.454, 43.051, 44.649, 45.000]

### Specific Enthalphy [ J/mol ] ###
H_alpha_gamma = 900
H_gamma_delta = 850
H_delta_liquid = 13810


### Energy of heating ###
Q = 0

# Alpha phase
for i in range(0, len(T_alpha) - 1):
    Q += M_IRON * C_alpha[i]/M_ATM_IRON * (T_alpha[i+1] - T_alpha[i])

# Gamma phase
for i in range(0, len(T_gamma) - 1):
    Q += M_IRON * C_gamma[i]/M_ATM_IRON * (T_gamma[i+1] - T_gamma[i])

# Delta phase
for i in range(0, len(T_delta) - 1):
    Q += M_IRON * C_delta[i]/M_ATM_IRON * (T_delta[i+1] - T_delta[i])

# Energy of phase changes
Q += (H_alpha_gamma / M_ATM_IRON) * M_IRON
Q += (H_gamma_delta / M_ATM_IRON) * M_IRON
Q += (H_delta_liquid/ M_ATM_IRON) * M_IRON

print("Energy required to melt steel in an electric arc furnace")
print("Total energy = %.3E J" % Q)

# REFERENCES
# [1] http://nist.gov/data/PDFfiles/jpcrd298.pdf