/* LB code for a simple ideal fluid in*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "lbschool.h"
#include "ran1.c"
#include "gasdev.c"
#include "nint.c"

static FLOAT pi=3.14159265358979323,ampl, wall_vel=0.0;

INT  max_x, max_y, max_z;
INT  num_x, num_y;
INT  siz_x, siz_y;

static INT x_dir[18] = {1, -1, 0, 0, 0, 0,
                           1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0};
static INT y_dir[18] = {0, 0, 1, -1, 0, 0,
                           1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1};
static INT z_dir[18] = {0, 0, 0, 0, 1, -1,
                           0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1};
static INT weighting[18] = {2, 2, 2, 2, 2, 2,
                           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
static INT l_opp[18] = {1, 0, 3, 2, 5, 4, 
                               7, 6, 9, 8, 11, 10, 13, 12, 15, 14, 17, 16};
static int  x_plus[5]  = {0, 6, 8, 10, 12};
static int  x_minus[5] = {1, 7, 9, 11, 13};
static int  y_plus[5]  = {2, 6, 9, 14, 16};
static int  y_minus[5] = {3, 7, 8, 15, 17};
static INT z_plus[5]  = {4,10,13,14,17};
static INT z_minus[5] = {5,11,12,15,16};

static long int seed= -16673;
static     FILE    *file_vcf;
static     FILE    *file_prof;
static     FILE    *file_temp;

int main()
{
  
  INT     num_step, n_step, link_dir;
  INT     x, y, z, velc_dir, index, index0,c0;
  
  FLOAT   *velcs_ptr[MAX_PTR][18];
  FLOAT   nu, rho0, m0_tot, vy0, vy,z1, kbt;
  char    buffer[MAX_BUFF];
  struct vector ptot, str;
  
  /* Set up here, for the time being */ 
  printf("init setup\n"); 
  rho0     = 24.0;
  nu       = 0.6;
  kbt      = 0.1;
  /* For simulation without noise set ampl  = 0.0 */
  ampl     = -sqrt(2.*rho0*nu*kbt)*(2.0/(6.0*nu-1.0))/8.0;
  nu       = (6.0*nu - 1.0)/(6.0*nu + 1.0);

# if DIM3 == 1
  max_x    = 20;
# else
  max_x    = 2;
