#! /usr/bin/python3
# Last edited on 2018-11-11 03:18:52 by stolfilocal

# Computes the ternary phase diagram MgO MgCl2 H2O and writes 
# gnuplot commands for the same.

from math import sqrt, log, hypot
from sys import stderr, stdout

def define_diagram_points():
  """Returns a pair {(g,c)} of dictionaries whose 
  keys are compound labels.  The values of {g} are molas masses
  in grams, and those of {c} are coordinate pairs (X,Y) in the phase diagram.
  The coordinates assigned in {c} assume that the triangle is isosceles 
  with horizontal base."""
  
  # Side length of triangle (arbitrary units).

  g = {}.copy()
  c = {}.copy()
  
  # Data for the diagram in the new reference "phase_diagram.png": */
  L = 678.0; H = -581; x0 = 191; y0 = 651
  
  # Diagram corners:
  g_c_set(g, c, "H_2O",    18.0,  (x0,y0))                     
  g_c_set(g, c, "MgCl_2",  95.2,  (x0+L,y0))
  g_c_set(g, c, "MgO",     40.3,  (x0+L/2, y0+H))

  # Main stable phases at 23 C:
  g_c_combine(g, c, True, "Mg(OH)_2",       ((1, "MgO"), (1, "H_2O")))
  g_c_combine(g, c, True, "MgCl_2.6(H_2O)", ((1, "MgCl_2"), (6, "H_2O")))
  g_c_combine(g, c, True, "P3",             ((3, "Mg(OH)_2"), (1, "MgCl_2"), (8, "H_2O")))
  g_c_combine(g, c, True, "P5",             ((5, "Mg(OH)_2"), (1, "MgCl_2"), (8, "H_2O")))
  
  # Saturated solutions in equilibrium with two phases:
  g_c_combine(g, c, False, "S1", ((0.008, "MgO"), (0.822, "H_2O"), (0.170, "MgCl_2"))) # Sat sol Mg(OH)_2:P5.
  g_c_combine(g, c, False, "S2", ((0.010, "MgO"), (0.768, "H_2O"), (0.222, "MgCl_2"))) # Sat sol P5:P3.
  g_c_combine(g, c, False, "S3", ((0.012, "MgO"), (0.643, "H_2O"), (0.345, "MgCl_2"))) # Sat sol P3:MgCl_2.6(H_2O).

  # Approximate corners of the gel region:
  g_c_combine(g, c, False, "G1", ((0.172, "MgO"), (0.691, "H_2O"), (0.137, "MgCl_2"))) # Gel point Mg(OH)_2:S1.
  g_c_combine(g, c, False, "G2", ((0.184, "MgO"), (0.597, "H_2O"), (0.219, "MgCl_2"))) # Gel Point P3:S2.

  return (g, c)
  
  ######################################################################
  
def define_diagram_lines():
  """Returns a list of pairs {(fi,fj)} of compound labels, whose points
  must be connected by straight lines in the diagram."""
  
  # Side length of triangle (arbitrary units).

  seg = [].copy()
  
  seg = seg + [ ("MgO", "MgCl_2.6(H_2O)") ]
  seg = seg + [ ("MgO", "P5") ]
  seg = seg + [ ("MgO", "P3") ]
  seg = seg + [ ("Mg(OH)_2", "P5") ]
  seg = seg + [ ("Mg(OH)_2", "S1") ]
  seg = seg + [ ("P5", "P3") ]
  seg = seg + [ ("P5", "S1") ]
  seg = seg + [ ("P5", "S2") ]
  seg = seg + [ ("P3", "S2") ]
  seg = seg + [ ("P3", "S3") ]
  seg = seg + [ ("P3", "MgCl_2.6(H_2O)") ]
  seg = seg + [ ("H_2O", "S1") ]
  seg = seg + [ ("S1", "S2") ]
  seg = seg + [ ("S2", "S3") ]
  seg = seg + [ ("S3", "MgCl_2.6(H_2O)") ]

  return seg
  
  ######################################################################
  
