//static void compute_wall_density_adami(particle* particles){
//  /*
//    Notes
//    -----
//    This is basically the implementation of Eq.(28) in the following paper:
//    "a generalized wall boundary condition for SPH"
//    by : adami et al. 2012
//   */
//
//  singleton_physics  *physics   = get_singleton_physics();
//  singleton_geometry *geometry  = get_singleton_geometry();
//
//  double X = physics->X;
//  double gamma_gas = physics->gamma_gas;
//  double gamma_liq = physics->gamma_liq;
//  double epsilon = 1e-12;
//
//  for (unsigned int b = 0; b < geometry->num_cell; b++) {
//      // find neighbors
//      int gi = 0;
//      int gj = 0;
//      grid_unhash(b,&gi,&gj);
//
//      int low_i  = gi-1 < 0 ? 0 : gi-1;
//      int low_j  = gj-1 < 0 ? 0 : gj-1;
//      int high_i = gi+2 > geometry->nx ? geometry->nx : gi+2;
//      int high_j = gj+2 > geometry->ny ? geometry->ny : gj+2;
//
//      for (unsigned int i = geometry->cell_start[b]; i < geometry->cell_end[b]; i++) {
//          if (!particles[i].bnd) continue;
//          double P0_fluid = particles[i].p0;
//          double rho0_fluid = particles[i].rho0;
//          double gamma = (particles[i].tag == AIR) ? gamma_gas : gamma_liq;
//
//          // loop over dummy neighbors
//          for (unsigned int ni = low_i; ni < high_i; ni++) {
//              for (unsigned int nj = low_j; nj < high_j; nj++) {
//                  for (unsigned int j = geometry->cell_start[ni*geometry->ny+nj]; j < geometry->cell_end[ni*geometry->ny+nj]; j++) {
//                      if (particles[j].bnd) continue;
//                      double Pwall = particles[j].p;
//
//                      // note: particle "i" from fluid phase is now interacting with particle "j" from dummy wall phase
//
//                      assert(fabs(P0_fluid) > epsilon);
//                      double frac = ((Pwall - X)/P0_fluid + 1.);
//                      assert(frac > epsilon);
//                      double rho_wall = rho0_fluid * pow(frac , (1./gamma));
//                      assert(rho_wall > epsilon);
//
//                      particles[i].rho = rho_wall;
//                  }
//              }
//          }
//      }
//  }
//}

static void compute_wall_density_adami(particle* particles){
    /*
    Notes
    -----
    This is basically the implementation of Eq.(28) in the following paper:
    "a generalized wall boundary condition for SPH"
    by : adami et al. 2012
     */

    singleton_physics  *physics   = get_singleton_physics();
    singleton_geometry *geometry  = get_singleton_geometry();

    double X = physics->X;
    double gamma_gas = physics->gamma_gas;
    double gamma_liq = physics->gamma_liq;
    double epsilon = 1e-12;

    // loop over dummy particles at wall only
    for (unsigned int i = 0; i < geometry->n; i++) {
        if (!particles[i].bnd) continue;

        // assign density/pressure of the interacting actual fluid particle
        double rho0_fluid = (particles[i].rho0==physics->rho_wat) ? physics->rho_wat : physics->rho_air;
        double P0_fluid = physics->p0L;
        double Pwall = particles[i].p;
        double gamma = (particles[i].tag == AIR) ? gamma_gas : gamma_liq;
        double frac = ((Pwall - X)/P0_fluid + 1.);
        assert(frac > epsilon);
        double rho_wall = rho0_fluid * pow(frac , (1./gamma));
        assert(rho_wall > epsilon);

        particles[i].rho = rho_wall;
    }
}