/** # This file include mathematical operator for coord struct */
typedef struct ten{
  coord x;
  coord y;
  #if dimension == 3
  coord z;
  #endif
}ten;
/** declaration of the null coord*/

const coord Coord_nul = {0};
const tens tens_nul = {{0}};
#define Identity\
          {{1,0,0},\
          {0,1,0},\
          {0,0,1}}
const tens I = Identity;

/** and operator on coord */
coord and_coord(coord a,coord b){
  coord c = {(int)a.x && (int)b.x,(int)a.y &&(int) b.y,(int)a.z &&(int) b.z};
  return c;
}
coord add_coord(coord a,coord b){
  #if dimension == 2
  coord c = {a.x+b.x, a.y+b.y};
  #elif dimension == 3
  coord c = {a.x+b.x, a.y+b.y, a.z+b.z};
  #endif
  return c;
}
coord diff_coord(coord a,coord b){
  #if dimension == 2
  coord c = {a.x-b.x,a.y-b.y};
  #elif dimension == 3
  coord c = {a.x-b.x, a.y-b.y, a.z-b.z};
  #endif
  return c;
}
coord mult_coord(coord a,double b){
  #if dimension == 2
  coord c = {a.x*b,  a.y*b};
  #elif dimension == 3
  coord c = {a.x*b,  a.y*b,  a.z*b};
  #endif
  return c;
} 
coord div_coord(coord a,double b){
  #if dimension == 2
  coord c = {a.x/b,  a.y/b};
  #elif dimension == 3
  coord c = {a.x/b,  a.y/b,  a.z/b};
  #endif
  return c;
}
tens div_tens(tens a,double b){
  einstein_sum(i,j){
     a.i.j /= b;
  }
  return a;
}
tens mult_tens(tens a,double b){
  einstein_sum(i,j){
     a.i.j *= b;
  }
  return a;
}
tens add_tens(tens A,tens B){
  einstein_sum(i,j){
     A.i.j += B.i.j;
  }
  return A;
}
double tens_norm(tens A){
  double norm = 0;
  einstein_sum(k,l){
     norm = A.k.l * A.k.l;
  }
  norm = sqrt(norm);
  return norm;
}
double normL2_coord(coord a){
  #if dimension == 2
  double c = sqrt(pow(a.x,2)+pow(a.y,2));
  #elif dimension == 3
  double c = sqrt(pow(a.x,2)+pow(a.y,2)+pow(a.z,2));
  #endif
  return c;
}
double normL2_coord_sq(coord a){
  #if dimension == 2
  double c = sq(a.x)+sq(a.y);
  #elif dimension == 3
  double c = sq(a.x)+sq(a.y)+sq(a.z);
  #endif
  return c;
}
double normL2_ten(ten a){
  #if dimension == 2
  double c = sqrt(pow(a.x.x,2)+pow(a.y.y,2)+pow(a.x.y,2)+pow(a.y.x,2));
  #elif dimension == 3
  double c = sqrt(pow(a.x.x,2)+pow(a.y.y,2)+pow(a.z.z,2)+pow(a.x.y,2)+pow(a.y.x,2)+pow(a.x.z,2)+pow(a.z.x,2)+pow(a.z.y,2)+pow(a.y.z,2));
  #endif
  return c;
}
//only in 3D for the cross product 
#if dimension == 3
coord cross_coord(coord a,coord b){
  coord c = {a.y*b.z-a.z*b.y, 
            a.z*b.x-a.x*b.z, 
            a.x*b.y-a.y*b.x};
#elif dimension == 2
double cross_coord(coord a,coord b){
  double c = a.x*b.y-a.y*b.x;
#endif
  return c;
}

coord EigenValue(ten A){
  double Delta = pow(A.x.x,2.) - 2.*A.x.x*A.y.y+pow(A.y.y,2.) +4.*A.x.y*A.y.x;
  coord Eig = {0.,0.};
  if(Delta > 0){
    Eig.x = 1./2.*(  A.x.x + A.y.y - sqrt( Delta ));
    Eig.y = 1./2.*(  A.x.x + A.y.y + sqrt( Delta ));
  }
  return Eig;
}

ten EigenVectors(ten A){
  ten V;
  if(A.y.x != 0){
    V.x.x = -1./(2.*A.y.x)*( - A.x.x + A.y.y + sqrt( pow(A.x.x,2.) - 2.*A.x.x*A.y.y+pow(A.y.y,2.) +4.*A.x.y*A.y.x ) ); 
    V.x.y = 1.;
    V.y.x = -1./(2.*A.y.x)*( - A.x.x + A.y.y - sqrt( pow(A.x.x,2.) - 2*A.x.x*A.y.y+pow(A.y.y,2.) +4.*A.x.y*A.y.x ) );
    V.y.y = 1.;
  }else if((A.x.x-A.y.y)!=0){
    V.x.x = 1.;
    V.x.y = 0.;
    V.y.x = -A.x.y/(A.x.x-A.y.y);
    V.y.y = 1.;
  }else{
    V.x.x = 1.;
    V.x.y = 0.;
    V.y.x = 0.;
    V.y.y = 1.;
  }
  V.x = div_coord(V.x,normL2_coord(V.x));
  V.y = div_coord(V.y,normL2_coord(V.y));
  return V;
}

/** Set the coordinate of a coord at the positive side of the domain */
coord POS_perio(coord pos,coord per){
  foreach_dimension(){
    if(per.x) pos.x = pos.x>0?pos.x:pos.x+Ls;
  }
  return pos;
}
/** compute the periodic distance between two coord */
coord dist_perio_coord(coord a,coord b){
  foreach_dimension(){
    if(a.x>0) b.x = (b.x<a.x-Ls/2)*Ls+b.x;
    else b.x = (b.x>a.x+Ls/2)*(-Ls)+b.x;
  }
  return diff_coord(a,b);
}


double dist_perio(coord a,coord b){
  foreach_dimension(){
    if(a.x>0) b.x = (b.x<a.x-Ls/2)*Ls+b.x;
    else b.x = (b.x>a.x+Ls/2)*(-Ls)+b.x;
  }
  return normL2_coord(diff_coord(a,b));
}
/** compute the tranpsoe of a tensor */ 
tens transpsoe(tens A){
  tens B;
  einstein_sum(i,j){
     B.i.j = A.j.i;
  }
  return B;
}
tens sym(tens A){
  tens B;
  einstein_sum(i,j){
  B. i.j = 0.5*(A.i.j + A.j.i);
  }
  return B;
}
tens coord_prod(coord A,coord B){
  tens C;
  einstein_sum(i,j){
     C.i.j = A.i*B.j;
  }
  return C;
}
/** compute the symetric part of a tensor */