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

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

  This function contains "bookkeeping" functions, such as setting up the
  initial triangle array of compositions, scaling things to reasonable x,y,z,
  etc.
  ***************************************/


#include "freenergy.h" /*+ Free energy prototypes, typedefs, etc. +*/
#include "gibbs.h" /*+ Ternary prototypes, typedefs, etc. +*/
#include <math.h>    /*+ For sqrt(), fmin()/fmax() +*/
#include <stdlib.h>  /*+ For malloc and free +*/


/*+ Macro returning the start of a row in a triangle array +*/
#define ROWSTART(row) ((row)*(resolution+1) - ((row)*((row)-1))/2)


/*++++++++++++++++++++++++++++++++++++++
  This initializes the compositions
  +latex+$C_2$ and $C_3$
  -latex-C2 and C3
  in an array of points filling a triangle between (A2,A3), (B2,B3) and (C2,C3)
  in composition space.  It also sets the
  +latex+{\tt efunc}, $T$ and $P$
  +html+ <tt>efunc</tt>, <i>T</i> and <i>P</i>
  fields of each point to -1 to indicate that free energies and derivatives are
  not yet calculated.

  int init_triangle_array It returns zero or an error code.

  int resolution Resolution of the discretization (number of intervals on a
  side).

  energy_point *points Array of energy points of size
  +latex+$(N+1)(N+2)/2$, where $N$ is the resolution.
  -latex-(N+1)(N+2)/2, where N is the resolution.

  int efunc Energy function index indicating which energy curve this point will
  be on.

  double A2 First corner C2 value.

  double A3 First corner C3 value.

  double B2 Second corner C2 value.

  double B3 Second corner C3 value.

  double C2 Third corner C2 value.

  double C3 Third corner C3 value.
  ++++++++++++++++++++++++++++++++++++++*/

