#! /usr/bin/python3
# Last edited on 2021-11-25 21:16:04 by stolfi

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

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

# The following procedures use two dictionaries {g} and {c}
# whose keys are compound labels (strings).
#
# The value of {g[lab]} is the molar mass of compound {lab} in grams.
# 
# The value {c[lab]} is the Cartesian coordinates {(X,Y)}
# of the corresponding point on the diagram. They assume that the triangle is isosceles 
# with horizontal base.

def g_c_set(g, c, f, g_f, c_f):
  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):
  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 map_coords(c, f, fa, fb, fc):
  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):
  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):
  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):
  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()
  # ----------------------------------------------------------------------