Editing Point-in-polygon

; Language: C++
; Purpose: Determines whether or not a point is in a polygon. Jumps in complex arguments are preserved.
; Author: James R. Craig
; GW Gourmet Functions Used: [[logZp1oZm1]]

 #include <math.h>
 #include <complex>
 #include <iostream>
 typedef std::complex<double> cmplex;
 const double CLOSE_TO_LINE=1.02;// ellipse radius around line which determines closeness
 const double REALSMALL=1e-10; 	//almost zero
 /**************************************************************************
 IsInsidePolygon:
 Returns if a point z is in polygon defined by N vertices stored in 
 array ze[]. The last polygon point is NOT repeated.
 Adapted from Paul Bourke <http://astronomy.swin.edu.au/pbourke/>
 Uses function logZm1oZp1(Z), available at www.analyticelements.org 
 "Groundwater Gourmet"
 Author: James R. Craig
 -------------------------------------------------------------------------*/
 bool IsInsidePolygon(const cmplex &z, const cmplex *zv, const int N)
 {
 	double y(z.imag());
 	double x(z.real()); 
 	int  i;
 	bool close(false); 
 	if (abs(zv[0]-zv[N])>REALSMALL){
 		cout<<"PolyIsInside : Bad polygon"<<endl;
 		exit(0);
 	} 
 	//check if point is close to boundary of polygon
 	for (i=0; i<N; i++){
 		if (((abs(z-zv[i])+abs(z-zv[i+1]))/(abs(zv[i+1]-zv[i])))<CLOSE_TO_LINE){
 			close=true;break;
 		}
 	}
 	//if not close, use faster ray-tracing algorithm
 	//adapted from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
 	if (!close)
 	{
 		bool	 in(false);
 		int 	 j;
 		for (i=0, j=N-1; i<N; j=i++) {
 			if ((((zv[i].imag()<=y) && (y<zv[j].imag())) ||
 					 ((zv[j].imag()<=y) && (y<zv[i].imag()))) &&
 				 (x<(zv[j].real()-zv[i].real())*(y-zv[i].imag())/
 				 (zv[j].imag()-zv[i].imag())+zv[i].real())){
 				in=!in;
 			}
 		}
 		return in;
 	}
 	//if close, use complex winding algorithm
        //this is consistent with the jump calculation for line elements
 	//sum ln (Z-1)/(Z+1) should be zero if outside
 	else
 	{
 		cmplex Z;
 		double sum(0);
 		for (int i=0;i<N; i++)
 		{
 			Z=(z-0.5*(zv[i]+zv[i+1]))/(0.5*(zv[i+1]-zv[i]));
 			sum+=logZm1oZp1(Z).imag();
 		}
 		return (fabs(sum)<REALSMALL ? false : true);
 	}
 }