plotc.c

一个很好的分子动力学程序

/*
************************************************************************
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.
*/

/*
************************************************************************
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



/*
************************************************************************
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;
   float bondlo, bondhi, dispscale;
   int orient[3];
   float *radius;
   int DEVflag;
	float PageSize[2];
   } 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_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 sflag);

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);

   /*  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_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 = FALSE;
	BOOLEAN ReadColorFlag;
	

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

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

	if (!strcmpi(TokenStr,"FILL"))
		{
		FillFlag  = TRUE;
		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;
	return t ;
	}




void draw_part
(
Particle_t *a,
int i,
PlotInfo_t *pl,
double xstart,
double ystart,
double scale,
int hor,