// timestep.c
double rk1(double u1, double u2, double u3, double u4, double u5, double Dt, double Dx, double c)
{
  double D = c*( flux(u2,u3,u4,c) - flux(u1,u2,u3,c) )/Dx; 
  return  u3 - Dt*D;
}


double rk1nd(int d, int prof, VFl pc[], double Dt, double ex, double xmin[], double xmax[], Fluxo f)
{
  double D = 0.0;
  int dir; 
  for(dir = 0; dir < d; dir++)
  {
    double dx = fabs(xmax[dir] - xmin[dir]);
    double Dx = /*ex == dir ? 0.5*(dx + 0.5*dx) :*/ dx;
    int k = dir == 0 ? dir : dir + 1;
    if((pc[k].p == NULL) || (pc[k].p->fil[0] == NULL))
    {
      D += (pc[1].p->flxsup[dir] - pc[k].fvi[dir])/Dx; 
      //D += (pc[1].p->flxsup[dir] - pc[1].p->flxinf[dir])/Dx;
    }
    else
    {
      assert( pc[k].p != NULL );
      double somafluxo = 0.0;
      double xml[d], xMl[d];
      int i;
      for(i=0;i<d;i++)
      {
        xml[i] = i == dir ? xmin[i] - (xmax[i]-xmin[i]) : xmin[i];
        xMl[i] = i == dir ? xmax[i] - (xmax[i]-xmin[i]) : xmax[i];
      }
      encontra_vizinhos_inferiores_esq_refinada( d, pc[k].p, prof, dir, xml, xMl, xmin, &somafluxo);
      D += (pc[1].p->flxsup[dir] - somafluxo)/Dx; 
      //D += (pc[1].p->flxsup[dir] - pc[1].p->flxinf[dir])/Dx;
    }
  }
  double s3 = pc[1].p->fval; /* solucao (nesta celula) antes do primeiro passo de rk */
  return s3 - Dt*D;       /* primeiro passo de rk */
}


double rk1_nd(int d, int prof, VNo pac[], double Dt, double xmin[], double xmax[], Fluxo f)
{
  double D = 0.0;
  int ex = prof%d;
  int ord = 5;
  int ctr = pow(ord,d)/2;
  int dir; 
  for(dir = 0; dir < d; dir++)
  {
    int corr = ipow(ord,dir);
    double dx = fabs(xmax[dir] - xmin[dir]);
    double Dx = /*ex == dir ? 0.5*(dx + 0.5*dx) :*/ dx;
    double u1 = pac[ctr-2*corr].p == NULL ? pac[ctr-2*corr].fv : pac[ctr - 2*corr].p->fval;
    double u2 = pac[ctr-1*corr].p == NULL ? pac[ctr-1*corr].fv : pac[ctr - 1*corr].p->fval;
    double u3 = pac[ctr+0*corr].p == NULL ? pac[ctr+0*corr].fv : pac[ctr + 0*corr].p->fval;
    double u4 = pac[ctr+1*corr].p == NULL ? pac[ctr+1*corr].fv : pac[ctr + 1*corr].p->fval;
    double u5 = pac[ctr+2*corr].p == NULL ? pac[ctr+2*corr].fv : pac[ctr + 2*corr].p->fval;
    double divf = (f(u4)-f(u3))/(u4-u3);;
    double c = divf == 1.0 ? divf : 1.0;
    D += c*( flux(u2, u3, u4, c) - flux(u1, u2, u3, c) )/Dx;
    //D += ( fluxo_roe(u2,u3,u4,u5, divf, f) - fluxo_roe(u1,u2,u3,u4,divf,f) )/Dx; //(f_{i+0.5} -f_{i-0.5})
  }
  double s3 = pac[ctr].fv;//.p->fval; /* solucao (nesta celula) antes do primeiro passo de rk */
  return s3 - Dt*D;       /* primeiro passo de rk */
}

void avancaTempoMeiov2(int d, int prof, VNo pc[], double  xmin[], double  xmax[], double Dt, Fluxo f, Preditor pred, bool_t op)
{
  int tp = 5;
  int tm = ipow(tp,d);
  int ct = tm/2;
  if(pc[ct].p == NULL){ return; }
  int ex = prof%d;
  int qual;
  for(qual = 0; qual <= 1; qual++)
  { 
    //calcula a posicao da celula filho
    double xminf[d], xmaxf[d];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    
    VNo pcf[tm];
    /* Constroi o pacote do filho {qual}: */
    divide_pacote_geral(d, pc, ex, pred, qual, pcf, op);
    avancaTempoMeiov2(d, prof + 1, pcf, xminf, xmaxf, Dt, f, pred,op);
  }
  assert((pc[ct].p->fil[0] == NULL) == (pc[ct].p->fil[1] == NULL));
  if(pc[ct].p->fil[0] == NULL)
  {
    /* primeiro passo de Runge-Kutta 2 */
//     VFl pcc[2];
//     
//     int dir;
//     int i;
//     for(i=0;i<2;i++){
//       for(dir=0; dir<d;dir++)
//       {
//         pcc[i].p = pc[ct].p;
//         pcc[i].fvs[dir] = pc[ct].p->flxsup[dir];
//         pcc[i].fvi[dir] = pc[ct].p->flxinf[dir];
//       }
//     }
    //pc[ct].p->fv[0] =  rk1ndv2(d, prof, pcc, Dt, ex, xmin, xmax, f);/*  */
    pc[ct].p->fv[0] = rk1_nd( d, prof, pc, Dt, xmin, xmax, f);
  }
}

