root/trunk/matml/src/ternary/Ternary.c

/***************************************
  $Header$

  This is a Glade front end to Ternary GUI.

  Copyright 2009 Adam C. Powell, IV
  ***************************************/


#include "gibbs.h"
#include <glade/glade.h>
#include <libgnomeui/libgnomeui.h>

GladeXML *xml;
FILE *isotherm_geomview_pipe=NULL, *phdiag_geomview_pipe=NULL;

/* eparams: name, R,T0, G1@T0,G2,G3, C1,C2,C3, S1,S2,S3,S4,S5,
   O12,O13,O23, O123,O234,O235,O245,O345,O2345 gauss ptr, gauss number */
energy_params allparams[2];
double T_iso=298., P=1., T_min=273., T_max=373.;
int P1_res=20, P2_res=20, T_steps=50, numpoints=462, refine=2, *verts=NULL;
int max_boundaries=0;
gboolean spinodal_toggle [2];

char P1_name[100], P2_name[100], gv_version [100], *entry_names[] = {
  "T_entry", "R_entry", "P_entry", "Tmin_entry", "Tmax_entry",
  "P1_T_entry", "P1_G1_entry", "P1_G2_entry", "P1_G3_entry",
  "P1_C1_entry", "P1_C2_entry", "P1_C3_entry",
  "P1_S1_entry", "P1_S2_entry", "P1_S3_entry",
  "P1_O12_entry", "P1_O13_entry", "P1_O23_entry", "P1_O123_entry",
  "P2_T_entry", "P2_G1_entry", "P2_G2_entry", "P2_G3_entry",
  "P2_C1_entry", "P2_C2_entry", "P2_C3_entry",
  "P2_S1_entry", "P2_S2_entry", "P2_S3_entry",
  "P2_O12_entry", "P2_O13_entry", "P2_O23_entry", "P2_O123_entry" };

char *hscale_names[] = {
  "T_hscale", "R_hscale", "P_hscale", "Tmin_hscale", "Tmax_hscale",
  "P1_T_hscale", "P1_G1_hscale", "P1_G2_hscale", "P1_G3_hscale",
  "P1_C1_hscale", "P1_C2_hscale", "P1_C3_hscale",
  "P1_S1_hscale", "P1_S2_hscale", "P1_S3_hscale",
  "P1_O12_hscale", "P1_O13_hscale", "P1_O23_hscale", "P1_O123_hscale",
  "P2_T_hscale", "P2_G1_hscale", "P2_G2_hscale", "P2_G3_hscale",
  "P2_C1_hscale", "P2_C2_hscale", "P2_C3_hscale",
  "P2_S1_hscale", "P2_S2_hscale", "P2_S3_hscale",
  "P2_O12_hscale", "P2_O13_hscale", "P2_O23_hscale", "P2_O123_hscale" };

GtkWidget *entry_widgets [33], *hscale_widgets [33];
double *pointers[] =
  { &T_iso, &(allparams[0].R), &P, &T_min, &T_max,
    &(allparams[0].T0), &(allparams[0].G1_T0), &(allparams[0].G2_T0), &(allparams[0].G3_T0),
    &(allparams[0].G1_C), &(allparams[0].G2_C), &(allparams[0].G3_C),
    &(allparams[0].S1), &(allparams[0].S2), &(allparams[0].S3),
    &(allparams[0].Omega12), &(allparams[0].Omega13), &(allparams[0].Omega23),
    &(allparams[0].Omega123),
    &(allparams[1].T0), &(allparams[1].G1_T0), &(allparams[1].G2_T0), &(allparams[1].G3_T0),
    &(allparams[1].G1_C), &(allparams[1].G2_C), &(allparams[1].G3_C),
    &(allparams[1].S1), &(allparams[1].S2), &(allparams[1].S3),
    &(allparams[1].Omega12), &(allparams[1].Omega13), &(allparams[1].Omega23),
    &(allparams[1].Omega123) };

void update_pointer (int i, double content)
{
  *pointers [i] = content;

  if (i==1)
    allparams[1].R = content;
}

