xmd Code

/*
xmd - molecular dynamics for metals and ceramics
By Jonathan Rifkin <jon.rifkin@uconn.edu>
Copyright 1995-2006 Jonathan Rifkin

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
/*
Set default vi settings:
vi: ts=3
*/

/*
************************************************************************
History
************************************************************************
*/

/*
Jan 14, 1999
  Have found that plot_coord() does not handle free surface
  case correctly when particles positions are less than
  zero or greater then a->bcur[].

  A quick hack is to have ClipAtom() ignore cases when
  both atoms are outside the box.

Mar 1, 2006
  - Add PLOT SYMBOL FILLALL option to plot filled atoms *without* the 
     black perimeter.
  - Add the PLOT LABEL {ON|OFF} command.
*/

/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "sortsub.h"
#include "strsub.h"
#include "parse.h"
#include "cdhouse.h"
#include "particle.h"
#include "grdevsub.h"
#include "plotc.h"
#include "cdsubs.h"
#include "cdboun.h"

/*
************************************************************************
Defines
************************************************************************
*/
#define  MAXUNS   0x8fff
#define  RADIUS     0.032
#define  MTOP       1.0
#define  MBOT       1.5
#define  MRIGHT     0.80
#define  MLEFT      1.00
#define  MARGIN     0.20
#define  HLETT      0.125

/*  Default page size (in inches)  */
#define  XPAGE      8.5
#define  YPAGE     11.0

#define FILL_OFF 0
#define FILL_ON  1
#define FILL_ALL 2



/*
************************************************************************
Type Definitions
************************************************************************
*/
typedef struct
   {
   int iflag, npl;
   int page, nlayerplot, nlayerbox;
   char *symbol;
   char *fill; 
   ColorTable_t *ColorTable;
   ColorModel_t *ColorModel;
   int *layerplot;
   BOOLEAN UseDisplacement;
   BOOLEAN UseBond;
   BOOLEAN WriteLabel;
   float bondlo, bondhi, dispscale;
   int orient[3];
   float *radius;
   int DEVflag;
   float PageSize[2];
   double MinCoord[NDIR];
   } PlotInfo_t;

enum symbols { NON, CIR, CRO, TRI, ITR, DIA, AST, SQU, NUM };


/*
************************************************************************
External Variables
************************************************************************
*/
extern Particle_t *a;



/*
************************************************************************
Module-Wide Variables
************************************************************************
*/
static PlotInfo_t *PlotInfo_m     = NULL;
static char         *OrientSymbol_m = "XYZ";
static ColorTable_t  BlackTable_m   = {0,0,0};
static ColorTable_t  WhiteTable_m   = {1,1,1};



/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void ReadPlot_Bond   (char *);
void ReadPlot_Clear  (char *);
void ReadPlot_Device (char *);
void ReadPlot_Disp   (char *);
void ReadPlot_Label  (char *);
void ReadPlot_Layer  (char *);
void ReadPlot_Orient (char *);
void ReadPlot_Page   (char *);
void ReadPlot_Size   (char *);
void ReadPlot_Symbol (char *);
void ReadPlot_Write  (char *);

void plot_coord (char *outname, Particle_t *a, PlotInfo_t *pl, int UseSelect);

PlotInfo_t *InitializePlotInfo (void);

void CalcDisp (float Disp[NDIR], Particle_t *a, int i);
void ClipAtom (double [NDIR], float [NDIR], double [NDIR], int, int);

void DrawCircle (int ipart, float radius, PlotInfo_t *p);

/*
************************************************************************
Exported Functions
************************************************************************
*/

