//Line number in the x direction. Range is from 0 to I-1.
 int i_vr(int l){
 return l%I; 
}

//Line number in the y direction
 int j_vr(int l){
 return (l%(I*J))/I; 
}

//Line number in the z direction
 int k_vr(int l){
 return l/(I*J); 
}

//Characteristic velocity vectors of the LV model. Opposite vector pairs have to be written in odd-even combination (important for calcLBmnew function)!
void initialiseE()
{
  E[0][1]= 0;		E[0][2]= 0;		E[0][3]= 0;
  
  E[1][1]= 1;		E[1][2]= 0;		E[1][3]= 0;
  E[2][1]= -1;		E[2][2]= 0;		E[2][3]= 0;
  E[3][1]= 0;		E[3][2]= 1;		E[3][3]= 0;
  E[4][1]= 0; 		E[4][2]= -1;		E[4][3]= 0;
  E[5][1]= 0;		E[5][2]= 0;		E[5][3]= 1;
  E[6][1]= 0;  		E[6][2]= 0;  		E[6][3]= -1;
  
  E[7][1]= 1;  		E[7][2]= 1;	  	E[7][3]= 0;
  E[8][1]= -1; 		E[8][2]= -1;		E[8][3]= 0;
  E[9][1]= -1;  	E[9][2]= 1;		E[9][3]= 0;
  E[10][1]= 1;		E[10][2]= -1;		E[10][3]= 0;
  
  E[11][1]= 1;  	E[11][2]= 0;		E[11][3]= 1;
  E[12][1]= -1;  	E[12][2]= 0;  		E[12][3]= -1;
  E[13][1]= -1;  	E[13][2]= 0;		E[13][3]= 1;
  E[14][1]= 1;		E[14][2]= 0;		E[14][3]= -1;
  
  E[15][1]= 0;  	E[15][2]= 1;		E[15][3]= 1;
  E[16][1]= 0;  	E[16][2]=-1;  		E[16][3]= -1;
  E[17][1]= 0;  	E[17][2]= -1;		E[17][3]= 1;
  E[18][1]= 0;		E[18][2]= 1;		E[18][3]= -1;
  
  //Constants for the model. Adapted from Thurey, Rude, Korner
  int i;
  WKONST[0]=1.f/3.f;
  for(i=1;i<=6; i++) WKONST[i]=1.f/18.f;
  for(i=7;i<=18; i++) WKONST[i]=1.f/36.f;
}



//Function that multiplies matrices Aik x Bkj = Cij t and returns the element Cij
void AMatrixB(float (*c)[4], float (*a)[4],float (*b)[4])
{
  for (int i = 1; i<=3; i++)
    {
      for (int j = 1; j<=3; j++)
	{
	  c[i][j]=a[i][1]*b[1][j]+a[i][2]*b[2][j]+a[i][3]*b[3][j];
	}
    }
}

// Function to take a little-endian encoded float and make it big-endian (or vice versa). Useful for reading to paraview
float FloatSwap( float f )
{
    union
    {
        float f;
        char b[4];
    } dat1, dat2;

    dat1.f = f;
    dat2.b[0] = dat1.b[3];
    dat2.b[1] = dat1.b[2];
    dat2.b[2] = dat1.b[1];
    dat2.b[3] = dat1.b[0];
    return dat2.f;
}