void recalculate_and_display ()
{
  int i, numbounds, hullnumfacets;
  energy_point *points=NULL;
  hull_facet *hullfacets = NULL;
  phase_boundary spinodals [2], *allbounds=NULL;
  visual_point *vpoints=NULL;

  // Reallocate points array, and set up triangles
  numpoints = (P1_res+1)*(P1_res+2)/2 + (P2_res+1)*(P2_res+2)/2;
  if (!(points = (energy_point *) realloc
        (points, numpoints * sizeof (energy_point))))
    { printf ("recalculate_and_display: "
              "Cannot (re-)allocate point coordinate array\n"); exit (1); }
  if (!(verts = (int *) realloc
        (verts, (P1_res*P1_res*3 + P2_res*P2_res*3) * sizeof (int))))
    { printf ("recalculate_and_display: "
              "Cannot (re-)allocate triangle vertex array\n"); exit (1); }
  if (i=init_triangle_array (P1_res, points, 0, 0.,0., 1.,0., 0.,1.))
    { printf ("recalculate_and_display: Error %d in init_triangle_array\n", i);
      exit (i); }
  if (i=init_triangle_array (P2_res, points+(P1_res+1)*(P1_res+2)/2,
                             0, 0.,0., 1.,0., 0.,1.))
    { printf ("recalculate_and_display: Error %d in init_triangle_array\n", i);
      exit (i); }
  if (i=init_triangle_vertices (P1_res, verts, 0))
    { printf ("recalculate_and_display: "
              "Error %d in init_triangle_vertices\n", i); exit (i); }
  if (i=init_triangle_vertices (P2_res, verts+P1_res*P1_res*3,
                                (P1_res+1)*(P1_res+2)/2))
    { printf ("recalculate_and_display: "
              "Error %d in init_triangle_vertices\n", i); exit (i); }

  // Calculate free energies and convex hull, and refine
  if (i=free_energies (points, numpoints/2, T_iso,P, allparams,0,
                       WITH_DERIVATIVES_NO_INFINITY))
    { printf ("recalculate_and_display: Error %d in free_energies\n", i);
      exit (i); }
  if (i=free_energies (points+numpoints/2, numpoints/2, T_iso,P, allparams,1,
                       WITH_DERIVATIVES_NO_INFINITY))
    { printf ("recalculate_and_display: Error %d in free_energies\n", i);
      exit (i); }

  if (i=hullCalculate (3, points, numpoints, allparams, 2, T_iso, P,
                       &hullfacets, &hullnumfacets))
    { printf ("recalculate_and_display: Error %d in hullCalculate\n", i);
      exit (i); }
  for (i=0; i<refine; i++)
    hullRefine (&points, &numpoints, &hullfacets, &hullnumfacets,
                allparams, 2, T_iso, P, 1e-10, 1e-7, 0);

  // Calculate spinodals if appropriate
  if (spinodal_toggle [0])
    if (i=calc_spinodal
        (&points, &numpoints, verts, P1_res*P1_res, T_iso,P, allparams,0,
         spinodals))
      { printf ("recalculate_and_display: Error %d in calc_spinodal\n", i);
        exit (i); }
  if (spinodal_toggle [1])
    if (i=calc_spinodal
        (&points, &numpoints, verts+P1_res*P1_res, P2_res*P2_res, T_iso,P,
         allparams,1, spinodals+1))
      { printf ("recalculate_and_display: Error %d in calc_spinodal\n", i);
        exit (i); }

  // Calculate phase boundaries
  if (i=hullReturnPhaseBoundaries (points,numpoints, hullfacets,hullnumfacets,
                                   allparams, T_iso, P, &allbounds, &numbounds))
    { printf ("recalculate_and_display: "
              "error %d in hullReturnPhaseBoundaries\n", i); exit (i); }

  // Display free energy surface, spinodals and boundaries using geomview
  if (isotherm_geomview_pipe)
    {
      if (i=energy_to_visual_array (numpoints, points, &vpoints, 0.,-1500.,
                                    1.,500., NULL,NULL, 0))
        { printf ("recalculate_and_display: "
                  "Error %d in energy_to_visual_array\n", i); exit (i); }

      if (i=GeomviewDisplayTriangleCOFF
          (isotherm_geomview_pipe, "Free Energy Function", "tfe", "shading smooth",
           numpoints, vpoints, P1_res*P1_res + P2_res*P2_res, verts))
        { printf ("recalculate_and_display: Error %d in Geomview Display\n", i);
          exit (i); }

      if (i=GeomviewDisplayPhaseBoundary
          (isotherm_geomview_pipe, "Spinodal 1", "spin1", "-face +edge",
           numpoints, vpoints, spinodal_toggle[0] ? spinodals : NULL))
        { printf ("recalculate_and_display: Error %d in Geomview Display\n", i);
          exit (i); }
      if (i=GeomviewDisplayPhaseBoundary
          (isotherm_geomview_pipe, "Spinodal 2", "spin2", "-face +edge",
           numpoints, vpoints, spinodal_toggle[1] ? spinodals+1 : NULL))
        { printf ("recalculate_and_display: Error %d in Geomview Display\n", i);
          exit (i); }

      max_boundaries = (numbounds>max_boundaries) ? numbounds : max_boundaries;
      for (i=0; i<max_boundaries; i++)
        {
          char name [20], label [10];
          int ret;

          snprintf (name, 18, "Boundary %d", i);
          snprintf (label, 8, "bd%d", i);
          if (ret=GeomviewDisplayPhaseBoundary
              (isotherm_geomview_pipe, name, label, "-face +edge", 
               numpoints, vpoints, (i<numbounds) ? allbounds+i : NULL))
            { printf ("recalculate_and_display: "
                      "Error %d in Geomview Display\n", ret); exit (ret); }
        }
    }

  free (vpoints);
  free (allbounds);
  free (hullfacets);
  free (points);

  //printf ("New parameters: \n");
  //for (i=0; i<33; i++)
  //  printf (" %s %g", entry_names[i], *(pointers[i]));
  //printf ("\n");
}