void read_plot  (char *InputString)
   {
   char *LeadingToken;

   /*  Initialize space for plot information  */
   if (PlotInfo_m==NULL)
      PlotInfo_m = InitializePlotInfo();

   /*  Get leading token  */
   LeadingToken = strhed (&InputString);

   /*  Bond command  */
   if      (!strcmpi (LeadingToken, "BOND"  ))
      ReadPlot_Bond   (InputString);

   /*  Device command  */
   else if (!strcmpi (LeadingToken, "DEVICE"))
      ReadPlot_Device (InputString);

   /*  Displacement command  */
   else if (!strcmpi (LeadingToken, "DISP"  ))
      ReadPlot_Disp   (InputString);

   /*  Label command  */
   else if (!strcmpi (LeadingToken, "LABEL" ))
      ReadPlot_Label (InputString);

   /*  Layer command  */
   else if (!strcmpi (LeadingToken, "LAYER" ))
      ReadPlot_Layer  (InputString);

   /*  Orient command  */
   else if (!strcmpi (LeadingToken, "ORIENT"))
      ReadPlot_Orient (InputString);

   /*  Page command  */
   else if (!strcmpi (LeadingToken, "PAGE"  ))
      ReadPlot_Page (InputString);

   /*  Symbol command  */
   else if (!strcmpi (LeadingToken, "SIZE"))
      ReadPlot_Size (InputString);

   /*  Symbol command  */
   else if (!strcmpi (LeadingToken, "SYMBOL"))
      ReadPlot_Symbol (InputString);

   /*  Write command  */
   else if (!strcmpi (LeadingToken, "WRITE" ))
      ReadPlot_Write (InputString);

   /*  else unknown command  */
   else
      {
      ERROR_PREFIX
      printf ("Unknown plot command.\n");
      CleanAfterError();
      }
   }

void ReadPlot_Bond  (char *instr)
   {
   char *tokstr = strhed (&instr);
   if (!strcmpi(tokstr,"ON"))
      {
      PlotInfo_m->UseBond = TRUE;
      PlotInfo_m->bondlo = dblstrf (&instr) * 1e-8;
      PlotInfo_m->bondhi = dblstrf (&instr) * 1e-8;
      }
   else if (!strcmpi(tokstr,"OFF"))
      PlotInfo_m->UseBond = FALSE;
   else
      {
      ERROR_PREFIX
      printf ("Unknown bond option (%s)\n", instr);
      CleanAfterError();
      }
   }


void ReadPlot_Disp  (char *instr)
   {
   char *tokstr = strhed (&instr);
   if      (!strcmpi(tokstr,"ON"))
      PlotInfo_m->UseDisplacement = TRUE;
   else if (!strcmpi(tokstr,"OFF"))
      PlotInfo_m->UseDisplacement = FALSE;
   else if (!strcmpi(tokstr,"SCALE"))
      PlotInfo_m->dispscale = dblstrf (&instr);
   else
      {
      ERROR_PREFIX
      printf ("Unknown option (%s)\n", tokstr);
      CleanAfterError();
      }
   }

void ReadPlot_Device(char *instr)
   {
   if     (!strcmpi(instr,"CANON"))
      PlotInfo_m->DEVflag = 1;
   else if (!strcmpi(instr,"POSTSCRIPT"))
      PlotInfo_m->DEVflag = 0;
   else
      {
      ERROR_PREFIX
      printf ("Unknown device option (%s)\n", instr);
      CleanAfterError();
      }
   }

void ReadPlot_Label (char *instr) {
   char *tokstr = strhed (&instr);
   if (!strcmpi(tokstr,"ON"))
      PlotInfo_m->WriteLabel = TRUE;
   else if (!strcmpi(tokstr,"OFF"))
      PlotInfo_m->WriteLabel = FALSE;
   else {
      ERROR_PREFIX
      printf ("Unknown label option (%s)\n", instr);
      CleanAfterError();
   }
}

void ReadPlot_Layer (char *instr)
   {
   int t;

   /*  READ NUMBER OF LAYERS  */
   PlotInfo_m->nlayerbox = intstrf (&instr);

   /*  TEST FOR ZERO LAYERS  */
   if (PlotInfo_m->nlayerbox<=0)
      {
      ERROR_PREFIX
      printf 
         (
         " Number of layers (%i) must be greater than zero.\n",
         PlotInfo_m->nlayerbox
         );
      CleanAfterError();
      }

   /*  ALLOCATE SPACE  */
   FREE  (PlotInfo_m->layerplot)
   ALLOCATE (PlotInfo_m->layerplot, int, PlotInfo_m->nlayerbox)
   if (*instr)
      {
      PlotInfo_m->nlayerplot = 0;
      while (*instr && PlotInfo_m->nlayerplot<PlotInfo_m->nlayerbox)
         {
         t = intstrf (&instr);
         if (t>PlotInfo_m->nlayerbox)
            {
            ERROR_PREFIX
            printf ("Invalid layer (%i of %i)\n",
               t, PlotInfo_m->nlayerbox);
            CleanAfterError();
            }
         else
            PlotInfo_m->layerplot[PlotInfo_m->nlayerplot++] = t;
         }
      }
   else
      {
      PlotInfo_m->nlayerplot = PlotInfo_m->nlayerbox;
      LOOP(t,  PlotInfo_m->nlayerbox)
         PlotInfo_m->layerplot[t] = t+1;
      }
   }

