/*
What's new:
C44 - [Martin Zoubek] first public release
*/

#include <math.h>
#include "Triangle.h"
#include "Types.h"

/* void *att_init(FILE *,int,int);  */

double Edge_Length(const double2 p1,const double2 p2){
  double2 A;
  double2_sub(A,p2,p1);
  return sqrt(double2_norm2(A)); 
  }

double Triangle_Volume(const double2 p1,const double2 p2,const double2 p3){
  double2 A,B;
  double c;
  if(!p1||!p2||!p3)
    return -1;
  double2_sub(A,p2,p1);
  double2_sub(B,p3,p1);
  double2_crossprod(c,A,B);
  return 0.5*(c);
  /*  return 0.5*sqrt(double2_norm2(C)); */
  }

double Triangle_Deter(const double2 p1,const double2 p2,const double2 p3){
  double c;
  if(!p1||!p2||!p3)
    return -1;
  c = p1[0]*(p2[1]-p3[1]) + p2[0]*(p3[1]-p1[1]) + p3[0]*(p1[1]-p2[1]);
  return 0.5*(c);
}


void Triangle_NormalVector(const double2 p1,const double2 p2,const double2 p3,
			   double2 n){
  double2 A,B,C;
  double v;
  if(!p1||!p2||!p3||!n)
    return;
  double2_sub(A,p2,p1);
  /*double2_sub(B,p3,p1);
  double2_crossprod(n,A,B);
  v=double2_norm2(n);
  if(v<1e-40){
    n[0]=n[1]=n[2]=0;
    return;
    }
  else
    v=1/sqrt(v);
  if(p4){
    double3_sub(C,p4,p1);
    if(double3_dotprod(n,C)>0)
      v*=-1;
    }
  n[0]*=v;
  n[1]*=v;
  n[2]*=v;*/

  n[0] = A[1];
  n[1] = -A[0];
  v=double2_norm2(n);
  if(v<1e-40){
    n[0]=n[1]=n[2]=0;
    printf("Zero outer normal\n");
    return;
    }
  else
    v=1/sqrt(v);
  n[0]*=v;
  n[1]*=v;

  }

void Triangle_Gravicity(const double2 p1,const double2 p2,const double2 p3,
    double2 n){

  if(!p1||!p2||!p3||!n)
    return;
  double2_gravicity(n,p1,p2,p3);
}

/* tSurface_Surface *Surface_Init(FILE *f,uint method,uint maxiter){
  return att_init(f,maxiter,method);
  } */

void Triangle_VectorRotation(const double2 A1,const double2 B1,
  const double2 C1,const double2 A2,const double2 B2,
  const double2 C2,const double2 v1,double2 v2){
  double2 a1,b1,a2,b2,n1,n2,c1,c2;
  double c,ab,sign=-1,f,g,ca,cb,cn;
  uint i;
  if(!A1||!B1||!C1||!A2||!B2||!C2||!v1||!v2)
    return;
  double2_sub(a1,A1,C1);
  double2_sub(b1,B1,C1);
  double2_sub(a2,A2,C2);
  double2_sub(b2,B2,C2);
  
  /*  double2_crossprod(n1,a1,b1);  //normal vectors
  double2_crossprod(n2,a2,b2); */
  /*  double2_sub(c1,D1,C1);  //need outer normal vectors
  double2_sub(c2,D2,C2); */
  if(double2_dotprod(n1,c1)*double2_dotprod(n2,c2)>0)
    sign=1;
  c=1/double2_norm2(n1);
  ab=double2_dotprod(a1,b1);
  f=double2_dotprod(a1,v1);
  g=double2_dotprod(b1,v1);
  ca=double2_norm2(b1)*f-ab*g;
  cb=ab*f-double2_norm2(a1)*g;
  cn=sign*double2_dotprod(n1,v1);
  for(i=0;i<2;i++)
    v2[i]=c*(ca*a2[i]-cb*b2[i]-cn*n2[i]);
  }

