/***************************************************************************
 *   Copyright (C) 2009 by Douglas Castro   *
 *   douglas@ime.unicamp.br   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "timestep.h"


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 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;
}

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

    double fa[] = {pac[ctr-2*corr].fv/*p->fval*/, pac[ctr-corr].fv/*p->fval*/, pac[ctr].fv/*p->fval*/, pac[ctr+corr].fv/*p->fval*/, pac[ctr+2*corr].fv/*p->fval*/};
    
    pac[ctr].p->fv[0] =  rk1(fa[0],fa[1],fa[2],fa[3],fa[4], Dt, Dx, c);

    if((pac[ctr].p->fil[0]!=NULL)&&(pac[ctr].p->fil[0]->virtual_leaf))
    {
      int ex = prof % d;
      int corr = pow(ord,ex);
      double x0 = pac[ctr-corr].fv;
      double x1 = pac[ctr].fv;
      double x2 = pac[ctr+corr].fv;
      double pe = pred(x0, x1, x2);
      double pd = 2.0*x1 - pe;
      pac[ctr].p->fil[0]->fv[0] = pe /*pac[ctr].p->fv[0]*/;
      pac[ctr].p->fil[1]->fv[0] = pd /*pac[ctr].p->fv[0]*/;
    }
  }
 
}

void avancaTempoUm(int d, int prof, VNo pac[], double xmin[], double xmax[], double Dt, double c, Preditor pred)
{
  int ord = 5;
  int ctr = pow(ord,d)/2;
  if(pac[ctr].p==NULL){return;}
  int ex = prof%d;
  int qual;
  for (qual = 0; qual <= 1; qual++)
  { 
    //calcula a posicao da celula filho
    int ex = prof%d;
    double xmin_fil[d], xmax_fil[d];
    
    limites_celula( d, ex, qual, xmin, xmax, xmin_fil, xmax_fil);
    
    int tam = pow(ord,d);
    VNo pac_fil[tam];
    
    /* Constrói o pacote do filho {qual}: */
    divide_pacote_geral(d, pac, ex, pred, qual, pac_fil);
    
    avancaTempoUm(d, prof, pac_fil , xmin_fil, xmax_fil, Dt,c, pred);
  }
  if(pac[ctr].p->leaf)
  {
    /* segundo passo de Runge-Kutta 2 */
    int corr = pow(ord,ex);
    double Dx = fabs(xmax[ex] - xmin[ex]);
    int i;
    for(i=-2;i<3;i++){
      if(pac[ctr+i*corr].p==NULL)
        { printf("erro - avancaTempoUm - pac[ctr %d].p == NULL\n",i);}
    }
    double t0 = pac[ctr-2*corr].p->fv[0];
    double t1 = pac[ctr-corr].p->fv[0];
    double t2 = pac[ctr].p->fv[0];
    double t3 = pac[ctr+corr].p->fv[0];
    double t4 = pac[ctr+2*corr].p->fv[0];
    double fa[] = {t0, t1, t2, t3, t4};
    double temp = pac[ctr].p->fval;
    pac[ctr].p->fval =  rk2(fa[0], fa[1], fa[2], fa[3], temp, Dt, Dx, c);
//     return;
  }
  
}