void ReadPlot_Orient(char *instr)
   {
   char t;
   int  orient[3];
   /*  CHECK FORMAT (must be "X/Z" or "Y/X" etc  */
   while (*instr && *instr==' ')
      instr++;
   t = instr[0];
   if (t>='x' && t<='z')
      t += 'X'-'x';
   orient[1] = t - 'X';

   t = instr[2];
   if (t>='x' && t<='z')
      t += 'X'-'x';
   orient[0] = t - 'X';

   orient[2] = 3 - orient[1] - orient[0];

   if (   orient[0]<0 || orient[1]<0 || orient[2]<0
       || orient[0]>2 || orient[1]>2 || orient[2]>2)
      {
      ERROR_PREFIX
      printf ("Invalid ORIENT specification (%s)\n",
    instr);
      CleanAfterError();
      }
   else
      {
      PlotInfo_m->orient[0] = orient[0];
      PlotInfo_m->orient[1] = orient[1];
      PlotInfo_m->orient[2] = orient[2];
      }
   }

void ReadPlot_Page  (char *instr)
   {
   PlotInfo_m->page = intstrf (&instr);
   }

void ReadPlot_Size (char *instr) {
   PlotInfo_m->PageSize[X] = dblstrf(&instr);
   PlotInfo_m->PageSize[Y] = dblstrf(&instr);
}


void ReadPlot_Symbol(char *instr) {
   ColorModel_t ColorModel;
   ColorTable_t ColorTable = {0, 0, 0, 0};
   int ipart;
   int icolor;
   int NumColor;
   float Radius;
   float CurrentColor;
   char  Symbol;
   char *TokenStr = NULL;
   BOOLEAN FillFlag = FILL_OFF;
   BOOLEAN ReadColorFlag;
   

   /*  Make sure some atoms are selected  */
   CheckForNoneSelected();

   TokenStr = strhed  (&instr);
   
   FillFlag = FILL_OFF;

	/*  Set atom fill with black outline  */
   if (!strcmpi(TokenStr,"FILL")) {
      FillFlag  = FILL_ON;
      TokenStr = strhed  (&instr);
		
	/*  Set atom fill with no outline  */
   } else if (!strcmpi(TokenStr,"FILLALL")) {
      FillFlag  = FILL_ON;
      TokenStr = strhed  (&instr);
   }

   /*  Parse color  */
   ReadColorFlag = TRUE;
   if (!strcmpi(TokenStr,"RGB")) {
      ColorModel = COLOR_RGB;
      NumColor  = 3;
   } else if (!strcmpi(TokenStr, "HSB")) {
      ColorModel = COLOR_HSB;
      NumColor  = 3;
#ifdef INCLUDE_CMYB_MODEL
	} else if (!strcmpi(TokenStr, "CMYB")) {
      ColorModel = COLOR_CMYB;
      NumColor  = 4;
#endif
	} else {
      ReadColorFlag = FALSE;
      ColorModel    = COLOR_RGB;
      NumColor      = 3;
      ColorTable[0] = 0;
      ColorTable[1] = 0;
      ColorTable[2] = 0;
   }

   /*  
   If  color command found, 
   read in appropriate number of color paramters  
   */
   if (ReadColorFlag) {
      LOOP (icolor, NumColor) {
         CurrentColor = dblstrf(&instr);
         ColorTable[icolor] = 
            CurrentColor < 0.0 ?   0 :
            CurrentColor > 1.0 ? 255 :
            CurrentColor  * 255;
      }
      TokenStr  = strhed (&instr);
   }


   /*  Convert  radius from  angs to cm  */
   Radius    = dblstrf (&instr) * 1E-8;  

   /*  SET SYMBOL TYPE  */
   if      (!strcmpi(TokenStr, "NON"))
      Symbol = NON;
   else if (!strcmpi(TokenStr, "CIR"))
      Symbol = CIR;
   else if (!strcmpi(TokenStr, "CRO"))
      Symbol = CRO;
   else if (!strcmpi(TokenStr, "TRI"))
      Symbol = TRI;
   else if (!strcmpi(TokenStr, "ITR"))
      Symbol = ITR;
   else if (!strcmpi(TokenStr, "DIA"))
      Symbol = DIA;
   else if (!strcmpi(TokenStr, "AST"))
      Symbol = AST;
   else if (!strcmpi(TokenStr, "SQU"))
      Symbol = SQU;
   else if (!strcmpi(TokenStr, "NUM"))
      Symbol = NUM;
   else {
      ERROR_PREFIX
      printf ("Symbol type (%s) not recognized.\n", TokenStr);
      CleanAfterError();
   }

   /*  Allocate storage  */
   if (PlotInfo_m->npl!=a->np) {
      FREE  (PlotInfo_m->symbol)
      FREE  (PlotInfo_m->radius)
      FREE  (PlotInfo_m->fill  )
      FREE  (PlotInfo_m->ColorTable)
      FREE  (PlotInfo_m->ColorModel)
      ALLOCATE (PlotInfo_m->symbol, char,         a->np)
      ALLOCATE (PlotInfo_m->radius, float,        a->np)
      ALLOCATE (PlotInfo_m->fill,   char,         a->np)
      ALLOCATE (PlotInfo_m->ColorTable,  ColorTable_t,  a->np)
      ALLOCATE (PlotInfo_m->ColorModel,  ColorModel_t,  a->np)
      PlotInfo_m->npl    = a->np;
   }

   /*  Store plot information  */
   LOOP (ipart, a->np)
      if (IS_SELECT(ipart)) {
         PlotInfo_m->symbol[ipart] = Symbol;
         PlotInfo_m->radius[ipart] = Radius;
         PlotInfo_m->fill  [ipart] = FillFlag;
         PlotInfo_m->ColorModel[ipart] = ColorModel;
         LOOP (icolor, NumColor) {
            PlotInfo_m->ColorTable[ipart][icolor] = ColorTable[icolor];
         }
      }
}