int init_triangle_array
(int resolution, energy_point *points, int efunc,
 double A2, double A3, double B2, double B3, double C2, double C3)
{
  int i,j, index;

  for (i=0; i<resolution; i++)
    {
      for (j=0; j<resolution-i; j++)
        {
          index = ROWSTART (i) + j;
          points[index].C2 = A2 + (B2-A2) * ((double) i/resolution) +
            (C2-A2) * ((double) j/resolution);
          points[index].C3 = A3 + (B3-A3) * ((double) i/resolution) +
            (C3-A3) * ((double) j/resolution);
          points[index].efunc = efunc;
          points[index].T = points[index].P = -1;
        }
      /* Make sure these are exactly on the C2+C3=1 boundary */
      index = ROWSTART (i) + j;
      points[index].C2 = A2 + (B2-A2) * ((double) i/resolution) +
        (C2-A2) * (1.-((double) i/resolution));
      points[index].C3 = A3 + (B3-A3) * ((double) i/resolution) +
        (C3-A3) * (1.-((double) i/resolution));
    }
  index = (resolution+1)*(resolution+2)/2-1;
  points[index].C2 = B2;
  points[index].C3 = B3;

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  This initializes the compositions
  +latex+$C_2$ and $C_3$
  -latex-C2 and C3
  in an array of points filling a rectangle between (A2,A3), (B2,B3) in
  composition space.  It also sets the
  +latex+{\tt efunc}, $T$ and $P$
  +html+ <tt>efunc</tt>, <i>T</i> and <i>P</i>
  fields of each point to -1 to indicate that free energies and derivatives are
  not yet calculated.  This is of course trivial compared with the triangle
  case.

  int init_rect_array It returns zero or an error code.

  int res2 Resolution of the discretization (number of intervals) in the
  +latex+$C_2$-direction.
  -latex-C2-direction.

  int res3 Resolution of the discretization (number of intervals) in the
  +latex+$C_3$-direction.
  -latex-C3-direction.

  energy_point *points Array of energy points of size
  +latex+({\tt res2}+1)$\times$({\tt res3}+1).
  -latex-(res2+1) * (res3+1).

  int efunc Energy function index indicating which energy curve this point will
  be on.

  double A2 Minimum value of
  +latex+$C_2$.
  -latex-C2.

  double A3 Minimum value of
  +latex+$C_3$.
  -latex-C3.

  double B2 Maximum value of
  +latex+$C_2$.
  -latex-C2.

  double B3 Maximum value of
  +latex+$C_3$.
  -latex-C3.
  ++++++++++++++++++++++++++++++++++++++*/

int init_rectangle_array
(int res2, int res3, energy_point *points, int efunc,
 double A2, double A3, double B2, double B3)
{
  int i,j;

  for (i=0; i<=res2; i++)
    for (j=0; j<=res3; j++)
      {
        points [j*(res2+1)+i].C2 = A2 + (B2-A2) * ((double) i/res2);
        points [j*(res2+1)+i].C3 = A3 + (B3-A3) * ((double) j/res3);
        points [j*(res2+1)+i].efunc = efunc;
        points [j*(res2+1)+i].T = points [j*res2+i].P = -1;
      }

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  This function fills an array with vertex indices.

  int init_triangle_vertices It returns zero or an error code.

  int resolution Resolution of the discretization (number of intervals on a
  side).

  int *vertex_array Array of triangle indices of size
  +latex+$N^2$, where $N$ is the resolution.
  +html+ <i>N</i><sup>2</sup>, where <i>N</i> is the resolution.

  int offset Constant to add to vertex indices, usually 0.
  ++++++++++++++++++++++++++++++++++++++*/

int init_triangle_vertices (int resolution, int *vertex_array, int offset)
{
  int i,j, index;

  for (i=0, index=0; i<resolution-1; i++)
    {
      for (j=0; j<resolution-i-1; j++)
        {
          vertex_array [3*index]   = offset + ROWSTART(i)+j;
          vertex_array [3*index+1] = offset + ROWSTART(i)+j+1;
          vertex_array [3*index+2] = offset + ROWSTART(i+1)+j;
          index++;

          vertex_array [3*index]   = offset + ROWSTART(i)+j+1;
          vertex_array [3*index+1] = offset + ROWSTART(i+1)+j+1;
          vertex_array [3*index+2] = offset + ROWSTART(i+1)+j;
          index++;
        }
      vertex_array [3*index]   = offset + ROWSTART(i)+resolution-i-1;
      vertex_array [3*index+1] = offset + ROWSTART(i)+resolution-i;
      vertex_array [3*index+2] = offset + ROWSTART(i+1)+resolution-i-1;
      index++;
    }
  vertex_array [3*index]   = offset + ROWSTART(resolution-1);
  vertex_array [3*index+1] = offset + ROWSTART(resolution-1)+1;
  vertex_array [3*index+2] = offset + ROWSTART(resolution);

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  This function fills an array with triangle vertex indices, two triangles per
  rectangular cell.

  int init_rectangle_vertices It returns zero or an error code.

  int res2 Resolution of the discretization (number of intervals) in the
  +latex+$C_2$-direction.
  -latex-C2-direction.

  int res3 Resolution of the discretization (number of intervals) in the
  +latex+$C_3$-direction.
  -latex-C3-direction.

  int *vertex_array Array of triangle indices of size
  +latex+2$\times${\tt res2$\times$res3}.
  -latex-2 * res2 * res3.

  int offset Constant to add to vertex indices, usually 0.
  ++++++++++++++++++++++++++++++++++++++*/

int init_rectangle_vertices (int res2, int res3, int *vertex_array, int offset)
{
  int i,j;

  for (i=0; i<res2; i++)
    for (j=0; j<res3; j++)
      {
        vertex_array [(j*res2+i)*6]   = offset + j*(res2+1)+i;
        vertex_array [(j*res2+i)*6+1] = offset + j*(res2+1)+i+1;
        vertex_array [(j*res2+i)*6+2] = offset + (j+1)*(res2+1)+i;

        vertex_array [(j*res2+i)*6+3] = offset + j*(res2+1)+i+1;
        vertex_array [(j*res2+i)*6+4] = offset + (j+1)*(res2+1)+i+1;
        vertex_array [(j*res2+i)*6+5] = offset + (j+1)*(res2+1)+i;
      }

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  This scales the free energy of an array of ternary points to the
  +latex+$y$-coordinate
  -latex-y-coordinate
  and colors the points according to free energy as well.

  int energy_to_visual_array It returns zero or an error code.

  int num_points Number of points in the array

  energy_point *epoints Array of energy points of size
  +latex+$(N+1)(N+2)/2$, where $N$ is the resolution.
  -latex-(N+1)(N+2)/2, where N is the resolution.

  visual_point **tpoints Pointer to array of ternary points of size
  +latex+$(N+1)(N+2)/2$, where $N$ is the resolution.
  -latex-(N+1)(N+2)/2, where N is the resolution.
  This is changed by realloc() each time the function is called, so it should
  be NULL or a valid array.

  double y0 Reference (minimum)
  +latex+$y$-coordinate value.
  -latex-y-coordinate value.

  double G0 Free energy corresponding to reference
  +latex+$y$-coordinate value.
  -latex-y-coordinate value.

  double y1 Second (maximum)
  +latex+$y$-coordinate value.
  -latex-y-coordinate value.

  double G1 Free energy corresponding to second
  +latex+$y$-coordinate value.
  -latex-y-coordinate value.

  double *Gmin Optional address to store the minimum free energy in the array
  on return (NULL if unwanted), used as G0 if both G0 and G1 are zero.

  double *Gmax Optional address to store the maximum free energy in the array
  on return (NULL if unwanted), used as G1 if both G0 and G1 are zero.

  int square Flag: zero for triangle composition representation, non-zero for
  square.
  ++++++++++++++++++++++++++++++++++++++*/

int energy_to_visual_array
(int num_points, energy_point *epoints, visual_point **tpoints, double y0,
 double G0, double y1, double G1, double *Gmin, double *Gmax, int square)
{
  int i;

  if (!(*tpoints = (visual_point *) realloc
        (*tpoints, num_points * sizeof (visual_point))))
    return -1;

  if (G0==0. && G1==0.)
    {
      G0 = G1 = epoints[0].G;
      for (i=0; i<num_points; i++)
        {
          G0 = fmin (epoints[i].G, G0);
          G1 = fmax (epoints[i].G, G1);
        }
      if (Gmin)
        *Gmin = G0;
      if (Gmax)
        *Gmax = G1;
    }
  else
    {
      if (Gmin)
        {
          *Gmin = epoints[0].G;
          for (i=0; i<num_points; i++)
            *Gmin = fmin (epoints[i].G, *Gmin);
        }
      if (Gmax)
        {
          *Gmax = epoints[0].G;
          for (i=0; i<num_points; i++)
            *Gmax = fmax (epoints[i].G, *Gmax);
        }
    }

#define Grel(G) (y0 + (y1-y0) * ((G)-G0) / (G1-G0))
#define RED(G) (((G)<.25) ? 1. : (((G)<.5) ? 2.-4.*(G) : 0.))
#define GREEN(G) (((G)<.25) ? 4.*(G) : (((G)<.75) ? 1. : 4.-4.*(G)))
#define BLUE(G) (((G)<.5) ? 0. : (((G)<.75) ? 4.*(G)-2. : 1.))

  for (i=0; i<num_points; i++)
    {
      if (square)
        {
          (*tpoints)[i].x = epoints[i].C2;
          (*tpoints)[i].z = 1.-epoints[i].C3;
        }
      else
        {
          (*tpoints)[i].x = epoints[i].C2 + 0.5*epoints[i].C3;
          (*tpoints)[i].z = (1. - epoints[i].C3) * sqrt(3)/2.;
        }
      (*tpoints)[i].y = Grel(epoints[i].G);
      (*tpoints)[i].red   = RED   (Grel (epoints[i].G));
      (*tpoints)[i].green = GREEN (Grel (epoints[i].G));
      (*tpoints)[i].blue  = BLUE  (Grel (epoints[i].G));
      (*tpoints)[i].alpha = 1.;
    }

  return 0;
}