void get_entry_update_hscale_and_pointer (int i)
{
  GtkWidget *entry = entry_widgets [i], *hscale = hscale_widgets [i];
  G_CONST_RETURN char *entrytext;
  gdouble content;

  entrytext = gtk_entry_get_text (GTK_ENTRY(entry));
  sscanf (entrytext, "%lf", &content);

  gtk_range_set_value (GTK_RANGE (hscale), content);
  update_pointer (i, content);
  if (isotherm_geomview_pipe || phdiag_geomview_pipe)
    recalculate_and_display ();
}

void on_entry_activate (GtkWidget *widget, gpointer user_data)
{
  int i;

  for (i=0; i<33 && widget != entry_widgets [i]; i++)
    ;
  if (i<33)
    get_entry_update_hscale_and_pointer (i);
  else
    printf ("Internal error, on_entry_activate got widget not in table\n");
}

void on_hscale_value_changed (GtkWidget *widget, gpointer user_data)
{
  int i;

  for (i=0; i<33 && widget != hscale_widgets [i]; i++)
    ;
  if (i<33)
    {
      GtkWidget *entry = entry_widgets [i];
      const gchar entrytext [100];
      gdouble content;

      content = gtk_range_get_value (GTK_RANGE (hscale_widgets [i]));
      snprintf (entrytext, 98, "%g", content);
      gtk_entry_set_text (GTK_ENTRY(entry), entrytext);

      update_pointer (i, content);
      if (isotherm_geomview_pipe || phdiag_geomview_pipe)
        recalculate_and_display ();
    }
  else
    printf ("Internal error, on_hscale_value_changed got widget 0x%lx not in table\n", widget);
}

void on_spinodal_togglebutton_toggled (GtkWidget *widget, gpointer user_data)
{
  spinodal_toggle [0] = gtk_toggle_button_get_active
    (GTK_TOGGLE_BUTTON(glade_xml_get_widget(xml, "P1_spinodal_togglebutton")));
  spinodal_toggle [1] = gtk_toggle_button_get_active
    (GTK_TOGGLE_BUTTON(glade_xml_get_widget(xml, "P2_spinodal_togglebutton")));
  if (isotherm_geomview_pipe || phdiag_geomview_pipe)
    recalculate_and_display ();
}