void ReadPlot_Write  (char *InputString) {
   char *LeadingToken;
   BOOLEAN UseSelect;

   LeadingToken = strhed (&InputString);
   UseSelect = !strcmpi (LeadingToken, "SEL");

   if (UseSelect)
      LeadingToken = strhed (&InputString);

   plot_coord (LeadingToken, a, PlotInfo_m, UseSelect);
}

/*
************************************************************************
Local Functions
************************************************************************
*/

PlotInfo_t *InitializePlotInfo (void)
   {
   PlotInfo_t *t = NULL;

   ALLOCATE(t, PlotInfo_t, 1)
   t->orient[0] = 0;
   t->orient[1] = 1;
   t->orient[2] = 2;
   t->page      = 1;
   t->nlayerplot= 1;
   t->nlayerbox = 1;
   ALLOCATE(t->layerplot, int, 1)
   t->layerplot[0] = 1;
   t->npl = 0;
   t->dispscale = 1;
   t->PageSize[X] = XPAGE;
   t->PageSize[Y] = YPAGE;
   t->UseDisplacement = FALSE;
   t->UseBond         = FALSE;
   t->WriteLabel      = TRUE;
   return t ;
   }




/*  Draw a single particle - called from plot_coord()  */
void draw_part
(
Particle_t *a,
int i,
PlotInfo_t *pl,
double xstart,
double ystart,
double scale,
int hor,
int ver
)
   {
   float xpos, ypos, srad, textsize;
   char tstr[64];
   float Disp [NDIR];
   float Image[NDIR];

   if (pl->symbol[i] || pl->UseDisplacement)
      {
      xpos = xstart+scale*(a->cur[hor + 3*i]-pl->MinCoord[hor]);
      ypos = ystart+scale*(a->cur[ver + 3*i]-pl->MinCoord[ver]);
      srad =        scale*pl->radius[i];
      grdev_moveto (xpos, ypos);
      grdev_setfill (pl->fill[i]);
      grdev_setcolor (pl->ColorModel[i], pl->ColorTable[i]);
      switch (pl->symbol[i])
         {
         case CIR : DrawCircle (i, srad, pl); break;
#if 0
         case CIR : grdev_circle     (srad); break;
#endif
         case CRO : grdev_cross      (srad); break;
         case TRI : grdev_triangle   (srad); break;
         case ITR : grdev_itriangle  (srad); break;
         case DIA : grdev_diamond    (srad); break;
         case AST : grdev_asterisk   (srad); break;
         case SQU : grdev_square     (srad); break;
         case NUM : textsize = grdev_gettextsize();
                    grdev_settextsize (srad);
                    sprintf (tstr, "%i", i+1);
                    grdev_writetext  (tstr);
                    grdev_settextsize (textsize);
                    break;
         }
      if (pl->UseDisplacement)
         {
         /*  Find displacment (include wrapping)  */
         CalcDisp (Disp, a, i);

         /*  Scale displacement  */
         Disp[hor] *= pl->dispscale;
         Disp[ver] *= pl->dispscale;

         /*  Calculate image position  */
         Image[hor] = a->cur[hor + 3*i] + Disp[hor];
         Image[ver] = a->cur[ver + 3*i] + Disp[ver];

         /*  Clip image position if repeating boundary  */
         ClipAtom (&a->cur[3*i], Image, a->bcur, hor, ver);

         /*  Draw tail to image position  */
         xpos = xstart + scale*Image[hor];
         ypos = ystart + scale*Image[ver];
         grdev_drawto (xpos, ypos);
         }
      }
   }