void avancaTempoUmv2(int d, int prof, VNo pc[], double xmin[], double xmax[], double Dt, Fluxo f, Preditor pred, bool_t op)
{
  int tp = 5;
  int tm = ipow(tp,d);
  int ct = tm/2;
  if(pc[ct].p == NULL){ return; }
  int ex = prof%d;
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    //calcula a posicao da celula filho
    double xminf[d], xmaxf[d];
    limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf);
    VNo pcf[tm];
    /* Constroi o pacote do filho {qual}: */
    divide_pacote_geral(d, pc, ex, pred, qual, pcf,op);
    
    avancaTempoUmv2(d, prof + 1, pcf, xminf, xmaxf, Dt, f, pred,op);
  }
  
  assert((pc[ct].p->fil[0] == NULL) == (pc[ct].p->fil[1] == NULL));
  if(pc[ct].p->fil[0] == NULL)
  {
    /* segundo passo de Runge-Kutta 2 */
//     VFl pcc[2];
//     
//     int dir;
//     int i;
//     for(i=0;i<2;i++){
//     for(dir=0; dir<d;dir++)
//     {
//       pcc[i].p = pc[ct].p;
//       pcc[i].fvs[dir] = pc[ct].p->flxsup[dir];
//       pcc[i].fvi[dir] = pc[ct].p->flxinf[dir];
//       //fprintf(stderr,"dir = %d, flxs = %f, flxi = %f\n",dir,pcc[i].fvs[dir],pcc[i].fvi[dir]);
//     }}
    //pc[ct].p->fv[1] = rk2ndv2( d, prof, pcc, Dt, ex, xmin, xmax, f);/*  */
    pc[ct].p->fv[1] = rk2_nd( d, prof, pc, Dt, xmin, xmax, f);
    //pc[ct].p->fv[1] = rk2nd( d, prof, pc, Dt, ex, xmin, xmax, f);/* adaptativo, mas ainda nao esta funcionando */
    pc[ct].p->fval = pc[ct].p->fv[1];
    //fprintf(stderr,"fv = %f, fval = %f\n", pc[ct].p->fv[1], pc[ct].p->fval);
    assert(pc[ct].p->fv[1]==pc[ct].p->fval);
  }
}

// void avancaTempoMeio(int d, int prof, VNo pac[], double  xmin[], double  xmax[], double Dt, double c, Preditor pred)
// {
//   int ord = 5;
//   int tam = pow(ord,d);
//   int ctr = tam/2;
//   if(pac[ctr].p == NULL){return;}
//   int ex = prof%d;
//   int qual;
//   for (qual = 0; qual < 2; qual++)
//   { 
//     VNo pac_fil[tam];
//     double xmin_fil[d], xmax_fil[d];
//     /* calcula limites da celula filho */
//     limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
//     
//     /* Constroi o pacote do filho {qual}: */
//     divide_pacote_geral(d, ord, pac, ex, pred, qual, pac_fil);
//     
//     avancaTempoMeio(d, prof, pac_fil , xmin_fil, xmax_fil, Dt, c, pred);
//   }
//   assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
//   if(pac[ctr].p->fil[0] == NULL)
//   {
//     /* primeiro passo de Runge-Kutta 2 */
//     pac[ctr].p->fv[0] =  rk1_nd( d, pac, Dt, xmin, xmax, c);
//   }
// }
// 
// void avancaTempoUm(int d, int prof, VNo pac[], double xmin[], double xmax[], double Dt, double c, Preditor pred)
// {
//   int ord = 5;
//   int tam = pow(ord,d);
//   int ctr = tam/2;
//   
//   if(pac[ctr].p == NULL){ return;}
//   int ex = prof%d;
//   int qual;
//   for (qual = 0; qual < d; qual++)
//   { 
//     //calcula a posicao da celula filho
//     double xmin_fil[d], xmax_fil[d];
//     VNo pac_fil[tam];
//     
//     limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
//     
//     /* Constrói o pacote do filho {qual}: */
//     divide_pacote_geral(d, ord, pac, ex, pred, qual, pac_fil);
//     
//     avancaTempoUm(d, prof, pac_fil , xmin_fil, xmax_fil, Dt, c, pred);
//   }
//   
//   assert((pac[ctr].p->fil[0] == NULL) == (pac[ctr].p->fil[1] == NULL));
//   if(pac[ctr].p->fil[0] == NULL)
//   {
//     /* segundo passo de Runge-Kutta 2 */
//     pac[ctr].p->fv[1] = rk2_nd( d, pac, Dt, xmin, xmax, c);
//     pac[ctr].p->fval = pac[ctr].p->fv[1] ;
//   }
// }

