xgpaint3d.c

xgrafix 是PTSG模拟程序中的图形截面库 改版本是最新版本

#include "xgrafixint.h"

#define deg2rad(x) (double)((x)*M_PI/180.0)
#define rad2deg(x) (double)((x)*180.0/M_PI)

#define  BACK         0
#define  FRONT        1

#define  ABOVE        0
#define  UPPER        1
#define  LOWER        2

#define  FIRST_QUAD   0
#define  SECOND_QUAD  1
#define  THIRD_QUAD   2
#define  FOURTH_QUAD  3

typedef enum position {NONE, TOP, BELOW} POSITION;

int CheckHiddenCrosshair(WindowType theWindow, int direction, 
			 int xIndex, int yIndex, int LeaveIndex,
			 double tan_phi, double theOtherIndexStep)
{
  int theCurrentIndex, theLastIndex, theIndexStep = 1;
  LabelType theLabel = theWindow->label;
  DataType theData = (DataType)malloc(sizeof(struct data));
  double theLineZStep, theLineZValue, theOtherIndex, theSurfZValue;
  double xmin, xmax, ymin, ymax, zmin, zmax, rho, phi;
  int xindex1, xindex2, yindex1, yindex2;
  POSITION linePos = NONE;
  double theta = theWindow->theta;

  phi = rad2deg(atan(tan_phi));
  xmin = theWindow->label->X_Min;
  xmax = theWindow->label->X_Max;
  ymin = theWindow->label->Y_Min;
  ymax = theWindow->label->Y_Max;
  zmin = theWindow->label->Z_Min;
  zmax = theWindow->label->Z_Max;
  theData->z = theWindow->data->z;
  theData->zArraySize = theWindow->data->zArraySize;
  theData->zArrayOffset = theWindow->data->zArrayOffset;

  switch(direction) {
  case 0:                    /* The X direction */
    theCurrentIndex = xIndex;
    theOtherIndex = yIndex;
    theLastIndex = LeaveIndex;
    if (xIndex > LeaveIndex) {
      theIndexStep = -1;
    }
    rho = 1.0/(sin(deg2rad(theta))*fabs(cos(deg2rad(phi))));
    theLineZStep = rho*cos(deg2rad(theta))*(zmax-zmin)/
      (double)(theWindow->mend-1-theWindow->mstart);
    break;
  case 1:                    /* the Y direction */
    theCurrentIndex = yIndex;
    theOtherIndex = xIndex;
    theLastIndex = LeaveIndex;
    if (yIndex > LeaveIndex) {
      theIndexStep = -1;
    }
    rho = 1.0/(sin(deg2rad(theta))*fabs(sin(deg2rad(phi))));
    theLineZStep = rho*cos(deg2rad(theta))*(zmax-zmin)/
      (double)(theWindow->nend-1-theWindow->nstart);
    break;
  }

  theLineZValue = SurfaceZ(theData,xIndex,yIndex)*theLabel->Z_Scale;

  theCurrentIndex += theIndexStep;
  theLastIndex += theIndexStep;
  for( ; theCurrentIndex != theLastIndex ; theCurrentIndex += theIndexStep) {
    theLineZValue += theLineZStep;

    theOtherIndex += theOtherIndexStep;
    if(direction == 0)  {
      xindex1 = theCurrentIndex;
      yindex1 = theOtherIndex + 1;
      xindex2 = theCurrentIndex;
      yindex2 = theOtherIndex;
    } else {
      xindex1 = theOtherIndex + 1;
      yindex1 = theCurrentIndex;
      xindex2 = theOtherIndex ;
      yindex2 = theCurrentIndex;
    }

    if ((int)theOtherIndex == theOtherIndex)
      theSurfZValue = SurfaceZ(theData,xindex2,yindex2)*theLabel->Z_Scale;
    else
      theSurfZValue = - (((int)theOtherIndex - theOtherIndex)*
			 SurfaceZ(theData,xindex1,yindex1))*theLabel->Z_Scale +
			   (((int)(fabs(theOtherIndex)+1)-theOtherIndex)*
			    SurfaceZ(theData,xindex2,yindex2))
			     *theLabel->Z_Scale;

    if(theSurfZValue > theLineZValue) {
      if (linePos == TOP)
	return 0;
      linePos = BELOW;
    } else if (theSurfZValue < theLineZValue) {
      if (linePos == BELOW)
	return 0;
      linePos = TOP;
    } else {
      return 0;
    }
  }

  return 1;
}

