#! /usr/bin/python3
# Last edited on 2022-07-20 13:55:26 by stolfi

# Computes the ternary phase diagram Na2O-B2O3-H2O and writes 
# gnuplot commands for the same.

import ternary_phase_diagram as tpd
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 molar 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:
  tpd.g_c_set(g, c, "H_2O",    18.0,  (x0+L/2,y0+H))  # Apex of diagram.                  
  tpd.g_c_set(g, c, "Na_2O",   62.0,  (x0,y0))        # Left base corner.                
  tpd.g_c_set(g, c, "B_2O_3",  69.6,  (x0+L,y0))      # Right base corner.                

  # Main stable phases at 23 C:
  tpd.g_c_combine(g, c, True, "B(OH)_3",       ((1, "B_2O_3"), (3, "H_2O")))
  tpd.g_c_combine(g, c, True, "NaB_5(OH)_3",       ((1, "B_2O_3"), (3, "H_2O")))
  ??

  tpd.g_c_combine(g, c, True, "Mg(OH)_2",       ((1, "MgO"), (1, "H_2O")))
  tpd.g_c_combine(g, c, True, "MgCl_2.6(H_2O)", ((1, "MgCl_2"), (6, "H_2O")))
  tpd.g_c_combine(g, c, True, "P3",             ((3, "Mg(OH)_2"), (1, "MgCl_2"), (8, "H_2O")))
  tpd.g_c_combine(g, c, True, "P5",             ((5, "Mg(OH)_2"), (1, "MgCl_2"), (8, "H_2O")))
  
  # Saturated solutions in equilibrium with two phases:
  tpd.g_c_combine(g, c, False, "S1", ((0.008, "MgO"), (0.822, "H_2O"), (0.170, "MgCl_2"))) # Sat sol Mg(OH)_2:P5.
  tpd.g_c_combine(g, c, False, "S2", ((0.010, "MgO"), (0.768, "H_2O"), (0.222, "MgCl_2"))) # Sat sol P5:P3.
  tpd.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:
  tpd.g_c_combine(g, c, False, "G1", ((0.172, "MgO"), (0.691, "H_2O"), (0.137, "MgCl_2"))) # Gel point Mg(OH)_2:S1.
  tpd.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 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
  # ----------------------------------------------------------------------
  
  
(g,c) = define_diagram_points()
seg = define_diagram_lines()
reg = define_diagram_regions()

(psz,off) = choose_point_sizes_and_offsets(c)

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