// Requires pcomplex
#include "math.h"

typedef struct
{
pcomplex z[3];
}
triangle;

triangle_scale(T, scale)
     triangle *T;
     double scale;
{
  int i;
  for (i=0; i<3; i++){
    T->z[i].x*=scale;
    T->z[i].y*=scale;
  }
}

pcomplex triangle_center(T)
     triangle *T;
{
  pcomplex z;
  z.x=(T->z[0].x+T->z[1].x+T->z[2].x)/3;
  z.y=(T->z[0].y+T->z[1].y+T->z[2].y)/3;
  return z;
}

triangle triangle_recenter(T)
     triangle *T;
{
  pcomplex z;
  int i;
  triangle ret;
  z=triangle_center(T);
  for (i=0; i<3; i++){
    ret.z[i].x=T->z[i].x-z.x;
    ret.z[i].y=T->z[i].y-z.y;
  }
  return ret;
}

void triangle_recenter2(T)
     triangle *T;
{
  pcomplex z;
  int i;
  z=triangle_center(T);
  for (i=0; i<3; i++){
    T->z[i].x-=z.x;
    T->z[i].y-=z.y;
  }
}

double triangle_radius(T)
     triangle *T;
{
  pcomplex z;  
  double a,b;
  z=triangle_center(T);
  a=sqrt((T->z[0].x-z.x)*(T->z[0].x-z.x)+
	 (T->z[0].y-z.y)*(T->z[0].y-z.y));
  b=sqrt((T->z[1].x-z.x)*(T->z[1].x-z.x)+
	 (T->z[1].y-z.y)*(T->z[1].y-z.y));
  if (b>a)
    a=b;
   b=sqrt((T->z[2].x-z.x)*(T->z[2].x-z.x)+
	  (T->z[2].y-z.y)*(T->z[2].y-z.y));
  if (b>a)
    a=b; 
  return a;
}

// a should be a triangle[4]
void triangle_subdivide(T, a)
     triangle *T, a[];
{
  int i,j;
  for (i=0; i<3;i++){
    for (j=0;j<3;j++){
      a[i].z[j].x=(T->z[i].x+T->z[j].x)/2;
      a[i].z[j].y=(T->z[i].y+T->z[j].y)/2;
    }
  }
  for (i=0; i<3;i++){
    a[3].z[i].x=(T->z[(i+1)%3].x+T->z[(i+2)%3].x)/2;
    a[3].z[i].y=(T->z[(i+1)%3].y+T->z[(i+2)%3].y)/2;
  }
}

double triangle_area(T)
     triangle *T;
{
  return 
    T->z[0].x*T->z[1].y-T->z[0].y*T->z[1].x+
    T->z[1].x*T->z[2].y-T->z[1].y*T->z[2].x+
    T->z[2].x*T->z[0].y-T->z[2].y*T->z[0].x;    
}

// Triangle in Rich's format
triangle triangle_make(z)
     pcomplex z;
{
  double x,y;
  triangle T;

  x=Pi*z.x/2.0;
  y=Pi*z.y/2.0;
  T.z[0].x=0.0;
  T.z[0].y=0.0;
  T.z[1].x=0.0;
  T.z[1].y=1.0;
  T.z[2].x=(tan(x)*tan(y))/(tan(x)+tan(y));
  T.z[2].y=tan(y)/(tan(x)+tan(y));
  return(T);
}