/*  Draw all particles and layers for a given step  */
void plot_coord (char *outname, Particle_t *a, PlotInfo_t *pl, int UseSelect) {
   time_t tstruct;

   FILE  *fout;
   int   *layerbot = NULL;
   int   *layertop = NULL;
   unsigned *ix = NULL;
   int      *vect = NULL;
   float xlimit, ylimit, ximage, yimage, xframe, yframe;
   float xgap, ygap , xorg, yorg;
   float xstart, ystart, scale, ratio;
   float tempbox, t;
   int   ipart, jpart, ilayerbox, ilayerplot;
   int   nx, ny, nbox;
   int   curpagelayer, npagelayer, curpage;
   int   iline, iind, jind;

   int    hor  = pl->orient[0];
   int    ver  = pl->orient[1];
   int    dep  = pl->orient[2];
   double MaxCoord[NDIR];
   double (*c)[3] = (double (*)[3]) a->cur;

   /*
   TEST DATA INTEGRITY
   */

   /*  TEST BOX  */
   if (a->bcur[0]<=0 || a->bcur[1]<=0 || a->bcur[2]<=0) {
      ERROR_PREFIX
      printf ("Cannot plot, box dimensions are less or equal to zero.\n");
      CleanAfterError();
   }

   /*  TEST FOR INITIALIZE OF PlotInfo_t  */
   if (pl->npl<a->np) {
      ERROR_PREFIX
      printf ("Cannot plot, not symbols specified.\n");
      CleanAfterError();
   }

   /*  WHEN PLOTTING DISPLACEMENTS, TEST FOR REFERENCE PARTICLES  */
   if (pl->UseDisplacement && a->ref==NULL) {
         ERROR_PREFIX
         printf ("Cannot plot displacements\n");
         printf ("   because no reference particle have been specified.\n");
         CleanAfterError();
      }

   /*  Wrap particles  */
   WrapParticles (a);

   /* set directions on plotted page */

   /*  room on page for drawing - in inches  */
   xlimit = pl->PageSize[X] - MLEFT - MRIGHT;  
   ylimit = pl->PageSize[Y] - MBOT  - MTOP  ;

	/*  Get min,max coordinates  (use Box if no free surface)  */
	GetMinMaxCoord (a, pl->MinCoord, MaxCoord, UseSelect, TRUE);

   /* size of image in simulation units   */
	ximage = MaxCoord[hor] - pl->MinCoord[hor];
	yimage = MaxCoord[ver] - pl->MinCoord[ver];

   /* Calculate maximum number of layers drawn on page */
   npagelayer = (pl->nlayerplot-1) / pl->page + 1;

   /*
     If all layers can't fit in one column then rescale xFrame and yFrame
     so that they fit into the least number of multiple columns
   */
   if (npagelayer==1)
      nx = ny = 1;
   else {
      ratio = (ylimit/yimage) / (xlimit/ximage);
      if (ratio>=1.0) {
          ny = sqrt (npagelayer * ratio) + 1;
          nx = (npagelayer-1) / ny + 1;
      } else {
         nx = sqrt (npagelayer / ratio) + 1;
         ny = (npagelayer-1) / nx + 1;
      }
   }


   /* conversion factor s.u. -> in.*/
   t      = (xlimit-2*nx*MARGIN)/ximage/nx;
   scale  = (ylimit-2*ny*MARGIN)/yimage/ny;
   if (t<scale)
      scale=t;
   xframe = scale*ximage + 2*MARGIN;  /* size of frame around float layer  */
   yframe = scale*yimage + 2*MARGIN;
   /* gap between frame and allotted space (inches) */
   xgap   = (xlimit/nx - xframe);
   ygap   = (ylimit/ny - yframe);
   /*x, y  coordiantes of first layer drawing origin (inches) */
   xorg   = MLEFT + xgap/2.0;
   yorg   = MBOT  + ygap/2.0;

	/* order atoms by layer */

   ALLOCATE(ix,   unsigned int, a->np)
   ALLOCATE(vect, int, a->np)
   tempbox = a->bcur[dep];
   if (tempbox==0) {
      printf ("INTERNAL ERROR (plot_coord): Box in depth direction is zero.\n");
      CleanAfterError();
   }

   /*  Wrap particles in depth direction  */
   LOOP (ipart, a->np) {
      t = fabs(c[ipart][dep])/tempbox;
      /*  Subtract integer number of boxes depths from depth coordinate  */
      if (t>=1) t -= (long) t;
      if (c[ipart][dep]<0)   /* RECOVER SIGN (do we really want this?) */ 
         t  = -t;
      c   [ipart][dep]  = tempbox * t;
      /*  Convert depth position to integer for ordering purposes  */
      vect[ipart]       = MAXUNS * t;
      ix  [ipart]       = ipart;
   }


   sortix (vect, ix, a->np);
   FREE(vect)


   /*  ASSIGN SPACE FOR LAYER TOP AND BOX INDICES  */
   ALLOCATE(layertop, int, pl->nlayerbox)
   ALLOCATE(layerbot, int, pl->nlayerbox)

   LOOP (ilayerbox,  pl->nlayerbox) {
      layertop[ilayerbox] = 1;
      layerbot[ilayerbox] = 0;
   }
   tempbox = a->bcur[dep];
   jind    = a->np-1;
   LOOP (iind, a->np) {
      ipart  = ix[iind];
      jpart  = ix[jind];
      ilayerbox =    pl->nlayerbox*c[ipart][dep]/tempbox;
      layerbot[ilayerbox] = iind;
      ilayerbox =    pl->nlayerbox*c[jpart][dep]/tempbox;
      layertop[ilayerbox] = jind;
      jind--;
   }


   /*  OPEN OUTPUT FILE  */
   if (*outname=='+')
      fout = fopen (outname+1, "at");
   else
      fout = fopen (outname  , "wt");

   /*  TEST FOR FILE OK  */
   if (fout==NULL) {
      ERROR_PREFIX
      printf ("Cannot open output file\n");
      CleanAfterError();
   }

   /*  SET DEVICE  */
   if (pl->DEVflag==1)
      grdev_canon();
   else if (pl->DEVflag==0)
      grdev_ps();

   /*  ASSIGN FILE TO PLOTTER  */
   grdev_assignfile (fout);

   /*  INITIALIZE PAGE  */
   grdev_init ();
   grdev_initfile();
   grdev_settextsize(0.11111111);  /*  1/9  */

          /*  loop over layers  */

   nbox = nx * ny;
   LOOP (ilayerplot, pl->nlayerplot) {
      ilayerbox = pl->layerplot[ilayerplot]-1;
      /*   determine current page */
      curpage      = ilayerplot / npagelayer;
      /*   determine position of this layer on current page */
      curpagelayer = ilayerplot % npagelayer;

      /*   eject current page if starting new page  */
      if ( (curpagelayer==0) && (curpage>0)  )
         grdev_writepage();

      /*   plot  identification   */
      if (pl->WriteLabel) {
	      if (curpagelayer==0) {
		 grdev_moveto (MLEFT, MBOT);
		 LOOP (iline, 8)
		 if (a->title[iline] && *(a->title[iline]) )
		    {
		    grdev_moveto (MLEFT, MBOT-(iline+1)*HLETT);
		    grdev_writetext (a->title[iline]);
		 }

		 grdev_moveto (MLEFT, MBOT-9*HLETT);
		 if (a->run>=0)
		    grdev_printf    ("RUN  %6li", a->run);
		 grdev_moverel (1.25, 0);
		 if (a->step>=0)
		    grdev_printf    ("STEP %6li", a->step);
		 grdev_moverel (1.25, 0);
		 grdev_printf
		    ("ORIENT %c/%c", OrientSymbol_m[ver], OrientSymbol_m[hor]);
		 grdev_moverel (1.25, 0);
		 time(&tstruct);
		 grdev_writetext (ctime(&tstruct));
		 grdev_moverel (1.75, 0);
		 grdev_printf    ("Page %i of %i", curpage+1, pl->page);
	      }
      }

      if (pl->WriteLabel) {
	      xstart = xorg +  (     curpagelayer  ) / ny * (xframe + xgap);
	      ystart = yorg +  (nbox-curpagelayer-1) % ny * (yframe + ygap);
	      /* draw frame around current layer */
	      grdev_moveto ( xstart,  ystart);
	      grdev_drawrel( xframe,  0     );
	      grdev_drawrel( 0,       yframe);
	      grdev_drawrel(-xframe,  0     );
	      grdev_drawrel( 0,      -yframe);
	      /* write layer number */
	      grdev_moveto (xstart+0.25*MARGIN, ystart+0.25*MARGIN);
	      grdev_printf ("%i of %i", ilayerbox+1, pl->nlayerbox);
      } else {
	      xstart = xorg +  (     curpagelayer  ) / ny * (xframe + xgap);
	      ystart = yorg +  (nbox-curpagelayer-1) % ny * (yframe + ygap);
      }

      xstart = xstart + MARGIN;
      ystart = ystart + MARGIN;

          /*  loop over particles in current layer  */

      for (iind=layertop[ilayerbox];iind<=layerbot[ilayerbox];iind++) {
         ipart= ix[iind];
         if ( !UseSelect || IS_SELECT(ipart))
            draw_part (a, ipart, pl, xstart, ystart, scale, hor, ver);
      }

      /*  Reset color to black after drawing particles  */
      grdev_setcolor (COLOR_RGB, BlackTable_m);

      /*  DRAW BONDS IN CURRENT LAYER */
      if (pl->UseBond) {
         float rsq,dt,dx,dy,dz, x1,y1,z1;
         float bond1,bond2,bond3,blosq,bhisq, xpos,ypos;

         bond1 = pl->bondhi;
         bond2 = 1.414215*bond1;  /*  sqrt(2)  */
         bond3 = 1.732052*bond1;  /*  sqrt(3)  */
         blosq = pl->bondlo * pl->bondlo;
         bhisq = pl->bondhi * pl->bondhi;
         for (iind=layertop[ilayerbox];iind<=layerbot[ilayerbox];iind++) {
            ipart = ix[iind];
            if ( !UseSelect || IS_SELECT(ipart)) {
               x1 = c[ipart][hor];
               y1 = c[ipart][ver];
               z1 = c[ipart][dep];
               for (jind=iind+1;jind<=layerbot[ilayerbox];jind++) {
                  jpart = ix[jind];
                  if ( !UseSelect || IS_SELECT(jpart)) {
                     dx = fabs(x1-c[jpart][hor]);
                     if (dx<bond1) {
                        dy = fabs(y1-c[jpart][ver]);
                        dt = dx+dy;
                        if (dt<bond2) {
                           dz = fabs(z1-c[jpart][dep]);
                           if (dt+dz<bond3) {
                              rsq = dx*dx + dy*dy + dz*dz;
                              /*  DRAW CURRENT BOND  */
                              if (rsq>=blosq && rsq<=bhisq) {
                                 xpos = xstart+scale*c[ipart][hor];
                                 ypos = ystart+scale*c[ipart][ver];
                                 grdev_moveto (xpos, ypos);
                                 xpos = xstart+scale*c[jpart][hor];
                                 ypos = ystart+scale*c[jpart][ver];
                                 grdev_drawto (xpos, ypos);
                              }
                           }
                        }
                     }
                  }
               }
            }
         }
      }
   }

   /*  PRINT FINAL PAGE ??  */
   grdev_writepage();

   /*  CLOSE FILE  */
   fclose (fout);

   /*  FREE STORAGE  */
   FREE(ix)
   FREE(layerbot)
   FREE(layertop)
}


