/*!\file IsInPoly.c
 * \brief:  from a contour, determine which nodes are in  domain.
 */

#include <math.h>
#include "../../toolkits/toolkits.h"
#include "./exp.h"


int IsInPoly(Vec in,double* xc,double* yc,int numgrids,double* x,double* y,int nods, int edgevalue){

	int i,j;
	double x0,y0;
	double value;

	/*Go through all grids of the mesh:*/
	for (i=MPI_Lowerrow(nods);i<MPI_Upperrow(nods);i++){
		/*pick up grid: */
		x0=x[i];
		y0=y[i];
		if (pnpoly(numgrids,xc,yc,x0,y0,edgevalue)){
			value=1;
		}
		else{
			value=0;
		}
		VecSetValues(in,1,&i,&value,INSERT_VALUES);
	}
}

int IsOutsidePoly(Vec out,double* xc,double* yc,int numgrids,double* x,double* y,int nods,int edgevalue){

	int i,j;
	double x0,y0;
	double value;

	/*Go through all grids of the mesh:*/
	for (i=MPI_Lowerrow(nods);i<MPI_Upperrow(nods);i++){
		/*pick up grid: */
		x0=x[i];
		y0=y[i];
		if (pnpoly(numgrids,xc,yc,x0,y0,edgevalue)){
			value=0;
		}
		else{
			value=1;
		}
		VecSetValues(out,1,&i,&value,INSERT_VALUES);
	}
}

int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue) {
	int i, j, c = 0;
	double n1, n2, normp, scalar;
	
	/*first test, are they colinear? if yes, is the point in the middle of the segment*/
	if (edgevalue != 2 ){
		for (i = 0, j = npol-1; i < npol; j = i++) {
			n1=pow(yp[i]-yp[j],2.0)+pow(xp[i]-xp[j],2.0);
			n2=pow(y-yp[j],2.0)+pow(x-xp[j],2.0);
			normp=pow(n1*n2,0.5);
			scalar=(yp[i]-yp[j])*(y-yp[j])+(xp[i]-xp[j])*(x-xp[j]);

			if (scalar == normp){
				if (n2<=n1){
					c = edgevalue;
					return c;
				}
			}
		}
	}
	/*second test : point is neither on a vertex, nor on a side, where is it ?*/
	for (i = 0, j = npol-1; i < npol; j = i++) {
		if ((((yp[i]<=y) && (y<yp[j])) ||
					((yp[j]<=y) && (y<yp[i]))) &&
				(x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])){
			c = !c;
		}
	}
	return c;
}
