text: #include"pluto.h"
  
  /* *********************************************************************** */
  double MP5_Reconstruct(double *F, double dx, int j)
  /*  
   *
   * 
   *  
   *  
   *
   *
   ************************************************************************* */
  {
    double f, d2, d2p, d2m; 
    double dMMm, dMMp;
    double scrh1,scrh2, fmin, fmax; 
    double fAV, fMD, fLC, fUL, fMP;
    static double alpha = 4.0, epsm = 1.e-12;
  
    f  = 2.0*F[j-2] - 13.0*F[j-1] + 47.0*F[j] + 27.0*F[j+1] - 3.0*F[j+2];
    f /= 60.0;   
  
    fMP = F[j] + MINMOD(F[j+1]-F[j],alpha*(F[j]-F[j-1]));
  
    if ((f - F[j])*(f - fMP) <= epsm) return f;
  
    d2m = F[j-2] + F[j  ] - 2.0*F[j-1];    /* -- Eq. (2.19) -- */
    d2  = F[j-1] + F[j+1] - 2.0*F[j];
    d2p = F[j  ] + F[j+2] - 2.0*F[j+1];    /* -- Eq. (2.19) -- */
  
    scrh1 = MINMOD(4.0*d2 - d2p, 4.0*d2p - d2);
    scrh2 = MINMOD(d2, d2p);
    dMMp  = MINMOD(scrh1,scrh2);   /* -- Eq. (2.27) -- */
  
    scrh1 = MINMOD(4.0*d2m - d2, 4.0*d2 - d2m);
    scrh2 = MINMOD(d2, d2m);
    dMMm  = MINMOD(scrh1,scrh2);   /* -- Eq. (2.27) -- */
  
    fUL = F[j] + alpha*(F[j] - F[j-1]);   /* -- Eq. (2.8) -- */
    fAV = 0.5*(F[j] + F[j+1]);        /* -- Eq. (2.16) -- */
    fMD = fAV - 0.5*dMMp; /* -- Eq. (2.28) -- */
    fLC = 0.5*(3.0*F[j] - F[j-1]) + 4.0/3.0*dMMm;  /* -- Eq. (2.29) -- */
  
    scrh1 = MIN(F[j], F[j+1]); scrh1 = MIN(scrh1, fMD);
    scrh2 = MIN(F[j], fUL);    scrh2 = MIN(scrh2, fLC);
    fmin  = MAX(scrh1, scrh2);  /* -- Eq. (2.24a) -- */
  
    scrh1 = MAX(F[j], F[j+1]); scrh1 = MAX(scrh1, fMD);
    scrh2 = MAX(F[j], fUL);    scrh2 = MAX(scrh2, fLC);
    fmax  = MIN(scrh1, scrh2);  /* -- Eq. (2.24b) -- */
  
    f = Median(f, fmin, fmax); /* -- Eq. (2.26) -- */
    return f;
  }
  
  /* ************************************************************* */
  double PPM_Reconstruct (double *v, double dx, int i)
  /*
   *
   *
   *
   *************************************************************** */
  {
    double dvm2, dvm1, dvp1, dvp2;
    double dvc, Sm1, Sp1, Sp2, SM;
    double ap, am, vp, vm;
  
    dvm2 = v[i-1] - v[i-2];
    dvm1 = v[i]   - v[i-1];
    dvp1 = v[i+1] - v[i];
    dvp2 = v[i+2] - v[i+1];
  
    dvc = 0.5*(dvm2 + dvm1);
    SM  = 2.0*MINMOD(dvm2, dvm1);
    Sm1 = MINMOD(dvc, SM);
  
    dvc = 0.5*(dvm1 + dvp1);
    SM  = 2.0*MINMOD(dvm1, dvp1);
    Sp1 = MINMOD(dvc, SM);
  
    dvc = 0.5*(dvp2 + dvp1);
    SM  = 2.0*MINMOD(dvp2, dvp1);
    Sp2 = MINMOD(dvc, SM);
  
    vp = 0.5*(v[i] + v[i+1]) - (Sp2 - Sp1)/6.0;
    vm = 0.5*(v[i] + v[i-1]) - (Sp1 - Sm1)/6.0;
  
    ap = vp - v[i];
    am = vm - v[i];
  
    if (ap*am >= 0.0) {
      ap = am = 0.0;
    }else{
      if (fabs(ap) >= 2.0*fabs(am)) ap = -2.0*am;
      if (fabs(am) >= 2.0*fabs(ap)) am = -2.0*ap;
    }
    vp = v[i] + ap;
    return vp;
  /*  vm = v[i] + am; */
  }
  
  /* ************************************************************* */
  double LIMO3_Reconstruct(double *v, double dx,  int i)
  /*
   *
   *
   *
   *************************************************************** */
  {
    double dvp, dvm, r = 1.;
    double a,b,c,q, th, lim;
    double eta, psi, eps = 1.e-12;
  
    dvp = v[i+1] - v[i];
    dvm = v[i] - v[i-1];
  
    th  = dvm/(dvp + 1.e-16);
  
    q = (2.0 + th)/3.0;
    a = MIN(1.5,2.0*th);
    a = MIN(q,a);
    b = MAX(-0.5*th,a);
    c = MIN(q,b);
    psi = MAX(0.0,c);
  
    eta = r*dx;
    eta = (dvm*dvm + dvp*dvp)/(eta*eta);
    if ( eta <= 1.0 - eps) {
      lim = q;
    }else if (eta >= 1.0 + eps){
      lim = psi;
    }else{
      psi =   (1.0 - (eta - 1.0)/eps)*q
            + (1.0 + (eta - 1.0)/eps)*psi;
      lim = 0.5*psi;
    }
    return (v[i] + 0.5*dvp*lim);
  }
  /* *********************************************************************** */
  double WENOZ_Reconstruct(double *F, double dx,  int j)
  /*  
   *
   * 
   *  Provide interface values F_{i+1/2} using 
   *  WENO-Z reconstruction.
   *   
   *
   *
   ************************************************************************* */
  #define eps 1.e-6
  {
    double IS0, IS1, IS2, IS3, IS4, IS5;
    double a0, a1, a2;
    double f0, f1, f2;
    double w0, w1, w2;
    double sum_a, f;
    static double thirteen_12 = 13.0/12.0;
    double b0, b1, b2,t5;
  
    a0 = F[j-2] - 2.0*F[j-1] +     F[j]; 
    a1 = F[j-2] - 4.0*F[j-1] + 3.0*F[j]; 
    b0 = thirteen_12*a0*a0 + 0.25*a1*a1;
  
    a0 = F[j-1] - 2.0*F[j] + F[j+1]; 
    a1 = F[j-1] - F[j+1]; 
    b1 = thirteen_12*a0*a0 + 0.25*a1*a1;
  
    a0 =     F[j] - 2.0*F[j+1] + F[j+2]; 
    a1 = 3.0*F[j] - 4.0*F[j+1] + F[j+2]; 
    b2 = thirteen_12*a0*a0 + 0.25*a1*a1;
  
    t5 = fabs(b0-b2);
  
    a0 = 1.0*(1.0 + t5/(b0 + 1.e-40));
    a1 = 6.0*(1.0 + t5/(b1 + 1.e-40));
    a2 = 3.0*(1.0 + t5/(b2 + 1.e-40));
  
    sum_a = 1.0/(a0 + a1 + a2);
    w0 = a0*sum_a;
    w1 = a1*sum_a;
    w2 = a2*sum_a;
  
    f0 = (2.0*F[j-2] - 7.0*F[j-1] + 11.0*F[j]);
    f1 = (   -F[j-1] + 5.0*F[j]   +  2.0*F[j+1]);
    f2 = (2.0*F[j]   + 5.0*F[j+1] -      F[j+2]);
  
    f = (w0*f0 + w1*f1 + w2*f2)/6.0;
    return(f);
  }
  
  /* *********************************************************************** */
  double WENO3_Reconstruct(double *F, double dx, int j)
  /*  
   *
   * 
   *  Provide interface values F_{i+1/2} using 
   *  third-order WENO reconstruction.
   *   
   *
   *
   ************************************************************************* */
  {
    double a0, b0, w0, f0, a1, b1, w1, f1;
    double tau, dx2;
  
    dx2 = dx*dx;
  
    b0 = F[j+1] - F[j];
    b1 = F[j-1] - F[j];
  
    b0 = b0*b0;
    b1 = b1*b1;
  
  /* -- traditional version -- */
  
  /*
    a0 = eps + b0;
    a1 = eps + b1;
    
    a0 = 2.0/(a0*a0);
    a1 = 1.0/(a1*a1);
  */
  /* -- improved version -- */
  
    tau = (F[j+1] - 2.0*F[j] + F[j-1]);
    tau = tau*tau;
  
    a0 = 2.0*(1.0 + tau/(dx2 + b0));
    a1 = 1.0*(1.0 + tau/(dx2 + b1));
  
    w0 = a0/(a0 + a1);
    w1 = a1/(a0 + a1);
  
    f0 =  0.5*( F[j+1] +     F[j]);
    f1 =  0.5*(-F[j-1] + 3.0*F[j]);
    return (w0*f0 + w1*f1);
  }
  #undef eps
  
  
  /* *********************************************************************** */
  double LIN_Reconstruct(double *F, double dx, int j)
  /*  
   *
   * 
   *  Provide interface values F_{i+1/2} using 
   *  simple slope-limited linear reconstruction.
   *   
   *
   *
   ************************************************************************* */
  {
    double dFp, dFm, dF;
  
    dFp = F[j+1] - F[j];
    dFm = F[j] - F[j-1];
  
    dF = MINMOD(dFp, dFm);
    return (F[j] + 0.5*dF);
  }
  
  
  /* ************************************************************* */
  double Median (double a, double b, double c)
  /*
   *
   *
   *
   *************************************************************** */
  {
    return (a + MINMOD(b-a,c-a));
  }
  