void DrawCrosshair(WindowType theWindow, int xindex, int yindex)
{
  Pixmap pm;
  XRectangle      theClientArea;
  Tk_Window tkwin = theWindow->tkwin;
  DataType theData = theWindow->data;
  LabelType theLabel = theWindow->label;
  MeshType **theMesh = theWindow->theMesh;
  double xtemp, ytemp;
  double xmin, xmax, ymin, ymax, phi, theta, tan_phi, gamma, tan_180minusphi
    , theOtherIndexStep1, theOtherIndexStep2, tan_360minusphi;
  int result = 0, LeaveIndex1, LeaveIndex2;

  theClientArea.x = 0;
  theClientArea.y = 0;
  theClientArea.width = Tk_Width(tkwin);
  theClientArea.height = Tk_Height(tkwin);
  XSetClipRectangles(theDisplay, theWindow->xwingc, 0, 0, &theClientArea, 
		     1, Unsorted);

  pm = XCreatePixmap(theDisplay, Tk_WindowId(tkwin), Tk_Width(tkwin),
		     Tk_Height(tkwin), Tk_Depth(tkwin));
  XCopyArea(theDisplay, theWindow->pixmap_buffer, pm, theWindow->xwingc, 0,
	    0, Tk_Width(tkwin), Tk_Height(tkwin), 0, 0);

  phi = theWindow->phi;
  theta = theWindow->theta;
  xmin = theLabel->X_Min;
  xmax = theLabel->X_Max;
  ymin = theLabel->Y_Min;
  ymax = theLabel->Y_Max;

  theWindow->data->xCoord = xindex;
  theWindow->data->yCoord = yindex;

  if ((xindex==theWindow->mstart && 90.0 < phi && phi < 270.0) ||
      (xindex==theWindow->mend-1 && (phi < 90.0 || phi > 270.0)) ||
      (yindex==theWindow->nend-1 && 0.0 < phi && phi < 180.0) ||
      (yindex==theWindow->nstart && 180.0 < phi && phi < 360.0) ||
      (theta == 0.0) ||
      (theta == 180.0)) {

    result = 1;
  } else {

    tan_phi = tan(deg2rad(phi))*(double)(theWindow->nend-1-theWindow->nstart)/(double)(theWindow->mend-1-theWindow->mstart);
    tan_180minusphi = tan(deg2rad(180.0-phi))*(double)(theWindow->nend-1-theWindow->nstart)/(double)(theWindow->mend-1-theWindow->mstart);
    tan_360minusphi = tan(deg2rad(360.0-phi))*(double)(theWindow->nend-1-theWindow->nstart)/(double)(theWindow->mend-1-theWindow->mstart);

    if (phi == 0.0) {
      result = CheckHiddenCrosshair(theWindow, 0, xindex, yindex,
				    theWindow->mend-1, tan_phi, 0);
    } else if (0.0 < phi && phi < 90.0) {
      gamma = rad2deg(atan((double)(theWindow->nend-1-yindex)/
			   (double)(theWindow->mend-1-xindex)*
			   (double)(theWindow->mend-1-theWindow->mstart)/
			   (double)(theWindow->nend-1-theWindow->nstart)));
      if (gamma >= phi) {
	LeaveIndex1 = theWindow->mend-1;
	theOtherIndexStep1 = tan_phi;

	LeaveIndex2 = tan_phi * (theWindow->mend-1 - xindex) + yindex;
	theOtherIndexStep2 = 1.0/tan_phi;
      } else {
	LeaveIndex1 = (theWindow->nend-1-yindex)/tan_phi + xindex;
	theOtherIndexStep1 = tan_phi;

	LeaveIndex2 = theWindow->nend-1;
	theOtherIndexStep2 = 1.0/tan_phi;
      }

      result = (CheckHiddenCrosshair(theWindow, 0, xindex, 
				     yindex, LeaveIndex1, tan_phi, 
				     theOtherIndexStep1) &&
		CheckHiddenCrosshair(theWindow, 1, xindex, 
				     yindex, LeaveIndex2, tan_phi,
				     theOtherIndexStep2));
    } else if (phi == 90.0) {
      result = CheckHiddenCrosshair(theWindow, 1, xindex, yindex,
				    theWindow->nend-1, tan_phi, 0);
    } else if (90.0 < phi && phi < 180.0) {
      gamma = rad2deg(atan((double)(theWindow->nend-1-yindex)/
			   (double)(xindex-theWindow->mstart)*
			   (double)(theWindow->mend-1-theWindow->mstart)/
			   (double)(theWindow->nend-1-theWindow->nstart)));
      if (gamma >= (180.0 - phi)) {
	LeaveIndex1 = theWindow->mstart;
	theOtherIndexStep1 = tan_180minusphi;

	LeaveIndex2 = (xindex-theWindow->mstart)*tan_180minusphi+yindex;
	theOtherIndexStep2 = -1.0/tan_180minusphi;
      } else {
	LeaveIndex1 = xindex-(theWindow->nend-1-yindex)/tan_180minusphi + 1;
	theOtherIndexStep1 = tan_180minusphi;

	LeaveIndex2 = theWindow->nend-1;
	theOtherIndexStep2 = -1.0/tan_180minusphi;
      }
      result = (CheckHiddenCrosshair(theWindow, 0, xindex,
				     yindex, LeaveIndex1, tan_phi,
				     theOtherIndexStep1) &&
		CheckHiddenCrosshair(theWindow, 1, xindex,
				     yindex, LeaveIndex2, tan_phi,
				     theOtherIndexStep2));
    } else if (phi == 180.0) {
      result = CheckHiddenCrosshair(theWindow, 0, xindex, yindex,
				    theWindow->mstart, tan_phi, 0);
    } else if (180.0 < phi && phi < 270.0) {
      gamma = rad2deg(atan((double)(yindex-theWindow->nstart)/
			   (double)(xindex-theWindow->mstart)*
			   (double)(theWindow->mend-1-theWindow->mstart)/
			   (double)(theWindow->nend-1-theWindow->nstart)));
      if (gamma >= (phi - 180.0)) {
	LeaveIndex1 = theWindow->mstart;
	theOtherIndexStep1 = 0.0 - tan_phi;

	LeaveIndex2 = yindex - (xindex - theWindow->mstart)*tan_phi + 1;
	theOtherIndexStep2 = -1.0/tan_phi;
      } else {
	LeaveIndex1 = xindex - (yindex - theWindow->nstart)/tan_phi + 1;
	theOtherIndexStep1 = 0.0 - tan_phi;

	LeaveIndex2 = theWindow->nstart;
	theOtherIndexStep2 = -1.0/tan_phi;
      }
      result = (CheckHiddenCrosshair(theWindow, 0, xindex,
				     yindex, LeaveIndex1, tan_phi,
				     theOtherIndexStep1) &&
		CheckHiddenCrosshair(theWindow, 1, xindex,
				     yindex, LeaveIndex2, tan_phi,
				     theOtherIndexStep2));
    } else if (phi == 270.0) {
      result = CheckHiddenCrosshair(theWindow, 1, xindex, yindex,
				    theWindow->nstart, tan_phi, 0);
    } else {
      gamma = rad2deg(atan((double)(yindex-theWindow->nstart)/
			   (double)(theWindow->mend-1-xindex)*
			   (double)(theWindow->mend-1-theWindow->mstart)/
			   (double)(theWindow->nend-1-theWindow->nstart)));
      if (gamma >= 360.0 - phi) {
	LeaveIndex1 = theWindow->mend-1;
	theOtherIndexStep1 = -tan_360minusphi;

	LeaveIndex2 = yindex - (theWindow->mend-1-xindex) * tan_360minusphi+1;
	theOtherIndexStep2 = 1.0/tan_360minusphi;
      } else {
	LeaveIndex1 = (yindex - theWindow->nstart)/tan_360minusphi + xindex;
	theOtherIndexStep1 = -tan_360minusphi;
	
	LeaveIndex2 = theWindow->nstart;
	theOtherIndexStep2 = 1.0/tan_360minusphi;
      }
      result = (CheckHiddenCrosshair(theWindow, 0, xindex,
				     yindex, LeaveIndex1, tan_phi,
				     theOtherIndexStep1) &&
		CheckHiddenCrosshair(theWindow, 1, xindex,
				     yindex, LeaveIndex2, tan_phi,
				     theOtherIndexStep2));
    }
  }
			 
  xtemp = theWindow->c1 * theMesh[xindex][yindex].x + theWindow->d1 - 4;
  ytemp = theWindow->c2 * theMesh[xindex][yindex].y + theWindow->d2 - 4;
  XSetForeground(theDisplay, theWindow->xwingc, theWhitePixel);

  if(result) {
    XFillArc(theDisplay, pm, theWindow->xwingc, (int)xtemp, 
	     (int)ytemp, 8, 8, 0, 23040);

  } else {
    XDrawArc(theDisplay, pm, theWindow->xwingc, (int)xtemp, 
	     (int)ytemp, 8, 8, 0, 23040);
  }

  XCopyArea(theDisplay, pm, Tk_WindowId(tkwin), 
	    theWindow->xwingc, 0, 0, Tk_Width(tkwin), Tk_Height(tkwin), 
	    0, 0);
  XFlush(theDisplay);
  XFreePixmap(theDisplay, pm);
  XSetClipMask(theDisplay, theWindow->xwingc, None);

  sprintf(TclCommand,"Update3dCrosshairValues \"%s\" %g %g %g\n",
	  theWindow->plot_title,SurfaceX(theData,xindex)*theLabel->X_Scale,
	  SurfaceY(theData,yindex)*theLabel->Y_Scale,
	  SurfaceZ(theData,xindex,yindex)*theLabel->Z_Scale);
  Tcl_Eval(interp,TclCommand);

}

