import os
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 500
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline 
%config InlineBackend.figure_format = 'retina'

"""
Explanation: Simulating CSTR
Here a CSTR model simulation is demoed with pure surface phase mechanism focussed on Ammonia dehydrogenation. The mechanims consists of 18 species and 7 reactions. In conventional terms, Ammonia dehydrogenation can written as 2NH3 -> N2 + 3H2 
In this demo, we showcase the effect of different operating conditions on the surface coverages and molar fractions. 
Note: The actual simulations are executed from command line. In the jupyter notebook, analysis of the resulting data is carried out.
End of explanation
"""


ls T-950K/

mole_df = pd.read_csv(os.path.join('T-950K','gas_mole_tr.csv'))
msdot_df = pd.read_csv(os.path.join('T-950K','gas_msdot_tr.csv'))
surf_cov_df = pd.read_csv(os.path.join('T-950K','surf_cov_tr.csv'))

"""
Explanation: Influence of temperature
Evaluate the effect of temperatures of 950 K and then 900 K on the NH3 dehydrogenation
End of explanation
"""


ax = plt.subplot(1, 1, 1)
ax.plot('t(s)', 'H2', color='r', data=mole_df, marker='^', markersize=0.5, label="$H_2$")
ax.plot('t(s)', 'N2', color='g', data=mole_df, marker='v', markersize=0.5, label="$N_2$")
ax.plot('t(s)', 'NH3', color='b', data=mole_df, marker='*', markersize=0.5, label="$NH_3$")
ax.set_xlabel('Time (s)')
ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax.set_ylabel('Mole Fraction')
#ax.set_xlim([0,50])
# Change time limit to 1
ax.set_xlim([0,0.04])
ax.legend(loc="upper left", bbox_to_anchor=(1,1))

plt.tight_layout()
plt.show()
#plt.savefig('Ammonia_gas_species_molefrac.png', dpi=500)

"""
Explanation: Mole fractions
End of explanation
"""


ax3 = plt.subplot(1, 1, 1)
ax3.plot('t(s)', 'H2', color='r', data=msdot_df, marker='^', markersize=0.5, label="$H_2$")
ax3.plot('t(s)', 'N2', color='g', data=msdot_df, marker='v', markersize=0.5, label="$N_2$")
ax3.plot('t(s)', 'NH3', color='b', data=msdot_df, marker='*', markersize=0.5, label="$NH_3$")
ax3.set_xlabel('Time (s)')
ax3.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax3.set_ylabel('Surface Production Rate')
#ax3.set_xlim([0,50])
ax3.set_xlim([0,0.001])
ax3.set_ylim([-1e-6,1e-6])
ax3.legend(loc="upper left", bbox_to_anchor=(1,1))

plt.tight_layout()
plt.show()
#plt.savefig('Ammonia_gas_species_production_rate.png', dpi=500)

print(surf_cov_df.columns)
ax1 = plt.subplot(1, 1, 1)
ax1.plot('t(s)', 'RU(S1)', color='r', data=surf_cov_df, marker='^', markersize=0.5, label="$Ru*$")
ax1.plot('t(s)', 'H(S1)', color='g', data=surf_cov_df, marker='v', markersize=0.5, label="$H*$")
ax1.plot('t(s)', 'N(S1)', color='b', data=surf_cov_df, marker='*', markersize=0.5, label="$N*$")
ax1.set_xlabel('Time (s)')
ax1.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax1.set_ylabel('Surface Coverages')
ax1.set_xlim([0,5])
#ax1.set_xlim([0,0.04])
ax1.legend(loc="upper left", bbox_to_anchor=(1,1))

plt.tight_layout()
plt.show()
#plt.savefig('Ammonia_surf_cov.png', dpi=500)

"""
Explanation: Surface Production Rates
End of explanation
"""


mole_df_950 = pd.read_csv(os.path.join('T-950K','gas_mole_tr.csv'))
msdot_df_950 = pd.read_csv(os.path.join('T-950K','gas_msdot_tr.csv'))
surf_cov_df_950 = pd.read_csv(os.path.join('T-950K','surf_cov_tr.csv'))

mole_df_900 = pd.read_csv(os.path.join('T-900K','gas_mole_tr.csv'))
msdot_df_900 = pd.read_csv(os.path.join('T-900K','gas_msdot_tr.csv'))
surf_cov_df_900 = pd.read_csv(os.path.join('T-900K','surf_cov_tr.csv'))

ax = plt.subplot(1, 1, 1)
ax.plot('t(s)', 'H2', '-', color='r', data=mole_df_950,  marker='^', markersize=0.5, label="$H_2$ 950K")
ax.plot('t(s)', 'H2', '--', color='r', data=mole_df_900, marker='^', markersize=0.5, label="$H_2$ 900K")

ax.plot('t(s)', 'N2', '-', color='g', data=mole_df_950,  marker='v', markersize=0.5, label="$N_2$ 950K")
ax.plot('t(s)', 'N2', '--', color='g', data=mole_df_900, marker='v', markersize=0.5, label="$N_2$ 900K")

ax.plot('t(s)', 'NH3', '-', color='b', data=mole_df_950,  marker='*', markersize=0.5, label="$NH_3$ 950K")
ax.plot('t(s)', 'NH3', '--', color='b', data=mole_df_900, marker='*', markersize=0.5, label="$NH_3$ 900K")


