# include <iostream.h>
# include <math.h>
#include"headers.h"
double dot_product1(xyz *a,xyz *b)
{
	double d;
	d=(a->x*b->x)+(a->y*b->y)+(a->z*b->z);
	return(d);
}
void cross_product1(xyz *o,xyz *a,xyz *b)
{
	double d;
	o->x = (a->y * b->z) - (a->z * b->y);
	o->y = (a->z * b->x) - (a->x * b->z);
	o->z = (a->x * b->y) - (a->y * b->x);
	d= sqrt( o->x*o->x+o->y*o->y+o->z*o->z);
}

double main()
{
	xy p1,p2,p3,p4;
	xy *p[3];
	//make plane
double d,max_coefficient,magu,magv,samples,cellu,cellv;
xyz pnta,pntb,pntc,pntd,pnts,u,v,n;
samples=3;
pnta.x=0;
pnta.y=50;
pnta.z=70;
pntb.x=0;
pntb.y=50;
pntb.z=-70;
pntc.x=0;
pntc.y=-10;
pntc.z=-70;
pntd.x=0;
pntd.y=-10;
pntd.z=70;


u.x=pntb.x-pnta.x;
u.y=pntb.y-pnta.y;
u.z=pntb.z-pnta.z;
v.x=pntc.x-pntb.x;
v.y=pntc.y-pntb.y;
v.z=pntc.z-pntb.z;

//cout<<"u x "<<u.x<<"\n";
//cout<<"u y "<<u.y<<"\n";
//cout<<"u z "<<u.z<<"\n\n";
//cout<<"v x "<<v.x<<"\n";
//cout<<"v y "<<v.y<<"\n";
//cout<<"v z "<<v.z<<"\n\n";

magu=u.normalize();
magv=v.normalize();

cellu=magu/samples;
cellv=magv/samples;

pnts.x=pnta.x+(cellu/2)*u.x+(cellv/2)*v.x;
pnts.y=pnta.y+(cellu/2)*u.y+(cellv/2)*v.y;
pnts.z=pnta.z+(cellu/2)*u.z+(cellv/2)*v.z;

cout<<"cell 11  x "<<v.x<<"\n";
cout<<"cell 11  y "<<v.y<<"\n";
cout<<"cell 11  z "<<v.z<<"\n\n";










//-----------------------------------------------------------------------------------------------
/* finds the normal of the plane containing the quardilatral */

cross_product1(&n,&u,&v);
cout<<"n x "<<n.x<<"\n";
cout<<"n y "<<n.y<<"\n";
cout<<"n z "<<n.z<<"\n";
/* calculates the coefficient d of the plane equation */

d=-(n.x*pnta.x+n.y*pnta.y+n.z*pnta.z);
cout<<"d "<<d<<"\n";

/* projects the the 3d points on to a 2d plane,by selecting 
the normal component that has the highest coefficient(absolute value)*/

if((fabs(n.x)>fabs(n.y))&&(fabs(n.x)>fabs(n.z)))
{
	max_coefficient=1;
	p1.x=pnta.y;
	p1.y=pnta.z;
	p2.x=pntb.y;
	p2.y=pntb.z;
	p3.x=pntc.y;
	p3.y=pntc.z;
	p4.x=pntd.y;
	p4.y=pntd.z;
}
else if((fabs(n.y)>fabs(n.x))&&(fabs(n.y)>fabs(n.z)))
{
	max_coefficient=2;
	p1.x=pnta.x;
	p1.y=pnta.z;
	p2.x=pntb.x;
	p2.y=pntb.z;
	p3.x=pntc.x;
	p3.y=pntc.z;
	p4.x=pntd.x;
	p4.y=pntd.z;
}
else
{
	max_coefficient=3;
	p1.x=pnta.x;
	p1.y=pnta.y;
	p2.x=pntb.x;
	p2.y=pntb.y;
	p3.x=pntc.x;
	p3.y=pntc.y;
	p4.x=pntd.x;
	p4.y=pntd.y;
}

cout<<"p1 x"<<p1.x<<"\n";
cout<<"p1 y"<<p1.y<<"\n";

// check for ray-plane intersection

int i,sign=0,next_sign=0,crossing=0;
xy hit2d;
xyz rayo,rayd,hit;
double vd,vo,t=-1;

rayo.x=2;
rayo.y=3;
rayo.z=4;
rayd.x=.577;
rayd.y=.577;
rayd.z=.577;
vd=dot_product1(&n,&rayd);
if(vd==0)
cout<<"distance = 0.0";
vo=-(dot_product1(&n,&rayo)+d);
t=vo/vd;
//cout<<"T "<<t<<"\n";
if(t<0)
return(0.0);
cout<<"t1 0\n";

/* find whether the point of intersection is inside or
outside the polygon */


hit.x=rayo.x+t*rayd.x;
hit.y=rayo.y+t*rayd.y;
hit.z=rayo.z+t*rayd.z;
cout<<"n x "<<hit.x<<"\n";
cout<<"n y "<<hit.y<<"\n";
cout<<"n z "<<hit.z<<"\n";
if(max_coefficient==1)
{
	hit2d.x=hit.y;
	hit2d.y=hit.z;

}
else if(max_coefficient==2)
{
	hit2d.x=hit.x;
	hit2d.y=hit.z;
}
else
{
	hit2d.x=hit.x;
	hit2d.y=hit.y;
}
//cout<<"hit2d x"<<hit2d.x<<"\n";
//cout<<"hit2d y"<<hit2d.y<<"\n";
/* after projecting the plane on to 2d...check whether 
the point lies inside or outside the polygon */

/*make the point of intersection as the origin and translate 
the vertices of the polygon relative to the point of intersection*/
cout<<"p1 x "<<p1.x<<"\n";
cout<<"p1 y "<<p1.y<<"\n";
cout<<"p2 x "<<p2.x<<"\n";
cout<<"p2 y "<<p2.y<<"\n";
p1.x=p1.x-hit2d.x;
p1.y=p1.y-hit2d.y;
p2.x=p2.x-hit2d.x;
p2.y=p2.y-hit2d.y;
p3.x=p3.x-hit2d.x;
p3.y=p3.y-hit2d.y;
p4.x=p4.x-hit2d.x;
p4.y=p4.y-hit2d.y;

// assign pointers to make looping easier

p[0]=&p1;
p[1]=&p2;
p[2]=&p3;
p[3]=&p4;
cout<<"2p1 x "<<p1.x<<"\n";
cout<<"2p1 y "<<p1.y<<"\n";
cout<<"2p2 x "<<p2.x<<"\n";
cout<<"2p2 y "<<p2.y<<"\n";
for(i=0;i<4;i++)
{
	if(p[i]->y>=0)
		sign=1;
	else
		sign=-1;
cout<<"sign "<<sign<<"\n";
	if(p[(i+1)%4]->y>=0)
		next_sign=1;
	else
		next_sign=-1;
cout<<"next sign "<<next_sign<<"\n";
	if(sign != next_sign)
	{
		if(((p[i]->x)>=0)&&(p[(i+1)%4]->x>=0))
		{
			crossing=crossing+1;
		}
		else if(((p[i]->x)<0)&&(p[(i+1)%4]->x<0))
		{}
		else
		{
			if(((p[i]->x - p[i]->y) *( p[(i+1)%4]->x - p[i]->x) / (p[(i+1)%4]->y - p[i]->y))>0)
				crossing=crossing+1;
		}
	}
}

if(crossing%2==0)
	return(0.0);

 else

   return(t);
}