double rk2(double u1, double u2, double u3, double u4, double u5, double Dt, double Dx, double c)
{
  double D = c*( flux(u2,u3,u4,c) - flux(u1,u2,u3,c) )/Dx;
  return (u5 + u3 - Dt*D)/2.0;
}


double rk2_nd(int d, int prof, VNo pac[], double Dt, double xmin[], double xmax[], Fluxo f)
{
  double D = 0.0;
  int ex = prof%d;
  int ord = 5;
  int ctr = pow(ord,d)/2;
  int dir; 
  for(dir = 0; dir < d; dir++)
  {
    int corr = ipow(ord,dir);
    double dx = fabs(xmax[dir] - xmin[dir]);
    double Dx = /*ex == dir ? 0.5*(dx + 0.5*dx) :*/ dx;
    double u1 = pac[ctr - 2*corr].p == NULL ? pac[ctr - 2*corr].fv : pac[ctr - 2*corr].p->fval;//fv[0];
    double u2 = pac[ctr - 1*corr].p == NULL ? pac[ctr - 1*corr].fv : pac[ctr - 1*corr].p->fval;//fv[0];
    double u3 = pac[ctr + 0*corr].p == NULL ? pac[ctr + 0*corr].fv : pac[ctr + 0*corr].p->fval;//fv[0];
    double u4 = pac[ctr + 1*corr].p == NULL ? pac[ctr + 1*corr].fv : pac[ctr + 1*corr].p->fval;//fv[0];
    double u5 = pac[ctr + 2*corr].p == NULL ? pac[ctr + 2*corr].fv : pac[ctr + 2*corr].p->fval;//fv[0];
    double divf = (f(u4)-f(u3))/(u4-u3);
    double c = divf==1.0 ? divf : 1.0;
    D += c*( flux(u2,u3,u4,c) - flux(u1,u2,u3,c) )/Dx;
    //D += ( fluxo_roe(u2,u3,u4,u5,divf,f) - fluxo_roe(u1,u2,u3,u4,divf,f) )/Dx; //(f_{i+0.5} -f_{i-0.5})
  }
  double s3 = pac[ctr].p->fval;//[0];  /* solucao (nesta celula) depois de aplicar rk1_nd */
  double s2 = pac[ctr].p->fv[0];//al;  /* solucao (nesta celula) antes de aplicar rk1_nd */
  return  (s2 + s3 - Dt*D)/2.0;  /* segundo passo de rk */
}

double rk2nd(int d, int prof, VFl pc[], double Dt, double ex, double xmin[], double xmax[], Fluxo f)
{
  double D = 0.0;
  int dir; 
  for(dir = 0; dir < d; dir++)
  {
    double dx = fabs(xmax[dir] - xmin[dir]);
    double Dx = /*ex == dir ? 0.5*(dx + 0.5*dx) :*/ dx;
    int k = dir == 0 ? dir : dir + 1;
    if((pc[k].p == NULL) || (pc[k].p->fil[0] == NULL))
    {
      D += (pc[1].p->flxsup[dir] - pc[k].fvi[dir])/Dx; ////??????????????????????
      //D += (pc[1].p->flxsup[dir] - pc[1].p->flxinf[dir])/Dx;
    }
    else
    {
      assert( pc[k].p != NULL );
      double somafluxo = 0.0;
      double xml[d], xMl[d];
      int i;
      for(i=0;i<d;i++)
      {
        xml[i] = i == dir ? xmin[i] - (xmax[i]-xmin[i]) : xmin[i];
        xMl[i] = i == dir ? xmax[i] - (xmax[i]-xmin[i]) : xmax[i];
      }
      encontra_vizinhos_inferiores_esq_refinada( d, pc[k].p, prof, dir, xml, xMl, xmin, &somafluxo);
      D += (pc[1].p->flxsup[dir] - somafluxo)/Dx; ////???????????????????????????
      //D += (pc[1].p->flxsup[dir] - pc[1].p->flxinf[dir])/Dx;
    }
  }
  double s3 = pc[1].p->fval;//[0];  /* solucao (nesta celula) depois de aplicar rk1_nd */
  double s2 = pc[1].p->fv[0];//al;  /* solucao (nesta celula) antes de aplicar rk1_nd */
  return  (s2 + s3 - Dt*D)/2.0;  /* segundo passo de rk */
}