[cig-commits] [commit] add_thermodynamic_potentials: Added example plot to test_process_solidsolution.py (8a557c9)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/burnman

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51

>---------------------------------------------------------------

commit 8a557c96432ceb1fed0b689f2acc18333bd77652
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Sat Aug 30 00:29:15 2014 +0200

    Added example plot to test_process_solidsolution.py


>---------------------------------------------------------------

8a557c96432ceb1fed0b689f2acc18333bd77652
 burnman/test_process_solidsolution.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/burnman/test_process_solidsolution.py b/burnman/test_process_solidsolution.py
index f4e66e6..dda3bff 100644
--- a/burnman/test_process_solidsolution.py
+++ b/burnman/test_process_solidsolution.py
@@ -153,3 +153,26 @@ print np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Wh,phi)), 'J/mol
 print np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Ws,phi)), 'J/K/mol'
 
 print np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Wv,phi)), 'm^3/mol'
+
+
+
+# Plot excess volumes for the pyrope-grossular join
+
+ppy= np.linspace(0, 1, 101)
+vex= np.empty(shape=(101))
+for p in range(len(ppy)):
+    a=ppy[p]
+    endmember_proportions = np.array([ a, 0.0, 1.-a, 0.0 ])
+    phi=np.array([alpha[i]*endmember_proportions[i] for i in range(n_endmembers)])
+    phi=np.divide(phi, np.sum(phi))
+
+    vex[p]=np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Wv,phi))
+
+
+import matplotlib.pyplot as plt
+plt.plot(ppy,vex,color='r',linestyle='-',marker='o',markerfacecolor='r',markersize=0)
+plt.xlim(min(ppy),max(ppy))
+plt.xlabel("p(pyrope)")
+plt.title("V excess (m^3/mol) for pyrope-grossular garnets")
+
+plt.show()