text: #ifndef traverse_eager_h
  #define traverse_eager_h
  #include "exafmm.h"
  #include "kernel.h"
  
  namespace exafmm {
    //! Recursive call to post-order tree traversal for upward pass
    void upwardPass(Cell * Ci) {
      for (Cell * Cj=Ci->child; Cj!=Ci->child+Ci->numChilds; Cj++) {
  #pragma omp task untied if(Cj->numBodies > 100)
        upwardPass(Cj);
      }
  #pragma omp taskwait
      if(Ci->numChilds==0) P2M(Ci);
      M2M(Ci);
    }
  
    //! Upward pass interface
    void upwardPass(Cells & cells) {
  #pragma omp parallel
  #pragma omp single nowait
      upwardPass(&cells[0]);
    }
  
    //! Upward pass to fill in missing numBodies and M
    void upwardPassLET(Cell * Ci) {
      for (Cell * Cj=Ci->child; Cj!=Ci->child+Ci->numChilds; Cj++) {
  #pragma omp task untied if(Cj->numBodies > 100)
        upwardPassLET(Cj);
      }
  #pragma omp taskwait
      real_t M = 0;
      for (int n=0; n<NTERM; n++) M += std::abs(Ci->M[n]);
      if (Ci->numChilds==0) {
        if (M < EPS) {
          P2M(Ci);
        }
      } else {
        for (int n=0; n<NTERM; n++) Ci->M[n] = 0;
        M2M(Ci);
      }
    }
  
    //! Upward pass LET interface
    void upwardPassLET(Cells & cells) {
  #pragma omp parallel
  #pragma omp single nowait
      upwardPassLET(&cells[0]);
    }
  
    //! Recursive call to dual tree traversal for horizontal pass
    void horizontalPass(Cell * Ci, Cell * Cj) {
      vec3 dX;
      for (int d=0; d<3; d++) dX[d] = Ci->X[d] - Cj->X[d] - IX[d] * CYCLE;
      real_t R2 = norm(dX) * THETA * THETA;
      if (R2 > (Ci->R + Cj->R) * (Ci->R + Cj->R)) {
        M2L(Ci, Cj);
      } else if (Ci->numChilds == 0 && Cj->numChilds == 0) {
        if (Cj->numBodies == 0) {
          assert((Ci->R+Cj->R)/std::sqrt(norm(dX)) < THETA*1.1);
          M2L(Ci, Cj);
        } else {
          P2P(Ci, Cj);
        }
      } else if (Cj->numChilds == 0 || (Ci->R >= Cj->R && Ci->numChilds != 0)) {
        for (Cell * ci=Ci->child; ci!=Ci->child+Ci->numChilds; ci++) {
  #pragma omp task untied if(ci->numBodies > 100) firstprivate(IX)
          horizontalPass(ci, Cj);
        }
      } else {
        for (Cell * cj=Cj->child; cj!=Cj->child+Cj->numChilds; cj++) {
          horizontalPass(Ci, cj);
        }
      }
  #pragma omp taskwait
    }
  
    //! Horizontal pass for periodic images
    void periodic(Cell * Ci0, Cell * Cj0) {
      Cells pcells(27);
      for (size_t c=0; c<pcells.size(); c++) {
        pcells[c].M.resize(NTERM, 0.0);
        pcells[c].L.resize(NTERM, 0.0);
      }
      Cell * Ci = &pcells.back();
      *Ci = *Cj0;
      Ci->child = &pcells[0];
      Ci->numChilds = 26;
      for (int level=0; level<IMAGES-1; level++) {
        for (int ix=-1; ix<=1; ix++) {
          for (int iy=-1; iy<=1; iy++) {
            for (int iz=-1; iz<=1; iz++) {
              if (ix != 0 || iy != 0 || iz != 0) {
                for (int cx=-1; cx<=1; cx++) {
                  for (int cy=-1; cy<=1; cy++) {
                    for (int cz=-1; cz<=1; cz++) {
                      IX[0] = ix * 3 + cx;
                      IX[1] = iy * 3 + cy;
                      IX[2] = iz * 3 + cz;
                      M2L(Ci0, Ci);
                    }
                  }
                }
              }
            }
          }
        }
        Cell * Cj = &pcells[0];
        for (int ix=-1; ix<=1; ix++) {
          for (int iy=-1; iy<=1; iy++) {
            for (int iz=-1; iz<=1; iz++) {
              if (ix != 0 || iy != 0 || iz != 0) {
                Cj->X[0] = Ci->X[0] + ix * CYCLE;
                Cj->X[1] = Ci->X[1] + iy * CYCLE;
                Cj->X[2] = Ci->X[2] + iz * CYCLE;
                Cj->M = Ci->M;
                Cj++;
              }
            }
          }
        }
        M2M(Ci);
        CYCLE *= 3;
      }
    }
  
    //! Horizontal pass interface
    void horizontalPass(Cells & icells, Cells & jcells) {
  #pragma omp parallel
  #pragma omp single
      if (IMAGES == 0) {
        horizontalPass(&icells[0], &jcells[0]);
      } else {
        for (IX[0]=-1; IX[0]<=1; IX[0]++) {
          for (IX[1]=-1; IX[1]<=1; IX[1]++) {
            for (IX[2]=-1; IX[2]<=1; IX[2]++) {
              horizontalPass(&icells[0], &jcells[0]);
            }
          }
        }
        real_t saveCycle = CYCLE;
        periodic(&icells[0], &jcells[0]);
        CYCLE = saveCycle;
      }
    }
  
    //! Recursive call to pre-order tree traversal for downward pass
    void downwardPass(Cell * Cj) {
      L2L(Cj);
      if (Cj->numChilds==0) L2P(Cj);
      for (Cell * Ci=Cj->child; Ci!=Cj->child+Cj->numChilds; Ci++) {
  #pragma omp task untied if(Ci->numBodies > 100)
        downwardPass(Ci);
      }
  #pragma omp taskwait
    }
  
    //! Downward pass interface
    void downwardPass(Cells & cells) {
  #pragma omp parallel
  #pragma omp single nowait
      downwardPass(&cells[0]);
    }
  
    //! Direct summation
    void direct(Bodies & bodies, Bodies & jbodies) {
      Cells cells(2);
      Cell * Ci = &cells[0];
      Cell * Cj = &cells[1];
      Ci->body = &bodies[0];
      Ci->numBodies = bodies.size();
      Cj->body = &jbodies[0];
      Cj->numBodies = jbodies.size();
      int prange = (std::pow(3,IMAGES) - 1) / 2;
      for (int ix=-prange; ix<=prange; ix++) {
        for (int iy=-prange; iy<=prange; iy++) {
          for (int iz=-prange; iz<=prange; iz++) {
            IX[0] = ix;
            IX[1] = iy;
            IX[2] = iz;
            P2P(Ci, Cj);
          }
        }
      }
    }
  }
  #endif
  