void DrawWall(WindowType theWindow, MeshType *theBox, double c1, double d1, 
	      double c2, double d2, int wall, int side)
{
  XPoint          thePoint[4];
  
  if (wall == 0) {
    thePoint[0].x = ToShort(c1 * theBox[0].x + d1);
    thePoint[0].y = ToShort(c2 * theBox[0].y + d2);
    thePoint[1].x = ToShort(c1 * theBox[1].x + d1);
    thePoint[1].y = ToShort(c2 * theBox[1].y + d2);
    thePoint[2].x = ToShort(c1 * theBox[2].x + d1);
    thePoint[2].y = ToShort(c2 * theBox[2].y + d2);
    thePoint[3].x = ToShort(c1 * theBox[3].x + d1);
    thePoint[3].y = ToShort(c2 * theBox[3].y + d2);
  } else if (wall == 1) {
    thePoint[0].x = ToShort(c1 * theBox[0].x + d1);
    thePoint[0].y = ToShort(c2 * theBox[0].y + d2);
    thePoint[1].x = ToShort(c1 * theBox[1].x + d1);
    thePoint[1].y = ToShort(c2 * theBox[1].y + d2);
    thePoint[2].x = ToShort(c1 * theBox[6].x + d1);
    thePoint[2].y = ToShort(c2 * theBox[6].y + d2);
    thePoint[3].x = ToShort(c1 * theBox[5].x + d1);
    thePoint[3].y = ToShort(c2 * theBox[5].y + d2);
  } else if (wall == 2) {
    thePoint[0].x = ToShort(c1 * theBox[1].x + d1);
    thePoint[0].y = ToShort(c2 * theBox[1].y + d2);
    thePoint[1].x = ToShort(c1 * theBox[2].x + d1);
    thePoint[1].y = ToShort(c2 * theBox[2].y + d2);
    thePoint[2].x = ToShort(c1 * theBox[7].x + d1);
    thePoint[2].y = ToShort(c2 * theBox[7].y + d2);
    thePoint[3].x = ToShort(c1 * theBox[6].x + d1);
    thePoint[3].y = ToShort(c2 * theBox[6].y + d2);
  } else if (wall == 3) {
    thePoint[0].x = ToShort(c1 * theBox[2].x + d1);
    thePoint[0].y = ToShort(c2 * theBox[2].y + d2);
    thePoint[1].x = ToShort(c1 * theBox[3].x + d1);
    thePoint[1].y = ToShort(c2 * theBox[3].y + d2);
    thePoint[2].x = ToShort(c1 * theBox[4].x + d1);
    thePoint[2].y = ToShort(c2 * theBox[4].y + d2);
    thePoint[3].x = ToShort(c1 * theBox[7].x + d1);
    thePoint[3].y = ToShort(c2 * theBox[7].y + d2);
  } else if (wall == 4) {
    thePoint[0].x = ToShort(c1 * theBox[3].x + d1);
    thePoint[0].y = ToShort(c2 * theBox[3].y + d2);
    thePoint[1].x = ToShort(c1 * theBox[0].x + d1);
    thePoint[1].y = ToShort(c2 * theBox[0].y + d2);
    thePoint[2].x = ToShort(c1 * theBox[5].x + d1);
    thePoint[2].y = ToShort(c2 * theBox[5].y + d2);
    thePoint[3].x = ToShort(c1 * theBox[4].x + d1);
    thePoint[3].y = ToShort(c2 * theBox[4].y + d2);
  } else if (wall == 5) {
    thePoint[0].x = ToShort(c1 * theBox[4].x + d1);
    thePoint[0].y = ToShort(c2 * theBox[4].y + d2);
    thePoint[1].x = ToShort(c1 * theBox[5].x + d1);
    thePoint[1].y = ToShort(c2 * theBox[5].y + d2);
    thePoint[2].x = ToShort(c1 * theBox[6].x + d1);
    thePoint[2].y = ToShort(c2 * theBox[6].y + d2);
    thePoint[3].x = ToShort(c1 * theBox[7].x + d1);
    thePoint[3].y = ToShort(c2 * theBox[7].y + d2);
  }
  if (side == BACK) {
    XSetForeground(theDisplay, theWindow->xwingc, theBlackPixel);
    XFillPolygon(theDisplay, theWindow->pixmap_buffer,
		 theWindow->xwingc, thePoint, 4, Complex, CoordModeOrigin);
  }
  XSetForeground(theDisplay, theWindow->xwingc, theWhitePixel);
  XDrawLines(theDisplay, theWindow->pixmap_buffer, theWindow->xwingc, thePoint,
	     4, CoordModeOrigin);
}