ax.set_xlabel('Time (s)')
ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax.set_ylabel('Mole Fraction')
#ax.set_xlim([0,50])
# Change time limit to 1
ax.set_xlim([0,5])
ax.legend(loc="upper left", bbox_to_anchor=(1,1))

plt.tight_layout()
plt.savefig('Ammonia_molefrac_Tcomp.png', dpi=500)

ax = plt.subplot(1, 1, 1)
ax.plot('t(s)', 'RU(S1)', '-', color='r', data=surf_cov_df_950, marker='^', markersize=0.5, label="$Ru*$ 950K")
ax.plot('t(s)', 'RU(S1)', '--', color='r', data=surf_cov_df_900, marker='^', markersize=0.5, label="$Ru*$ 900K")

ax.plot('t(s)', 'H(S1)', '-', color='g', data=surf_cov_df_950, marker='v', markersize=0.5, label="$H*$ 950K")
ax.plot('t(s)', 'H(S1)', '--', color='g', data=surf_cov_df_900, marker='v', markersize=0.5, label="$H*$ 900K")

ax.plot('t(s)', 'N(S1)', '-', color='b', data=surf_cov_df_950, marker='*', markersize=0.5, label="$N*$ 950K")
ax.plot('t(s)', 'N(S1)', '--', color='b', data=surf_cov_df_900, marker='*', markersize=0.5, label="$N*$ 900K")

ax.set_xlabel('Time (s)')
ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax.set_ylabel('Surface Coverages')
#ax.set_xlim([0,0.5])
ax.set_xlim([0,0.04])
ax.legend(loc="upper left", bbox_to_anchor=(1,1))

plt.tight_layout()
plt.savefig('Ammonia_surf_cov_Tcomp.png', dpi=500)

"""
Explanation: Comparison between temperatures
End of explanation
"""


mole_df_10 = pd.read_csv(os.path.join('AbyV-1000','gas_mole_tr.csv'))
msdot_df_10 = pd.read_csv(os.path.join('AbyV-1000','gas_msdot_tr.csv'))
surf_cov_df_10 = pd.read_csv(os.path.join('AbyV-1000','surf_cov_tr.csv'))

mole_df_15 = pd.read_csv(os.path.join('AbyV-1500','gas_mole_tr.csv'))
msdot_df_15 = pd.read_csv(os.path.join('AbyV-1500','gas_msdot_tr.csv'))
surf_cov_df_15 = pd.read_csv(os.path.join('AbyV-1500','surf_cov_tr.csv'))

ax = plt.subplot(1, 1, 1)
ax.plot('t(s)', 'H2', '-', color='r', data=mole_df_15, marker='^', markersize=0.5, label="$H_2$ 1.5k")
ax.plot('t(s)', 'H2', '--', color='r', data=mole_df_10, marker='^', markersize=0.5, label="$H_2$ 1.0K")

ax.plot('t(s)', 'N2', '-', color='g', data=mole_df_15, marker='v', markersize=0.5, label="$N_2$ 1.5k")
ax.plot('t(s)', 'N2', '--', color='g', data=mole_df_10, marker='v', markersize=0.5, label="$N_2$ 1.0k")

ax.plot('t(s)', 'NH3', '-', color='b', data=mole_df_15, marker='*', markersize=0.5, label="$NH_3$ 1.5k")
ax.plot('t(s)', 'NH3', '--', color='b', data=mole_df_10, marker='*', markersize=0.5, label="$NH_3$ 1.0k")


ax.set_xlabel('Time (s)')
ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax.set_ylabel('Mole Fraction')
#ax.set_xlim([0,50])
# Change time limit to 1
ax.set_xlim([0,5])
ax.legend(loc="upper left", bbox_to_anchor=(1,1))

plt.tight_layout()
plt.show()
#plt.savefig('Ammonia_molefrac_catload_cmp.png', dpi=500)

ax = plt.subplot(1, 1, 1)
ax.plot('t(s)', 'RU(S1)', '-', color='r', data=surf_cov_df_15, marker='^', markersize=0.5, label="$Ru*$ 1.5k")
ax.plot('t(s)', 'RU(S1)', '--', color='r', data=surf_cov_df_10, marker='^', markersize=0.5, label="$Ru*$ 1.0k")

ax.plot('t(s)', 'H(S1)', '-', color='g', data=surf_cov_df_15, marker='v', markersize=0.5, label="$H*$ 1.5k")
ax.plot('t(s)', 'H(S1)', '--', color='g', data=surf_cov_df_10, marker='v', markersize=0.5, label="$H*$ 1.0k")

ax.plot('t(s)', 'N(S1)', '-', color='b', data=surf_cov_df_15, marker='*', markersize=0.5, label="$N*$ 1.5k")
ax.plot('t(s)', 'N(S1)', '--', color='b', data=surf_cov_df_10, marker='*', markersize=0.5, label="$N*$ 1.0k")

ax.set_xlabel('Time (s)')
ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax.set_ylabel('Surface Coverages')
#ax.set_xlim([0,0.5])
ax.set_xlim([0,0.04])
ax.legend(loc="upper left", bbox_to_anchor=(1,1))

plt.tight_layout()
plt.show()
#plt.savefig('Ammonia_surf_cov_catload-comp.png', dpi=500)

"""
Explanation: Effect of catalyst loading
Two different catalyst loadings 1500/cm A/V and 1000/cm A/V are evaluated.
End of explanation
"""