# endif
  max_y    = 200;  
  max_z    = 200;
  num_x = max_x; num_y = max_y;  
  siz_x    = num_x + 2; siz_y = num_y + 2;
  num_step = 900;
  printf("IINIcio %d %d\n", num_step,max_x*max_y+max_y);
  /* Get the space and set the pointers */
  
  file_vcf=fopen("vcf.dat","w");
  file_prof=fopen("profile.dat","w");
  file_temp=fopen("temperature.dat","w");
  for( x = 1; x <= num_x; x++ ){
    for( y = 1;y<= num_y; y++){ 
      index = x*siz_y + y;
      for( velc_dir = 0; velc_dir < 18; velc_dir++){
	velcs_ptr[index][velc_dir] = 
	  (FLOAT *)malloc( max_z*sizeof( FLOAT ));
      }
    }  
  }  
   	
  lbe_initvel(velcs_ptr);
  /* init loop*/
  for(n_step=0;n_step<=num_step;n_step++){
    /* Compute vacf where momentum is applied */
    vy=0;
    index0=1*siz_y+max_y/2;
    for(velc_dir = 0; velc_dir < 5; velc_dir++){
      link_dir = y_plus[velc_dir];
      vy += *(velcs_ptr[index0][link_dir]+max_z/2);
      link_dir = y_minus[velc_dir];
      vy -= *(velcs_ptr[index0][link_dir]+max_z/2);
    }
    if(n_step==0)vy0=vy;
    fprintf(file_vcf,"%d %e %e\n",n_step,vy/vy0,vy*vy/rho0);
    printf("%d %e\n",n_step,vy);
    /* ENDOF Compute vacf where momentum is applied */

    /* For periodic boundary conditions comment these lines out */
    wallplus_bc(velcs_ptr);	  
    wallminus_bc(velcs_ptr);			  
    /* ENDOF For periodic boundary conditions comment these lines out */
    
    /* Apply peridic boundary conditions */
    lbe_bconds(velcs_ptr); 
    /* Advection step */   
    lbe_movexy(velcs_ptr);
    
    m0_tot = 0.0;        
    for(x=1; x <= num_x; x++){
      for(y=1; y <= num_y; y++){
	index = x*siz_y + y;
	/* Collision step */
	do_z_column( velcs_ptr[index], 
			       nu,
			       x,
			       y,
			       &ptot,
			       &str);
	m0_tot += (ptot.x+ptot.y+ptot.z)/3.0;
      }
    }
    /* Print velocity profiles */
    fprintf(file_temp,"%d %e\n",n_step,m0_tot/(FLOAT)(max_x*max_y*max_z*rho0));
    if((n_step%50)==0){
      for(z=1;z<max_z;z++){
	fprintf(file_prof,"%d %e\n",z,get1_momentum(velcs_ptr,3,&z));
      }
    }
  }
  fclose(file_vcf);
  fclose(file_prof);
  fclose(file_temp);
  return(1);
}
/*---------------------------------------------------------------------*/
void             lbe_initvel(  FLOAT   *velcs_ptr[][18])
/*---------------------------------------------------------------------*/
{

  INT      x, y, z, velc_dir, link_dir, index;  
  FLOAT    uy, vy;
  
  vy=0.0;
  /* set initial velocity to study vacf */
  uy = 0.0;
  for(x = 1; x <= num_x; x++){
    for(y = 1; y <= num_y; y++){
      index = x*siz_y + y;
      for(velc_dir = 0; velc_dir < 18; velc_dir++){
	for(z = 0; z < max_z; z++ ){ 
	  *(velcs_ptr[index][velc_dir]+z) = 1.0;                  
	}
      }
      for(velc_dir = 0; velc_dir < 5; velc_dir++){
	/* initial velocity profile in the y direction */
	for(z = 0; z < max_z; z++ ){ 
	  if((x==1)&&(z==max_z/2)&&(y==max_y/2)){
	    link_dir = y_plus[velc_dir];
	    *(velcs_ptr[index][link_dir]+z) += uy/10.0;   
	    vy +=*(velcs_ptr[index][link_dir]+z);
	    link_dir = y_minus[velc_dir];
	    *(velcs_ptr[index][link_dir]+z) -= uy/10.0;  
	    vy -= *(velcs_ptr[index][link_dir]+z);
	  }     
	} 
      }    
    }
  }  
  for(x = 1; x <= num_x; x++){
    for(y = 1; y <= num_y; y++){
      index = x*siz_y + y;
      for(velc_dir = 0; velc_dir < 6; velc_dir++){
	for(z = 0; z < max_z; z++ ){ 
	  *(velcs_ptr[index][velc_dir]+z) = 
	    *(velcs_ptr[index][velc_dir]+z)*2.0;
	}
      }
    }
  }
}
/*---------------------------------------------------------------------*/
void             lbe_bconds(  FLOAT   *velcs_ptr[][18] )
/*---------------------------------------------------------------------*/
{
  INT   x, y, xy_new, xy_old, link_dir, velc_dir;
  
  for(velc_dir = 0; velc_dir < 5; velc_dir ++){
    link_dir = x_plus[velc_dir];
    for( y = 1; y <= max_y; y++){
      xy_new = y;
      xy_old = max_x*siz_y + y;
      velcs_ptr[xy_new][link_dir] = velcs_ptr[xy_old][link_dir];
    }
    link_dir = x_minus[velc_dir];
    for( y = 1; y <= max_y; y++){
      xy_new = (max_x +1)*siz_y + y;
      xy_old = siz_y + y;
      velcs_ptr[xy_new][link_dir] = velcs_ptr[xy_old][link_dir];
    }
    link_dir = y_plus[velc_dir];
    for( x = 0; x < siz_x; x++){
      xy_new = x*siz_y;
      xy_old = x*siz_y + max_y;
      velcs_ptr[xy_new][link_dir] = velcs_ptr[xy_old][link_dir];
      
    }
    link_dir = y_minus[velc_dir];
    for( x = 0; x < siz_x; x++){
      xy_new = x*siz_y+max_y+1;
      xy_old = x*siz_y + 1;
      velcs_ptr[xy_new][link_dir] = velcs_ptr[xy_old][link_dir];
    }
  }
}
/*---------------------------------------------------------------------*/
void             lbe_movexy(  FLOAT   *velcs_ptr[][18] )
/*---------------------------------------------------------------------*/
{
  INT x, y, xy_new, xy_old, link_dir, velc_dir;
  
  for(velc_dir = 0; velc_dir < 5; velc_dir++){
    link_dir = x_plus[velc_dir];
    for(x = max_x; x >= 1; x--){
      for(y = 0; y < siz_y; y++){
	xy_new = x*siz_y + y;
	xy_old = xy_new - siz_y;
	velcs_ptr[xy_new][link_dir] =
	  velcs_ptr[xy_old][link_dir];
      }
    }
    link_dir = x_minus[velc_dir];
    for(x = 1; x <= max_x ; x++){
      for(y = 0; y < siz_y; y++){
	xy_new = x*siz_y + y;
	xy_old = xy_new + siz_y;
	velcs_ptr[xy_new][link_dir] =
	  velcs_ptr[xy_old][link_dir];
             }
    }
    
    link_dir = y_plus[velc_dir];
    for(x = 1; x <= max_x; x++){
      for(y = max_y; y >= 1; y--){
	xy_new = x*siz_y + y;
	xy_old = xy_new - 1;
	velcs_ptr[xy_new][link_dir] = velcs_ptr[xy_old][link_dir];
      }
    }
    link_dir = y_minus[velc_dir];
    for(x = 1; x <= max_x; x++){
      for(y = 1; y <= max_y; y++){
	xy_new = x*siz_y + y;
	xy_old = xy_new + 1;
	velcs_ptr[xy_new][link_dir] =
	  velcs_ptr[xy_old][link_dir];
      }
    }
  }
}
/*--------------------------------------------------------------------*/
FLOAT           do_z_column(FLOAT   *velcz_ptr[18], 
                            FLOAT   nu,
                            INT     x,
                            INT     y,
			    struct vector  *ptot,
                            struct vector  *str)