/****************************************************************/

void Paint_ThreeD_Window(WindowType theWindow)
{
  DataType        data_Ptr;
  LabelType       label_Ptr;
  MeshType      **theMeshPoint=theWindow->theMesh, theBox[8];
  StringType      theString[8];
  XPoint          thePolygon[5];
  Tk_Window       tkwin;
  XRectangle      theClientArea;

  double           n1, n2, n3, v1, v2, v3;
  double           u1, u2, u3, up1, up2, up3;
  double           a1, a2, a3, b1, b2, b3;
  double           xprojmax, xprojmin;
  double           yprojmax, yprojmin;
  double           mag, c1, c2, d1, d2;
  double           xtemp, ytemp, ztemp;

  int             m, n, mstart, mend, nend, nstart;
  int             i, j, k, x1, y1, x2, y2, colornumber;
  int             phi_flag, theta_flag, i1, j1, k1, maxlen;
  int             xstart, xend, xstep, ystart, yend, ystep;

  tkwin = theWindow->tkwin;

  data_Ptr = theWindow->data;
  label_Ptr = theWindow->label;

  if (!Is_Scatter(data_Ptr)) {
    m = *(data_Ptr->mpoints);
  } 
  else {
    m = *(data_Ptr->npoints);
  }
  n = *(data_Ptr->npoints);
  
  if(m<3 || n<2) return;

  if (theWindow->theMesh == NULL) {
    theWindow->theMesh = (MeshType **) malloc(m * sizeof(MeshType *));
    theMeshPoint = theWindow->theMesh;
    for (i = 0; i < m; i++) {
      theMeshPoint[i] = (MeshType *) malloc(n * sizeof(MeshType));
    }

    theWindow->xSize = m;
    theWindow->ySize = n;
  } else if ((m != theWindow->xSize) || n != (theWindow->ySize)) {
    for (i=0; i<theWindow->xSize; i++) {