'''
MagmaForge: Mantle Melting Pseudosection
========================================

Demonstrates how to run a fractional (isentropic) decompression routine using MagmaForge.

Open this code in an executable MyBinder instance (MyBinder links may be slow to load-- please be patient!):

.. image:: https://mybinder.org/badge_logo.svg
  :target: https://mybinder.org/v2/gl/swmatthews-research%2FThermoEngineLite/main?urlpath=%2Fdoc%2Ftree%2F.%2Fdoc%2Fsource%2Fauto_examples%2F_2_magmaforge%2Fplot_magmaforge_decompression_solidus.ipynb
'''

# %% 
# Initialization
# --------------
#
# Import necessary packages:

from thermoengine import magmaforge

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# %%
# Define a bulk composition
# -------------------------

DMM_comp_adj = {
    'SiO2': 44.7354,
    'TiO2': 0.130074,
    'Al2O3': 3.98226,
    'Fe2O3': 0.275456,
    'Cr2O3': 0.570323,
    'FeO': 8.01254,
    'MgO': 38.752,
    'CaO': 3.1718,
    'Na2O': 0.130074
}

# %%
# Run crystallization at different P slices
# -----------------------------------------

state_histories = []
P_vals  = np.linspace(2,0.5,10)
for P_GPa in P_vals:
    # Start at liquidus
    sys = magmaforge.System(comp=DMM_comp_adj, 
                        T_liquidus=True,
                        P_GPa=P_GPa, 
                        database='MELTS_pMELTS',
                        )
    # Crystallize to solidus
    sys.crystallize(method='equil',
                    T_step=5,
                    fix_fO2=False,
                    calc_args={'debug':0},
                    #melt_frac_cutoff=0.001
                    )
    
    state_histories.append(sys.history)

# %%
# Calculate Liquidus, Solidus, and Melt Isopleths
# -----------------------------------------------

T_liquidus = np.zeros_like(P_vals)
T_solidus = np.zeros_like(P_vals)
F_vals = np.linspace(0.10,0.90,9)
T_F_array = np.zeros([len(F_vals), len(P_vals)])
for P_idx,state_history in enumerate(state_histories):
    T_liquidus[P_idx] = state_history.get_temperatures('C')[0]
    T_solidus[P_idx] = state_history.get_temperatures('C')[-1]

    F = state_history.get_melt_fractions()
    for F_idx, F_val in enumerate(F_vals):
        T_idx = np.abs(F - F_val).argmin()
        T_F_array[F_idx,P_idx] = state_history.get_temperatures('C')[T_idx]


# %%
# Create Plot
# -----------

fig, ax = plt.subplots(nrows=1, figsize=(4,4))

ax.plot(T_liquidus, P_vals, lw=2, color='r')
ax.plot(T_solidus, P_vals, lw=2, color='k')
ax.plot(T_F_array.T, P_vals, lw=1, ls=':', color='k')

ax.fill_betweenx(P_vals,1100,T_solidus,color='black', alpha=0.2)
ax.fill_betweenx(P_vals,1800,T_liquidus,color='red', alpha=0.2)

ax.set_xlabel('T (°C)')
ax.set_xlim([1450-273.15,2000-273.15])

ax.set_ylabel('P (GPa)')
ax.set_ylim([P_vals.max(), P_vals.min()])

ax.text(1550-273.15,1.5,'solid')
ax.text(1900-273.15,1.0,'liquid', color='r')
ax.text(1290,0.48,'melt fraction isopleths', color='k', fontsize=8)

ax.set_title('Partial Melting of DMM', pad=20)

plt.tight_layout()
plt.show()
# %%