void CalcDisp (float Disp[NDIR], Particle_t *a, int i)
   {
   int idir;

   LOOP (idir, NDIR)
      {
      Disp[idir] = a->ref[idir + 3*i] - a->cur[idir + 3*i];
      if (!a->surf[idir])
         {
         if      (2*Disp[idir] < -a->bcur[idir])
            Disp[idir] += a->bcur[idir];
         else if (2*Disp[idir] >  a->bcur[idir])
            Disp[idir] -= a->bcur[idir];
         }
      }
   }


/*
Assume position 1 is inside box which runs from 0 to Box[]
Move   position 2 inside box along line joining 1 and 2
*/
void ClipAtom
(
double  Pos1 [NDIR],
float   Pos2 [NDIR],
double  Box [NDIR],
int     hor,
int     ver
)
   {
   double hDel;
   double vDel;
   double Slope;

   /*  
   If both atoms are outside the box, don't clip  
   This is a quick-dirty hack to correct for the
   fact that plot_coord() does not handle free
   surface simulations well when atoms are outside
   the nominal box.
   */
   if 
   (
      (
      Pos1[hor] < 0 || Pos1[hor]>=Box[hor] ||
      Pos1[ver] < 0 || Pos1[ver]>=Box[ver]
      ) 
      &&
      (
      Pos2[hor] < 0 || Pos2[hor]>=Box[hor] ||
      Pos2[ver] < 0 || Pos2[ver]>=Box[ver]
      )
   )
      return;
   

   /*  Displacement between two positions  */
   hDel = Pos2[hor] - Pos1[hor];
   vDel = Pos2[ver] - Pos1[ver];

   /*  Two points are coincident - no need to clip  */
   if (hDel==0 && vDel==0)
      return;

   /*  Horizontal Line  */
   if (vDel==0.0)
      {
      if (Pos2[hor] < 0.0)
         Pos2[hor] = 0.0;
      else if (Pos2[hor] > Box[hor])
         Pos2[hor] = Box[hor];
      return ;
      }

   /*  Vertical Line  */
   if (hDel==0)
      {
      if (Pos2[ver] < 0.0)
         Pos2[ver] = 0.0;
      else if (Pos2[ver] > Box[ver])
         Pos2[ver] = Box[ver];
      return;
      }


   /*  General case  */
   Slope = vDel / hDel;

   if (Pos2[hor] < 0.0)
      {
      Pos2[ver] -= Slope * Pos2[hor];
      Pos2[hor]  = 0.0;
      }
   else if (Pos2[hor] > Box[hor])
      {
      Pos2[ver] += Slope * (Box[hor] - Pos2[hor]);
      Pos2[hor]  = Box[hor];
      }

   if (Pos2[ver] < 0.0)
      {
      Pos2[hor] -= Pos2[ver] / Slope;
      Pos2[ver]  = 0.0;
      }
   else if (Pos2[ver] > Box[ver])
      {
      Pos2[hor] += (Box[ver] - Pos2[ver]) / Slope;
      Pos2[ver]  = Box[ver];
      }

   }


