// requires triangle.c


// Matrix struct

typedef struct
{ 
  double m[2][3];
}
matrix;

typedef struct
{
  matrix m[3];
} 
reflection;

//-------------------------------------------------------
// Matrix multiplication functions

// returns with m1 storing m1*m2 (euclidean transforms)
void timesEquals(m1,m2)
     matrix *m1, *m2;
{
  double row[3];
  row[0]=m1->m[0][0]*m2->m[0][0]+m1->m[0][1]*m2->m[1][0];
  row[1]=m1->m[0][0]*m2->m[0][1]+m1->m[0][1]*m2->m[1][1];
  row[2]=m1->m[0][0]*m2->m[0][2]+m1->m[0][1]*m2->m[1][2]+m1->m[0][2];
  m1->m[0][0]=row[0];
  m1->m[0][1]=row[1];
  m1->m[0][2]=row[2];
  row[0]=m1->m[1][0]*m2->m[0][0]+m1->m[1][1]*m2->m[1][0];
  row[1]=m1->m[1][0]*m2->m[0][1]+m1->m[1][1]*m2->m[1][1];
  row[2]=m1->m[1][0]*m2->m[0][2]+m1->m[1][1]*m2->m[1][2]+m1->m[1][2];
  m1->m[1][0]=row[0];
  m1->m[1][1]=row[1];
  m1->m[1][2]=row[2];
}

pcomplex eucTimes(m1,z)
     matrix *m1; 
     pcomplex z;
{
  pcomplex ret;
  ret.x=m1->m[0][0]*z.x+m1->m[0][1]*z.y+m1->m[0][2];
  ret.y=m1->m[1][0]*z.x+m1->m[1][1]*z.y+m1->m[1][2];
  return ret;
}


// returns an array of reflections from a triangle
/* void reflect_from_ptriangle(T,r) */
/*      ptriangle T;  */
/*      reflection *r; */
/* { */
/*   int i,j; */
/*   matrix temp1,temp2; */
/*   double x,y,det; */
/*   { */
/*     temp1.m[0][0]=1; temp1.m[0][1]=0; temp1.m[0][2]=-T.z2.x; */
/*     temp1.m[1][0]=0; temp1.m[1][1]=1; temp1.m[1][2]=-T.z2.y; */
/*     x=T.z3.x-T.z2.x; */
/*     y=T.z3.y-T.z2.y; */
/*     temp2.m[0][0]=(x*x-y*y); temp2.m[0][1]=2*x*y; temp2.m[0][2]=0; */
/*     temp2.m[1][0]=2*x*y; temp2.m[1][1]=-(x*x-y*y); temp2.m[1][2]=0; */
/*     det=sqrt(-(temp2.m[0][0]*temp2.m[1][1]- */
/* 	       temp2.m[0][1]*temp2.m[1][0])); */
/*     for (i=0; i<2; i++) */
/*       for (j=0; j<3; j++) */
/* 	temp2.m[i][j]*=1/det; */
/*     r->m[0].m[0][0]=1; r->m[0].m[0][1]=0; r->m[0].m[0][2]=T.z2.x;  */
/*     r->m[0].m[1][0]=0; r->m[0].m[1][1]=1; r->m[0].m[1][2]=T.z2.y;  */
/*     timesEquals(&(r->m[0]),&temp2); */
/*     timesEquals(&(r->m[0]),&temp1); */
/*   }     */
/*   { */
/*     temp1.m[0][0]=1; temp1.m[0][1]=0; temp1.m[0][2]=-T.z3.x; */
/*     temp1.m[1][0]=0; temp1.m[1][1]=1; temp1.m[1][2]=-T.z3.y; */
/*     x=T.z1.x-T.z3.x; */
/*     y=T.z1.y-T.z3.y; */
/*     temp2.m[0][0]=(x*x-y*y); temp2.m[0][1]=2*x*y; temp2.m[0][2]=0; */
/*     temp2.m[1][0]=2*x*y; temp2.m[1][1]=-(x*x-y*y); temp2.m[1][2]=0; */
/*     det=sqrt(-(temp2.m[0][0]*temp2.m[1][1]- */
/* 	       temp2.m[0][1]*temp2.m[1][0])); */
/*     for (i=0; i<2; i++) */
/*       for (j=0; j<3; j++) */
/* 	temp2.m[i][j]*=1/det; */
/*     r->m[1].m[0][0]=1; r->m[1].m[0][1]=0; r->m[1].m[0][2]=T.z3.x;  */
/*     r->m[1].m[1][0]=0; r->m[1].m[1][1]=1; r->m[1].m[1][2]=T.z3.y;  */
/*     timesEquals(&(r->m[1]),&temp2); */
/*     timesEquals(&(r->m[1]),&temp1); */
/*   }     */
/*   { */
/*     temp1.m[0][0]=1; temp1.m[0][1]=0; temp1.m[0][2]=-T.z1.x; */
/*     temp1.m[1][0]=0; temp1.m[1][1]=1; temp1.m[1][2]=-T.z1.y; */
/*     x=T.z2.x-T.z1.x; */
/*     y=T.z2.y-T.z1.y; */
/*     temp2.m[0][0]=(x*x-y*y); temp2.m[0][1]=2*x*y; temp2.m[0][2]=0; */
/*     temp2.m[1][0]=2*x*y; temp2.m[1][1]=-(x*x-y*y); temp2.m[1][2]=0; */
/*     det=sqrt(-(temp2.m[0][0]*temp2.m[1][1]- */
/* 	      temp2.m[0][1]*temp2.m[1][0])); */
/*     for (i=0; i<2; i++) */
/*       for (j=0; j<3; j++) */
/* 	temp2.m[i][j]*=1/det; */
/*     r->m[2].m[0][0]=1; r->m[2].m[0][1]=0; r->m[2].m[0][2]=T.z1.x;  */
/*     r->m[2].m[1][0]=0; r->m[2].m[1][1]=1; r->m[2].m[1][2]=T.z1.y;  */
/*     timesEquals(&(r->m[2]),&temp2); */
/*     timesEquals(&(r->m[2]),&temp1); */
/*   }     */
/* } */