/*---------------------------------------------------------------------*/	{
  static INT     num_mode, initialize_chk=0;
  static FLOAT   nu_factor, *modes_ptr;
  FLOAT   m, m0, m1, m2, m3, m4, m5, m6, m7, m8, m9, rho_inverse;
  FLOAT   m09_a, m09_b, m45_a, m45_b, m42, m52;
  FLOAT   p0_tot,ran_1, ran_2, ran_3, ran_4, ran_5, ran_6, sum1, sum2, sum3;
  INT     z, z_plus, z_minus, z_ptr;

  if(initialize_chk == 0){
    num_mode = 10;
    nu_factor = nu;
    modes_ptr = (FLOAT *) malloc( max_z*num_mode*sizeof(FLOAT));
    initialize_chk = 1;
  }
  p0_tot = 0.0;
  (*ptot).x = (*ptot).y=(*ptot).z =0.0;
  for(z = 0, z_ptr = 0; z < max_z; z++, z_ptr += num_mode){
    z_plus  = mod( z+1, max_z );
    z_minus = mod( z-1+max_z, max_z );

    m  = *(velcz_ptr[0]+z);
    m0 = m; m1 = m; m4 = m; m5 = m;
    m  = *(velcz_ptr[1]+z);
    m0 += m; m1 -= m; m4 += m; m5 += m;
    m  = *(velcz_ptr[2]+z);
    m0 += m; m2 = m; m5 -= m;
    m  = *(velcz_ptr[3]+z);
    m0 += m; m2 -= m; m5 -= m;
    m  = *(velcz_ptr[4]+z_minus);
    m0 += m; m3 = m; m4 -= m;
    m  = *(velcz_ptr[5]+z_plus);
    m0 += m; m3 -= m; m4 -= m;
    m  = *(velcz_ptr[6]+z);
    m0 += m; m1 += m; m2 += m; m4 += m; m6  = m;
    m  = *(velcz_ptr[7]+z);
    m0 += m; m1 -= m; m2 -= m; m4 += m; m6 += m;
    m  = *(velcz_ptr[8]+z);
    m0 += m; m1 += m; m2 -= m; m4 += m; m6 -= m;
    m  = *(velcz_ptr[9]+z);
    m0 += m; m1 -= m; m2 += m; m4 += m; m6 -= m;
    m  = *(velcz_ptr[10]+z_minus);
    m0 += m; m1 += m; m3 += m; m5 += m; m7  = m;
    m  = *(velcz_ptr[11]+z_plus);
    m0 += m; m1 -= m; m3 -= m; m5 += m; m7 += m;
    m  = *(velcz_ptr[12]+z_plus);
    m0 += m; m1 += m; m3 -= m; m5 += m; m7 -= m;
    m  = *(velcz_ptr[13]+z_minus);
    m0 += m; m1 -= m; m3 += m; m5 += m; m7 -= m;
    m  = *(velcz_ptr[14]+z_minus);
    m0 += m; m2 += m; m3 += m; m4 -= m; m5 -= m; m8  = m;
    m  = *(velcz_ptr[15]+z_plus);
    m0 += m; m2 -= m; m3 -= m; m4 -= m; m5 -= m; m8 += m;
    m  = *(velcz_ptr[16]+z_plus);
    m0 += m; m2 += m; m3 -= m; m4 -= m; m5 -= m; m8 -= m;
    m  = *(velcz_ptr[17]+z_minus);
    m0 += m; m2 -= m; m3 += m; m4 -= m; m5 -= m; m8 -= m;
    m4 += m5; m5 = m4 - m5*2;
    m4 *= nu; m5 *= nu; m6 *= nu;
    m7 *= nu; m8 *= nu; m9  = 0.0;
#       if FLUID_TYPE == 2
    rho_inverse = 1.0/(m0);
        
    m4 += (2.0*m1*m1 - m2*m2 - m3*m3)*rho_inverse*(1.0 - nu_factor); 
    m5 += (            m2*m2 - m3*m3)*rho_inverse*(1.0 - nu_factor);
    m6 += (m1*m2                    )*rho_inverse*(1.0 - nu_factor);
    m7 += (m1*m3                    )*rho_inverse*(1.0 - nu_factor);
    m8 += (m2*m3                    )*rho_inverse*(1.0 - nu_factor);
    m9 += (    m1*m1 + m2*m2 + m3*m3)*rho_inverse*2.0;
#       endif 

    *(modes_ptr + z_ptr    ) = m0;
    *(modes_ptr + z_ptr + 1) = m1;
    *(modes_ptr + z_ptr + 2) = m2;
    *(modes_ptr + z_ptr + 3) = m3;
    
    p0_tot +=    m2;
    (*ptot).x += m1*m1;
    (*ptot).y += m2*m2;
    (*ptot).z += m3*m3;
    
    *(modes_ptr + z_ptr + 4) = m4;
    *(modes_ptr + z_ptr + 5) = m5;
    *(modes_ptr + z_ptr + 6) = m6;
    *(modes_ptr + z_ptr + 7) = m7;
    *(modes_ptr + z_ptr + 8) = m8;
    *(modes_ptr + z_ptr + 9) = m9;
  }  /* End loop over z column */

  for( z = 0, z_ptr = 0; z < max_z; z++, z_ptr += num_mode ){
    m0 = *(modes_ptr + z_ptr    )/24.0;
    m1 = *(modes_ptr + z_ptr + 1)/12.0;
    m2 = *(modes_ptr + z_ptr + 2)/12.0;
    m3 = *(modes_ptr + z_ptr + 3)/12.0;
    m4 = *(modes_ptr + z_ptr + 4)/48.0;
    m5 = *(modes_ptr + z_ptr + 5)/16.0;
    m6 = *(modes_ptr + z_ptr + 6)/ 4.0;
    m7 = *(modes_ptr + z_ptr + 7)/ 4.0;
    m8 = *(modes_ptr + z_ptr + 8)/ 4.0;
    m9 = *(modes_ptr + z_ptr + 9)/24.0;

    m09_a = (m0 - m9)*2.0; m09_b = m0 + m9;
    m45_a = m4 + m5; m45_b = m4 - m5;
    m42   = 2.0*m4; m52 = 2.0*m5;

    /* here moise implementation */    
    ran_1=gasdev(&seed);
    ran_2=gasdev(&seed);
    ran_3=gasdev(&seed);
    ran_4=gasdev(&seed);
    ran_5=gasdev(&seed);
    ran_6=gasdev(&seed);

    sum1 = ( 2.0*ran_1-    ran_2-    ran_3)/3.0;
    sum2 = (-    ran_1+2.0*ran_2-    ran_3)/3.0;
    sum3 = (-    ran_1-    ran_2+2.0*ran_3)/3.0;

    *(velcz_ptr[0]+z) = m09_a + m1*2.0 + m4*4.0             +ampl*sum1*2.0;
    *(velcz_ptr[1]+z) = m09_a - m1*2.0 + m4*4.0             +ampl*sum1*2.0;
    *(velcz_ptr[2]+z) = m09_a + m2*2.0 - m42    +m52        +ampl*sum2*2.0;
    *(velcz_ptr[3]+z) = m09_a - m2*2.0 - m42    +m52        +ampl*sum2*2.0;
    *(velcz_ptr[4]+z) = m09_a + m3*2.0 - m42    -m52        +ampl*sum3*2.0;
    *(velcz_ptr[5]+z) = m09_a - m3*2.0 - m42    -m52        +ampl*sum3*2.0;
    *(velcz_ptr[6]+z) = m09_b + m1     + m2     +m45_a + m6 +ampl*( 2*ran_4-sum3);
    *(velcz_ptr[7]+z) = m09_b - m1     - m2     +m45_a + m6 +ampl*( 2*ran_4-sum3);
    *(velcz_ptr[8]+z) = m09_b + m1     - m2     +m45_a - m6 +ampl*(-2*ran_4-sum3);
    *(velcz_ptr[9]+z) = m09_b - m1     + m2     +m45_a - m6 +ampl*(-2*ran_4-sum3);
    *(velcz_ptr[10]+z)= m09_b + m1     + m3     +m45_b + m7 +ampl*( 2*ran_5-sum2);
    *(velcz_ptr[11]+z)= m09_b - m1     - m3     +m45_b + m7 +ampl*( 2*ran_5-sum2);
    *(velcz_ptr[12]+z)= m09_b + m1     - m3     +m45_b - m7 +ampl*(-2*ran_5-sum2);
    *(velcz_ptr[13]+z)= m09_b - m1     + m3     +m45_b - m7 +ampl*(-2*ran_5-sum2);
    *(velcz_ptr[14]+z)= m09_b + m2     + m3     -m42   + m8 +ampl*( 2*ran_6-sum1);
    *(velcz_ptr[15]+z)= m09_b - m2     - m3     -m42   + m8 +ampl*( 2*ran_6-sum1);
    *(velcz_ptr[16]+z)= m09_b + m2     - m3     -m42   - m8 +ampl*(-2*ran_6-sum1);
    *(velcz_ptr[17]+z)= m09_b - m2     + m3     -m42   - m8 +ampl*(-2*ran_6-sum1);
  } /* End loop over column */
  return(p0_tot);   
}
/*---------------------------------------------------------------------*/
FLOAT       get1_momentum(  FLOAT      *velcs_ptr[][18], INT dir, INT *ival )
     /*---------------------------------------------------------------------*/
{
  INT x, y, z, index, l;
  FLOAT py_tot,size;
  
  py_tot = 0.0;
  if(dir==1){
    x=(*ival);
    size=(FLOAT)(max_y*max_z);
    for( y = 1; y <= max_y; y++ ){
      index = x*siz_y + y;
      for( z = 0; z < max_z; z++ ){
	for( l = 0; l < 18; l++ ){
	  py_tot += *(velcs_ptr[index][l] + z )*y_dir[l];
	}
      }
    }
  }else{
    if(dir==2){
      size=(FLOAT)(max_x*max_z);
      y=(*ival);
      for( x = 1; x <= max_x; x++){
	index = x*siz_y + y;
	for( z = 0; z < max_z; z++ ){
	  for( l = 0; l < 18; l++ ){
	    py_tot += *(velcs_ptr[index][l] + z )*y_dir[l];
	  }
	}
      }
    }else{
      size=(FLOAT)(max_y*max_x);
      z=(*ival);
      for( x = 1; x <= max_x; x++){
	for( y = 1; y <= max_y; y++ ){
	  index = x*siz_y + y;
	  for( l = 0; l < 18; l++ ){
	    py_tot += *(velcs_ptr[index][l] + z )*y_dir[l];
	  }
	}
      }
    }
  }
  return(py_tot/size);
}
/*---------------------------------------------------------------------*/
void      wallplus_bc(  FLOAT      *velcs_ptr[][18] )
     /*---------------------------------------------------------------------*/
{
  INT x,y,lk,l1,l2,x1,y1,z1,x2,y2,z2,xy1,xy2;
  FLOAT f1,f2,uw,ampl;
  
  uw=wall_vel;
  for (x = 1; x <= num_x; x++ ){
    for (y = 1; y <= num_y; y++ ){
      x1 = x;
      y1 = y;
      z1 = max_z-1;
      for (lk = 0; lk < 5; lk++ ){
	l1=z_minus[lk];
	x2     = x1 + x_dir[l1];
	y2     = y1 + y_dir[l1];
	z2     = z1 + z_dir[l1]; z2 = mod(z2 + max_z, max_z);
	
	x2 = mod(x2 + max_x -1, max_x) + 1;
	y2 = mod(y2 + max_y -1, max_y) + 1;
	
	l2     = l_opp[l1];
	
	xy1    = x1*siz_y + y1;
	xy2    = x2*siz_y + y2;
	ampl=1.0;
	if(l1<6) ampl=2.0;
	f1 = *(velcs_ptr[xy1][l1] + z1)+ 2.0*(24./12.)*uw*y_dir[l1]*ampl;
	f2 = *(velcs_ptr[xy2][l2] + z2)- 2.0*(24./12.)*uw*y_dir[l2]*ampl;
	
	*(velcs_ptr[xy1][l1] + z1) = f2;
	*(velcs_ptr[xy2][l2] + z2) = f1;
      }
    }
  }
}
/*---------------------------------------------------------------------*/
void      wallminus_bc(  FLOAT      *velcs_ptr[][18] )
     /*---------------------------------------------------------------------*/
{
  INT x,y,lk,l1,l2,x1,y1,z1,x2,y2,z2,xy1,xy2;
  FLOAT f1,f2,uw,ampl;
  
  uw = -wall_vel;
  for (x = 1; x <= num_x; x++ ){
    for (y = 1; y <= num_y; y++ ){
      x1 = x;
      y1 = y;
      z1 = 0;
      for (lk = 0; lk < 5; lk++ ){
	l1=z_plus[lk];
	x2     = x1 + x_dir[l1];
	y2     = y1 + y_dir[l1];
	z2     = z1 + z_dir[l1]; z2 = mod(z2 + max_z, max_z);
	
	x2 = mod(x2 + max_x -1, max_x) + 1;
	y2 = mod(y2 + max_y -1, max_y) + 1;
	
	l2     = l_opp[l1];
	
	xy1    = x1*siz_y + y1;
	xy2    = x2*siz_y + y2;
	
	ampl=1.0;
	if(l1<6) ampl=2.0;
	f1 = *(velcs_ptr[xy1][l1] + z1)+ 2.0*(24./12.)*uw*y_dir[l1]*ampl;
	f2 = *(velcs_ptr[xy2][l2] + z2)- 2.0*(24./12.)*uw*y_dir[l2]*ampl;
	
	*(velcs_ptr[xy1][l1] + z1) = f2;
	*(velcs_ptr[xy2][l2] + z2) = f1;
      }
    }
  }
}
/*---------------------------------------------------------------------*/