def define_diagram_regions():
  """Returns a dictionary that maps each region label to a list of compound labels, 
  which are the corners of the region in the diagram."""
  
  # Side length of triangle (arbitrary units).

  reg = {}.copy()
  
  reg["Sol"] = ( "H_2O", "S1", "S2", "S3", "MgCl_2.6(H_2O)", "H_2O")
  reg["Gel"] = ( "H_2O", "G1", "G2", "MgCl_2.6(H_2O)", "S3", "S2", "S1", "H_2O") 

  return reg
  ######################################################################
  
def g_c_set(g, c, f, g_f, c_f):
  """Assumes {g,c} as in {define_diagram_points}. 
  Defines the molar mass of compound label {f} as {g_f} and its coordinates
  as the 2-tuple {c_f}. Also prints to stderr."""
  
  g[f] = g_f;
  c[f] = c_f;
  stderr.write("%s  M = %.1f  at (%.1f,%.1f)\n" % (f, g[f], c[f][0], c[f][1]))
  ######################################################################
 
def g_c_combine(g, c, molar, f, L):
  """Assumes {g,c} as in {define_diagram_points}. 
  Defines the molar mass and coordinates of compound label {f} by linear combination of the 
  components in list {L}.  Each component is a pair {(mi,fi)} where
  {mi} is a numeric multiplier and {fi} is the label of a compound
  whose molar mass and coordinates are already in {g,c}.  
  
  The numeric multiplier is assumed to be moles if {molar} is true, 
  or a mass fraction if {molar} is false. In the second case, the 
  molar mass {g[f]} is set to {None}. Also prints to stderr."""
  
  c_tot = [0.0, 0.0]
  g_tot = 0.0 # Total mass in grams if {molar} is true, total mass fraction if false.
  for mfi in L:
    mi = mfi[0]
    fi = mfi[1]
    if molar:
      gi = mi*g[fi]  # Assume {mi} is amount of moles.
    else:
      gi = mi  # Assume {mi} is a mass or mass fraction.
    ci = c[fi]
    g_tot = g_tot + gi
    for j in range(2):
      c_tot[j] = c_tot[j] + gi*ci[j]

  if molar:
    g[f] = g_tot
  else:
    g[f] = None
  c[f] = (c_tot[0]/g_tot, c_tot[1]/g_tot)
  stderr.write("%s" % f)
  if g[f] != None:
    stderr.write("  M = %.1f g" % g[f])
  stderr.write("  at (%.1f,%.1f)\n" % (c[f][0], c[f][1]))
  ######################################################################

def choose_point_sizes_and_offsets(c): 
  """Assumes {c} is a point coordinate dict as in {define_diagram_points}.
  
  Returns {psz,off} where {psz} maps compound labels to point sizes,
  {off} maps compound labels to label coordinate offset pairs."""
  
  psz = {}.copy()
  off = {}.copy()
  
  for fi in c.keys():
    if fi[0] == "P":
      psz[fi] = 4.0
    elif fi[0] == "S":
      psz[fi] = 3.0
    elif fi[0] == "G":
      psz[fi] = 0.0
    else:
      psz[fi] = 3.0
   
  off["MgO"] = (+0.15, -0.05)
  off["H_2O"] = (-0.03, +0.14)
  off["MgCl_2"] = (-0.12, -0.10)

  off["MgCl_2.6(H_2O)"] = (+0.15, +0.04)
  off["Mg(OH)_2"] = (-0.12,00.00)

  off["P3"] = (+0.05,00.00)
  off["P5"] = (-0.05,00.00)

  off["S1"] = (-0.04,+0.03)
  off["S2"] = (-0.04,+0.03)
  off["S3"] = (-0.04,+0.03)

  off["G1"] = (-0.04,+0.04)
  off["G2"] = (+0.08,+0.04)

  return psz, off
  ######################################################################
  