// returns an array of reflections from a triangle
void reflect_from_triangle(T,r)
     triangle T; 
     reflection *r;
{
  int i,j,k;
  matrix temp1,temp2;
  double x,y,det;
  for (k=0;k<3;k++){
    temp1.m[0][0]=1; temp1.m[0][1]=0; temp1.m[0][2]=-T.z[(k+1)%3].x;
    temp1.m[1][0]=0; temp1.m[1][1]=1; temp1.m[1][2]=-T.z[(k+1)%3].y;
    x=T.z[(k+2)%3].x-T.z[(k+1)%3].x;
    y=T.z[(k+2)%3].y-T.z[(k+1)%3].y;
    temp2.m[0][0]=(x*x-y*y); temp2.m[0][1]=2*x*y; temp2.m[0][2]=0;
    temp2.m[1][0]=2*x*y; temp2.m[1][1]=-(x*x-y*y); temp2.m[1][2]=0;
    det=sqrt(-(temp2.m[0][0]*temp2.m[1][1]-
	       temp2.m[0][1]*temp2.m[1][0]));
    for (i=0; i<2; i++)
      for (j=0; j<3; j++)
	temp2.m[i][j]*=1/det;
    r->m[k].m[0][0]=1; r->m[k].m[0][1]=0; r->m[k].m[0][2]=T.z[(k+1)%3].x; 
    r->m[k].m[1][0]=0; r->m[k].m[1][1]=1; r->m[k].m[1][2]=T.z[(k+1)%3].y; 
    timesEquals(&(r->m[k]),&temp2);
    timesEquals(&(r->m[k]),&temp1);
  }    
}

void matrix_identity(m)
  matrix *m;
{
  m->m[0][0]=m->m[1][1]=1;
  m->m[0][1]=m->m[0][2]=m->m[1][0]=m->m[1][2]=0;
}

/* int det5(a,b,c) */
/*      pcomplex *a, *b, *c; */
/* { */
/*   doubleList p[6]; */
/*   doubleList dl; */
/*   int ret; */

/*   p[0]=storeProduct(a->x, b->y); */
/*   p[1]=storeProduct(b->x, c->y); */
/*   p[2]=storeProduct(c->x, a->y); */
/*   p[3]=storeProduct(-a->x, c->y); */
/*   p[4]=storeProduct(-b->x, a->y); */
/*   p[5]=storeProduct(-c->x, b->y); */
/*   dl=DLsumDestroy(p,6); */
/*   DLmag(&dl); */
/*   DLsort(&dl); */
/*   ret=DLsign(&dl); */
/*   DLdestroy(&dl); */
/*   return ret; */
/* }     */

// returns the sign of the determinant, that is the orientation of the triangle
/* int det(x,y,z) */
/*      pcomplex *x, *y, *z; */
/* { */
/*   double det=0, */
/*     error=0, */
/*     sum; */
/*   int exp; */

/*   sum=x->x*y->y; */
/*   if (sum!=0){ */
/*     frexp(sum,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   det+=sum; */
/*   sum=x->y*z->x; */
/*   if (sum!=0){ */
/*     frexp(sum,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   det+=sum; */
/*   if (det!=0){ */
/*     frexp(det,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   sum=y->x*z->y; */
/*   if (sum!=0){ */
/*     frexp(sum,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   det+=sum; */
/*   if (det!=0){ */
/*     frexp(det,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   sum=-x->x*z->y; */
/*   if (sum!=0){ */
/*     frexp(sum,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   det+=sum; */
/*   if (det!=0){ */
/*     frexp(det,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   sum=-x->y*y->x; */
/*   if (sum!=0){ */
/*     frexp(sum,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   det+=sum; */
/*   if (det!=0){ */
/*     frexp(det,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   sum=-y->y*z->x; */
/*   if (sum!=0){ */
/*     frexp(sum,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   det+=sum; */
/*   if (det!=0){ */
/*     frexp(det,&exp); */
/*     error+=ldexp(1.,exp-DIGITS); */
/*   } */
/*   if (det>0){ */
/*     if (det<error) */
/*       return det5(x,y,z); */
/*     else  */
/*       return 1; */
/*   }   */
/*   if (det<0){ */
/*     if (det>-error) */
/*       return det5(x,y,z); */
/*     else  */
/*       return -1; */
/*   } */
/*   return det5(x,y,z); */
/* } */