/************************************************************************************************
* User-defined function for arbitrary polygonal holes
*
* Called from inside a groundplane definition as:
*   +  hole user1 (x1,y1,x2,y2,x3,y3,...)
*
* Defines a polygon by a set of three or more tuples giving (x,y) coordinates for the polygon
* vertices. Removes all nodes from the groundplane that are inside or on the edge of the polygon.
*
* The function uses a standard ray-casting algorithm (http://en.wikipedia.org/wiki/Point_in_polygon)
* to determine if each groundplane node is inside or outside the polygon: for each node, extend a ray
* in any direction (here in the -x direction) and count the number of polygon edges it crosses. If
* the count is odd, the point is inside the polygon, if the count is even, the point is outside the
* polygon. This algorithm handles self-crossing polygons, like hourglasses, or pentagrams, but is not
* recommended for them as it creates groundplane islands, disconnected or linked by a single segment.
*
* Note: This user function assumes ground plane and hole are both at constant z. Feel free to
*       generalize it if you want to do the work.
*
*************************************************************************************************/
hole_user1(holep, gp, relx, rely, relz, units)
GROUNDPLANE *gp;
HoleList *holep;
double relx,rely,relz,units;
{
  int i, j, e;

  NODES *origin;
  int gn1, gn2;
  double gx0, gy0;
  double gd1, gd2;
  double ux1, uy1, ux2, uy2;

  int nedges;
  double ex0, ey0, ex1, ey1;

  int inPolygon;


  if ( ((holep->numvals)%2 != 0) || (holep->numvals < 6) )    // Check for correct number of inputs (expect even number greater than 5)
    hole_error("Need an even number of inputs for at least three (x,y) pairs.","(x,y)",holep);

  nedges = holep->numvals / 2;    // Number of edges of hole polygon

  // Extract useful variables from the groundplane structure
  origin = gp->pnodes[0][0];
  gx0 = origin->x / units;        // Origin corner of groundplane (x-coord)
  gy0 = origin->y / units;        // Origin corner of groundplane (y-coord)
  gn1 = gp->num_nodes1;            // Number of nodes along first side of groundplane
  gn2 = gp->num_nodes2;            // Number of nodes along second side of groundplane
  gd1 = gp->d1 / units;            // Groundplane segment length along first side
  gd2 = gp->d2 / units;            // Groundplane segment length along second side
  ux1 = gp->ux1;                // Unit vector between segments along first side (x-coord)
  uy1 = gp->uy1;                // Unit vector between segments along first side (y-coord)
  ux2 = gp->ux2;                // Unit vector between segments along second side (x-coord)
  uy2 = gp->uy2;                // Unit vector between segments along second side (y-coord)

  for(j = 0; j < gn2; j++) {    // For each ground plane node
    for(i = 0; i < gn1; i++) {
      inPolygon = FALSE;
      for(e = 0; e < nedges; e++) {    // For each edge of the polygon
        // Get the endpoints of the edge in the groundplane coordinate system (ex0,ey0) and (ex1,ey1)
        // (Possible confusion: I use a different convention from that of the groundplane struct which has axis 1 and axis 2
        //  and use x and y in the groundplane coordinate system for later thinking convenience.)
        ex0 = (holep->vals[0+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux1 + (holep->vals[1+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy1;
        ey0 = (holep->vals[0+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux2 + (holep->vals[1+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy2;
        if(e+1 == nedges) {
          ex1 = (holep->vals[0] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux1 + (holep->vals[1] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy1;
          ey1 = (holep->vals[0] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux2 + (holep->vals[1] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy2;
        } else {
          ex1 = (holep->vals[2+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux1 + (holep->vals[3+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy1;
          ey1 = (holep->vals[2+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux2 + (holep->vals[3+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy2;
        }
        // Toggle inPolygon if the leftward ray from the groundplane point intersects the polygon edge
        if(ey0*ey1 < 0) {
          if((ex0 - ey0*(ex1-ex0)/(ey1-ey0)) <= 0) {
            // If the polygon edge crosses the horizontal line and does it to the left of the node, toggle inside/outside status
            inPolygon = !inPolygon;
          }
        } else if( ((ey0 == 0) && (ey1 < 0) && (ex0 <= 0)) || ((ey1 == 0) && (ey0 < 0) && (ex1 <= 0))) {
          // If one endpoint of a polygon edge touches the horizontal to the left of the node and the other node is under the line, toggle inside/outside status
          inPolygon = !inPolygon;
        }
      }
      if(inPolygon) {
        // Delete node if inside polygon
        delete_node(gp->pnodes[i][j]);
      }
    }
  }

}