// Fields of Quantities2 contain the primitive quantities.
// QtoU transforms to conserved quantities as Vector4,
// and UtoQ transforms back.

public immutable class Quantities2 {

  private static final int dimensions = 2;
  private static final int arity = dimensions + 2;
  private static final RectDomain<1> Ddirs = [-1:1:2];

  // MUST SET THESE with gammaSet and setThresholds
  private static double gamma;
  public static double minDensity, minPressure;

  // relative tolerance used in Riemann solver
  private static final double relativeTolerance = 1e-8;

  private final double density, pressure, velX, velY;

  public Quantities2() {
  }

  public final static void gammaSet(double g) {
    gamma = g;
  }

  public final static double gammaGet() {
    return gamma;
  }

  public Quantities2(double density,
		     double velocityX, double velocityY, double pressure) {
    this.density = density;
    this.velX = velocityX;
    this.velY = velocityY;
    this.pressure = pressure;
  }

  public Quantities2(int dim,
		     double density,
		     double velocityN, double velocityT, double pressure) {
    // Initializer with normal and tangential components.
    // Argh.
    this.density = density;
    if (dim == 1) { this.velX = velocityN; this.velY = velocityT; }
    else { this.velY = velocityN; this.velX = velocityT; }
    this.pressure = pressure;
  }

  public final Vector4 Vector4order(int dim, double a, double b, double c, double d) {
    return new Vector4(a,
		       dim, b,
		       3-dim, c,
		       d);
  }

  public final Vector4 CastVector4() {
    return new Vector4(density, velX, velY, pressure);
  }

  public final static Vector4 [2d] CastVector4(Quantities2 [2d] q) {
    RectDomain<2> D = q.domain();
    Vector4 [2d] v = new Vector4[D];
    foreach (p in D) v[p] = q[p].CastVector4();
    return v;
  }

  public final static Quantities2 CastQuantities2(Vector4 v) {
    return new Quantities2(v.elem(1), v.elem(2), v.elem(3), v.elem(4));
  }

  public final static Quantities2 [2d] CastQuantities2(Vector4 [2d] v) {
    RectDomain<2> D = v.domain();
    Quantities2 [2d] q = new Quantities2[D];
    foreach (p in D) q[p] = CastQuantities2(v[p]);
    return q;
  }

  // Retrieval of components

  public final double density() {
    return density;
  }

  public final double velocity(int dim) {
    if (dim == 1) return velX;
    else return velY;
  }

  public final double velocityN(int dim) {
    return velocity(dim);
  }

  public final double velocityT(int dim) {
    return velocityN(dimensions + 1 - dim);
  }

  public final double pressure() {
    return pressure;
  }

  // Derived components

  public final double soundspeed2() {
    return gamma * Math.abs(pressure() / density());
  }

  public final double soundspeed() {
    return Math.sqrt(soundspeed2());
  }

  public final static void setThresholds(double density, double pressure) {
    minDensity = density;
    minPressure = pressure;
  }

  public final Quantities2 thresholded() {
    return new Quantities2(Math.max(density, minDensity), velX, velY,
			   Math.max(pressure, minPressure));
  }

  // Arithmetic operators

  public final Quantities2 op+(Quantities2 v2) {
    return new Quantities2(density + v2.density, velX + v2.velX,
		       velY + v2.velY, pressure + v2.pressure);
  }

  public final Quantities2 op-(Quantities2 v2) {
    return new Quantities2(density - v2.density, velX - v2.velX,
		       velY - v2.velY, pressure - v2.pressure);
  }

  public final Quantities2 op*(Quantities2 v2) {
    return new Quantities2(density * v2.density, velX * v2.velX,
		       velY * v2.velY, pressure * v2.pressure);
  }

  public final Quantities2 op*(double scalar) {
    return new Quantities2(density * scalar, velX * scalar,
		       velY * scalar, pressure * scalar);
  }

  public final Quantities2 op/(double scalar) {
    return new Quantities2(density / scalar, velX / scalar,
		       velY / scalar, pressure / scalar);
  }

  // Other math functions

  public final double norm() {
    return Math.sqrt(density * density + velX * velX + velY * velY + pressure * pressure);
  }

  public final static Quantities2 abs(Quantities2 v) {
    return new Quantities2(Math.abs(v.density), Math.abs(v.velX),
		       Math.abs(v.velY), Math.abs(v.pressure));
  }

  public final static Quantities2 sgn(Quantities2 v) {
    return new Quantities2(sgn(v.density), sgn(v.velX),
		       sgn(v.velY), sgn(v.pressure));
  }

  private final static double sgn(double x) {
    double y;
    if (x > 0.0) {
      y = 1.0;
    } else if (x < 0.0) {
      y = -1.0;
    } else {
      y = 0.0;
    }
    return y;
  }

  public final static Quantities2 min(Quantities2 v1, Quantities2 v2) {
    return new Quantities2(Math.min(v1.density, v2.density), Math.min(v1.velX, v2.velX),
		       Math.min(v1.velY, v2.velY), Math.min(v1.pressure, v2.pressure));
  }

  public final static Quantities2 UtoQ(Vector4 U) {
    double density = Math.max(U.elem(1), minDensity);
    double velocityX = U.elem(2) / density;
    double velocityY = U.elem(3) / density;
    double pressure = Math.max(minPressure, (gamma - 1.0) * (U.elem(4) -
       density * 0.5 * (velocityX * velocityX + velocityY * velocityY)));

    return new Quantities2(density, velocityX, velocityY, pressure);
  }

  public final static Quantities2 [2d] UtoQ(Vector4 [2d] U) {
    RectDomain<2> D = U.domain();
    Quantities2 [2d] q = new Quantities2[D];
    foreach (p in D) q[p] = UtoQ(U[p]);
    return q;
  }

  public final Vector4 QtoU() {
    double dens = Math.max(density(), minDensity);
    return new
      Vector4(dens,
	      dens * velocity(1),
	      dens * velocity(2),
	      pressure() / (gamma - 1.0) + 
	      dens * 0.5 *
	      (velocity(1)*velocity(1) + velocity(2)*velocity(2)));
  }

  public final static Vector4 [2d] QtoU(Quantities2 [2d] q) {
    RectDomain<2> D = q.domain();
    Vector4 [2d] U = new Vector4[D];
    foreach (p in D) U[p] = q[p].QtoU();
    return U;
  }

  public final Quantities2 qeig0(int dim, double c) {
    // eigenvectors, used with deltaq
    return new
      Quantities2(dim,
                  density() - pressure() / (c * c),
                  0.0,
                  velocityT(dim),
                  0.0);
  }

  public final Quantities2 qeig(int sdim, Quantities2 q, double c) {
    // eigenvectors, used with deltaq
    int dim = (sdim > 0) ? sdim : -sdim;
    double sdensity = (sdim > 0) ? q.density() : -q.density();
    return new
      Quantities2(dim,
                  sdensity / c * velocityN(dim) +
                  pressure() / (c * c),
                  velocityN(dim) + pressure() / (sdensity * c),
                  0.0,
                  sdensity * c * velocityN(dim) + pressure()) * 0.5;
  }

  public final Vector4 getFlux(int dim) {
    double vN = velocity(dim);
    double vT = velocity(3 - dim);
    double bigterm = (pressure() / (1.0 - 1.0 / gamma) +
		      density() * 0.5 * (vN*vN + vT*vT));
    return Vector4order(dim, density() * vN,
			density() * vN * vN + pressure(),
			density() * vN * vT,
			vN * bigterm);
  }

  public final Vector4 solidBndryFlux(int signdim) {
    int dim = (signdim > 0) ? signdim : -signdim;
    double dir = (signdim > 0) ? 1.0 : -1.0;
    double pstar = pressure() +
      dir * density() * soundspeed() * velocityN(dim);
    return Vector4order(dim, 0.0, pstar, 0.0, 0.0);
  }


  private static double [1d] eigval = new double[-1:1];

  public final static void
    extrapolate(Quantities2 [2d] q, Quantities2 [2d] deltaq,
		int dim, double dtdx,
		Quantities2 [1d][2d] qhext) {

    RectDomain<2> D = q.domain();
    double halfdtdx = 0.5 * dtdx;

    foreach (p in D) {
      double soundspeed = q[p].soundspeed();

      foreach (ps in [-1:1])
	eigval[ps] = q[p].velocityN(dim) + (double)(ps[1]) * soundspeed;

      foreach (pd in Ddirs) {

	int ipd = pd[1]; Point<1> pdc = [-ipd]; double dpd = (double) ipd;

	qhext[pd][p] = q[p] + deltaq[p] *
	  ((0.5 - halfdtdx * Math.max(dpd * eigval[pd], 0.0)) * dpd);

	if (dpd * eigval[0] > 0.0) {
	  qhext[pd][p] = qhext[pd][p] +
	    deltaq[p].qeig0(dim, soundspeed) *
	    (halfdtdx * (eigval[pd] - eigval[0]));

	  if (dpd * eigval[pdc] > 0.0) {
	    qhext[pd][p] = qhext[pd][p] +
	      deltaq[p].qeig(-ipd * dim, q[p], soundspeed) *
	      (halfdtdx * (eigval[pd] - eigval[pdc]));
	  }
	}
	qhext[pd][p] = qhext[pd][p].thresholded();
	if (qhext[pd][p].norm() > 1e10) {
	  System.out.println("huge qhext[" + pd[1] + "][" + Format.toString(p) + "]:  " +
			     qhext[pd][p].density() + "  " +
			     qhext[pd][p].velocity(1) + "  " +
			     qhext[pd][p].velocity(2) + "  " +
			     qhext[pd][p].density());
	  System.out.println("huge qhext from q[" + Format.toString(p) + "]:  " +
			     q[p].density() + "  " +
			     q[p].velocity(1) + "  " +
			     q[p].velocity(2) + "  " +
			     q[p].density());
	  System.out.println("huge qhext soundspeed:  " + soundspeed);
	}

      }
    }
  }

  // for slopes

  private final static Quantities2 getdlim(Quantities2 v1, Quantities2 v2) {
    return new Quantities2(getdlim(v1.density, v2.density), getdlim(v1.velX, v2.velX),
			   getdlim(v1.velY, v2.velY), getdlim(v1.pressure, v2.pressure));
  }

  private final static double getdlim(double a, double b) {
    if (a * b > 0) {
      return Math.min(Math.abs(a), Math.abs(b)) * 2.0;
    } else {
      return 0.0;
    }
  }









  public final static Quantities2 vanLeer(Quantities2 a, Quantities2 b, Quantities2 c) {
    // Calculate van Leer slope
    Quantities2 dlim = getdlim(c - b, b - a);
    Quantities2 dif2 = c - a;
    Quantities2 retval = min(abs(dif2) * 0.5, dlim) * sgn(dif2);
    return retval;
  }
    
  public final static void Slopes2(Quantities2 [2d] v, int dim,
				   Quantities2 [2d] df) {
    Slopes2(v, dim, true, true, df);
  }

  public final static void
    Slopes2(Quantities2 [2d] v, int dim, boolean lowExt, boolean highExt,
	    Quantities2 [2d] df) {
    /*
      2nd-order slopes, with limiting in interior
      lowExt:  true iff low end is at exterior boundary
      highExt:  true iff high end is at exterior boundary
      
      If an end is at an exterior boundary, then you need to
      do one-sided differencing there.
    */

    // foreach (p in df.domain()) df[p] = new Quantities2(0.0, 0.0, 0.0, 0.0);

    //  dim is 1 or 2, giving unit [1, 0] or [0, 1]
    Point<2> unit = Point<2>.direction(dim);

    RectDomain<2> D = v.domain();
    RectDomain<2> Dinner1 = D.shrink(1, dim).shrink(1, -dim);

    RectDomain<2> Dget = Dinner1;
    if (lowExt) Dget = Dget.accrete(1, -dim);
    if (highExt) Dget = Dget.accrete(1, +dim);

    // Quantities2 [2d] df = new Quantities2[Dget];

    foreach (p in Dinner1) df[p] = vanLeer(v[p - unit], v[p], v[p + unit]);

    // one-sided differencing at ends where exterior boundaries lie;
    // here the min/max are same for D, Dget
    if (lowExt) foreach (p in Dget.border(1, -dim, 0))
      df[p] = v[p] * (-1.5) + v[p + unit] * 2.0 - v[p + unit + unit] * 0.5;

    if (highExt) foreach (p in Dget.border(1, dim, 0))
      df[p] = v[p] * 1.5 - v[p - unit] * 2.0 + v[p - unit - unit] * 0.5;
  }


  private static double [1d] sidespeed = new double[Ddirs];
  private static double [1d] w = new double[Ddirs];
  private static Quantities2 [1d] qsides = new Quantities2[Ddirs];

  public final static Quantities2
    Riemann0(int dim, Quantities2 qleft, Quantities2 qright) {
    // Solve a Riemann problem
    
    Quantities2 qdiff = qleft - qright;
    if (qdiff.norm() <= relativeTolerance * qleft.norm()) return qleft;

    // Assumes density and pressure of input quantities are positive.

    qsides[-1] = qleft; qsides[+1] = qright;

    foreach (pd in Ddirs) {
      sidespeed[pd] = qsides[pd].soundspeed();
      w[pd] = sidespeed[pd] * qsides[pd].density();
    }

    double p = (qsides[-1].pressure() * w[+1] +
		qsides[+1].pressure() * w[-1] +
		w[-1] * w[+1] *
		(qsides[-1].velocityN(dim) - qsides[+1].velocityN(dim))) /
      (w[-1] + w[+1]);
    p = Math.max(p, minPressure);

    double u = (qsides[+1].velocityN(dim) * w[+1] +
                qsides[-1].velocityN(dim) * w[-1] +
                qsides[-1].pressure() - qsides[+1].pressure()) /
      (w[-1] + w[+1]);

    int side = (u >= 0.0) ? -1 : +1;
    double usgn = (float) (-side);

    double density = qsides[side].density() +
      (p - qsides[side].pressure()) /
      (sidespeed[side] * sidespeed[side]);
    density = Math.max(density, minDensity);
  
    Quantities2 qstar = new
      Quantities2(dim, density, u, qsides[side].velocityT(dim), p);

    double ustar = u - usgn * qstar.soundspeed();
    double v = qsides[side].velocityN(dim) - usgn * sidespeed[side];
    if (p > qsides[side].pressure()) {
      // (A) or (C), shock at this side
      double shockspeed = (ustar + v) * 0.5;
      if (shockspeed * usgn >= 0.0) {
        return qsides[side];
      } else {
        return qstar;
      }
    } else {
      // (B) or (D), rarefaction at this side
      if (v * usgn >= 0.0) {
        return qsides[side];
      } else if (ustar * usgn <= 0.0) {
        return qstar;
      } else {
        // -usgn * v > 0, usgn * ustar > 0
        // hence -usgn * (v - ustar) > 0
        double vustardif = v - ustar;
        double frac;
        if (-usgn * vustardif < relativeTolerance * (-usgn * v)) {
          frac = 1.0;
        } else {
          frac = v / vustardif;
          frac = Math.max(0.0, Math.min(frac, 1.0));      
        }
        return qstar * frac + qsides[side] * (1.0 - frac);
      }
    }
  }

    
}