void DrawCircle (int ipart, float radius, PlotInfo_t *p) {

   /*  Outline only  */
   if (p->fill[ipart]==FILL_OFF) {
		/*  Draw outline  */
      grdev_setfill  (FALSE);
      grdev_setcolor (COLOR_RGB, BlackTable_m);
      grdev_circle   (radius);

	/*  Fill only  */
	} else if (p->fill[ipart]==FILL_ON) {
      /*  Draw fill  */
      grdev_setfill  (TRUE);
      grdev_setcolor (p->ColorModel[ipart], p->ColorTable[ipart]);
      grdev_circle   (radius);

      /* Reset color to black */
      grdev_setcolor (COLOR_RGB, BlackTable_m);

	/*  Both outline and fill  */
	} else if (p->fill[ipart]==FILL_ALL) {
      /*  Draw fill  */
      grdev_setfill  (TRUE);
      grdev_setcolor (p->ColorModel[ipart], p->ColorTable[ipart]);
      grdev_circle   (radius);

      /* Draw outline  */
      grdev_setfill  (FALSE);
      grdev_setcolor (COLOR_RGB, BlackTable_m);
      grdev_circle   (radius);
	/*  This shouldn't happen  */
	} else {
      printf ("INTERNAL ERROR (plot_coord): Fill value is out of range.\n");
      CleanAfterError();
	}
}