void on_spinbutton_value_changed (GtkWidget *widget, gpointer user_data)
{
  T_steps = (int) gtk_spin_button_get_value
    (GTK_SPIN_BUTTON (glade_xml_get_widget(xml, "phase_steps_spinbutton")));
  refine = (int) gtk_spin_button_get_value
    (GTK_SPIN_BUTTON (glade_xml_get_widget(xml, "refine_spinbutton")));
  P1_res = (int) gtk_spin_button_get_value
    (GTK_SPIN_BUTTON (glade_xml_get_widget(xml, "P1_res_spinbutton")));
  P2_res = (int) gtk_spin_button_get_value
    (GTK_SPIN_BUTTON (glade_xml_get_widget(xml, "P2_res_spinbutton")));
  numpoints = (P1_res+1)*(P1_res+2)/2 + (P2_res+1)*(P2_res+2)/2;
  recalculate_and_display ();
}

void on_name_entry_activate (GtkWidget *widget, gpointer user_data)
{
  strncpy (allparams[0].name, (char *) gtk_entry_get_text
           (GTK_ENTRY (glade_xml_get_widget (xml, "P1_name_entry"))), 99);
  strncpy (allparams[1].name, (char *) gtk_entry_get_text
           (GTK_ENTRY (glade_xml_get_widget (xml, "P2_name_entry"))), 99);
}

void show_isotherm (GtkWidget *widget, gpointer user_data)
{
  int ierr;
  gboolean active = gtk_toggle_button_get_active (GTK_TOGGLE_BUTTON(widget));

  if (active && isotherm_geomview_pipe == NULL)
    {
      ierr = GeomviewBegin (&isotherm_geomview_pipe, gv_version);
      recalculate_and_display ();
    }
  else if ((!active) && isotherm_geomview_pipe != NULL)
    GeomviewEnd (&isotherm_geomview_pipe);
  else
    printf ("Error: toggle is %d but geomview was 0x%lx\n", active,
            isotherm_geomview_pipe);
}

void save_isotherm_snapshot (GtkWidget *widget, gpointer user_data)
{
  G_CONST_RETURN gchar *entrytext;

  entrytext = gtk_entry_get_text (GTK_ENTRY (widget));

  if (isotherm_geomview_pipe)
    {
      printf ("Saving snapshot %s...\n", entrytext);
      fprintf (isotherm_geomview_pipe, "(snapshot c0 %s)\n", entrytext);
    }
}

void show_phasediag (GtkWidget *widget, gpointer user_data)
{} // Not yet implemented

void save_phase_snapshot (GtkWidget *widget, gpointer user_data)
{} // Not yet implemented

int main (int argc, char *argv[])
{
  int i;
  char filename [200];

  // Initialize the program, and load and connect the Glade XML interface
  gnome_program_init ("Ternary", VERSION, LIBGNOMEUI_MODULE, argc, argv, NULL);
  strncpy (filename, GLADE_DIRECTORY, 186);
  strcat (filename, "/Ternary.glade");
  printf ("Loading Glade XML user interface file %s\n", filename);
  xml = glade_xml_new (filename, NULL, NULL);
  glade_xml_signal_autoconnect (xml);
  // Fill the entry_widgets and hscale_widgets tables
  allparams[0].name = P1_name;
  allparams[1].name = P2_name;
  for (i=32; i>=0; i--)
    {
      entry_widgets [i] = glade_xml_get_widget (xml, entry_names [i]);
      hscale_widgets [i] = glade_xml_get_widget (xml, hscale_names [i]);
      get_entry_update_hscale_and_pointer (i);
    }
  on_spinbutton_value_changed (NULL, NULL);
  on_name_entry_activate (NULL, NULL);
  on_spinodal_togglebutton_toggled (NULL, NULL);

  // Run main loop
  if (i=GeomviewBegin (&isotherm_geomview_pipe, gv_version))
    { printf ("main: Error %d in GeomviewBegin\n", i); exit (i); }
  gtk_widget_show (glade_xml_get_widget (xml, "main_window"));
  recalculate_and_display();
  gtk_main();
  if (isotherm_geomview_pipe)
    GeomviewEnd (&isotherm_geomview_pipe);

  return 0;
}