def map_coords(c, f, fa, fb, fc):
  """Assumes {c} is a point coordinate dict as in {define_diagram_points}.
  
  Computes the XY plot coordinates of the compoind with label {f}.
  Assumes that {fa,fb,fc} are the formulas of the compounds at
  bottom left, top center, and bottm right corners
  of the diagram, respectively. 
  
  The coordinates in {c} are mapped so that the triangle is 
  equilateral with base corners {(0,0), (1,0)}
  and the apex up."""

  ca = c[fa];
  cb = c[fb];
  cc = c[fc];
  
  x0 = ca[0]; y0 = ca[1]; H = cb[1] - y0; L = cc[0] - x0;
  htri = sqrt(3)/2
  
  x = (c[f][0] - x0)/L
  y = (c[f][1] - y0)/H*htri
  return (x,y)
  ######################################################################


def write_XY_point_file(fname, c, psz, off, moff, fa, fb, fc):
  """Assumes {c} is a point coordinate dict as in {define_diagram_points}.
  
  Assumes {psz} maps compound labels to point sizes,
  {off} maps compound labels to label coordinate offset pairs.
  The {moff} is a scale factor applied to each coordinate of
  each offset in {off}
  
  Writes to "{fname}-pt.txt" the XY point data file of 
  the tri-component phase diagram for all compounds in the dictionary.
  Assumes that {fa,fb,fc} are the formulas of the compounds at
  bottom left, top center, and bottm right corners
  of the diagram, respectively."""
  
  wr = open(fname + "-pt.txt", "w")
  for fi in c.keys():
    (xi,yi) = map_coords(c, fi, fa, fb, fc)
    offi = off[fi]
    pszi = psz[fi]
    xlabi = xi + moff*offi[0]
    ylabi = yi + moff*offi[1]
    wr.write("%6.3f %6.3f  %4.1f" % (xi, yi, pszi))
    wr.write("  %6.3f %6.3f" % (xlabi, ylabi)) 
    wr.write("  %s" % fi)
    wr.write("\n")
  wr.close()
  ######################################################################

def write_XY_line_file(fname, c, seg, fa, fb, fc):
  """Assumes {c} is a point coordinate dict as in {define_diagram_points}.
  
  Assumes that {seg} is a list of pairs of compound labels to be connected 
  by straight lines.
  
  Writes to "{fname}-pt.txt" the  point data file of 
  the tri-component phase diagram for all compounds in the dictionary.
  Assumes that {fa,fb,fc} are the formulas of the compounds at
  bottom left, top center, and bottm right corners
  of the diagram, respectively."""
  
  wr = open(fname + "-ln.txt", "w")
  for fij in seg:
    fi = fij[0]
    fj = fij[1]
    (xi, yi) = map_coords(c, fi, fa, fb, fc)
    (xj, yj) = map_coords(c, fj, fa, fb, fc)
    
    wr.write("%6.3f %6.3f\n" % (xi, yi))
    wr.write("%6.3f %6.3f\n" % (xj, yj))
    wr.write("\n")
  wr.close()
  ######################################################################

def write_XY_region_files(fname, c, reg, fa, fb, fc):
  """Assumes {c} is a point coordinate dict as in {define_diagram_points}.
  
  Assumes that {reg} is a dict that maps region labels to vertex lists,
  as in {define_diagram_regions}.
  
  For each region {ri} in {reg}, writes to "{fname}-rg-{ri}.txt" the 
  coordinates of the vertices of the region {ri}. Assumes that {fa,fb,fc}
  are the formulas of the compounds at bottom left, top center, and 
  bottm right corners of the diagram, respectively."""
  
  for rk in reg.keys():
    wr = open(fname + "-rg-" + rk + ".txt", "w")
    vk = reg[rk] # Vertex label list.
    for fi in vk:
      (xi, yi) = map_coords(c, fi, fa, fb, fc)
      wr.write("%6.3f %6.3f\n" % (xi, yi))
    wr.close()
  ######################################################################
  
(g,c) = define_diagram_points()
seg = define_diagram_lines()
reg = define_diagram_regions()

(psz,off) = choose_point_sizes_and_offsets(c)

write_XY_point_file("main", c, psz, off, 0.5, "H_2O", "MgO", "MgCl_2")
write_XY_line_file("main", c, seg, "H_2O", "MgO", "MgCl_2")
write_XY_region_files("main", c, reg, "H_2O", "MgO", "MgCl_2")