public immutable class Vector4 {

  private static int arity = 4;

  protected double val0, val1, val2, val3;

  public Vector4() {
  }

  public Vector4(double a, double b, double c, double d) {
    val0 = a; val1 = b; val2 = c; val3 = d;
  }

  public Vector4(double a,
		 int i2, double a2,
		 int i3, double a3,
		 double d) {
    /* Argh. */
    val0 = a;
    if (i2 == 1) { val1 = a2; val2 = a3; }
    else { val1 = a3; val2 = a2; }
    val3 = d;
  }

  public final double elem(int j) {
    // elements are indexed  1..arity
    if (j == 1) return val0;
    if (j == 2) return val1;
    if (j == 3) return val2;
    if (j == 4) return val3;
    return 0.0;
  }

  public final Vector4 op+(Vector4 v2) {
    return new Vector4(val0 + v2.val0, val1 + v2.val1,
		       val2 + v2.val2, val3 + v2.val3);
  }

  public final Vector4 op-(Vector4 v2) {
    return new Vector4(val0 - v2.val0, val1 - v2.val1,
		       val2 - v2.val2, val3 - v2.val3);
  }

  public final Vector4 op*(Vector4 v2) {
    return new Vector4(val0 * v2.val0, val1 * v2.val1,
		       val2 * v2.val2, val3 * v2.val3);
  }

  public final Vector4 op*(double scalar) {
    return new Vector4(val0 * scalar, val1 * scalar,
		       val2 * scalar, val3 * scalar);
  }

  public final Vector4 op/(double scalar) {
    return new Vector4(val0 / scalar, val1 / scalar,
		       val2 / scalar, val3 / scalar);
  }

  // Other math functions

  public final double norm() {
    return Math.sqrt(val0 * val0 + val1 * val1 + val2 * val2 + val3 * val3);
  }

  public final static Vector4 abs(Vector4 v) {
    return new Vector4(Math.abs(v.val0), Math.abs(v.val1),
		       Math.abs(v.val2), Math.abs(v.val3));
  }

  public final static Vector4 sgn(Vector4 v) {
    return new Vector4(sgn(v.val0), sgn(v.val1),
		       sgn(v.val2), sgn(v.val3));
  }

  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 Vector4 min(Vector4 v1, Vector4 v2) {
    return new Vector4(Math.min(v1.val0, v2.val0), Math.min(v1.val1, v2.val1),
		       Math.min(v1.val2, v2.val2), Math.min(v1.val3, v2.val3));
  }

  // average over a RectDomain<2>

  public final static Vector4 Average(Vector4 [2d] v) {
    Vector4 sum = new Vector4(0.0, 0.0, 0.0, 0.0);
    RectDomain<2> D = v.domain();
    foreach (p in D) sum = sum + v[p];
    return (sum / (double) D.size());
  }
  
  // routines for slopes

  private final static Vector4 getdlim(Vector4 v1, Vector4 v2) {
    return new Vector4(getdlim(v1.val0, v2.val0), getdlim(v1.val1, v2.val1),
		       getdlim(v1.val2, v2.val2), getdlim(v1.val3, v2.val3));
  }

  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 Vector4 vanLeer(Vector4 a, Vector4 b, Vector4 c) {
    // Calculate van Leer slope
    Vector4 dlim = getdlim(c - b, b - a);
    Vector4 dif2 = c - a;
    Vector4 retval = min(abs(dif2) * 0.5, dlim) * sgn(dif2);
    return retval;
  }
    
  public final static Vector4 [2d] Slopes2(Vector4 [2d] v, int dim) {
    return Slopes2(v, dim, true, true);
  }

  public final static Vector4 [2d]
    Slopes2(Vector4 [2d] v, int dim, boolean lowExt, boolean highExt) {
    /*
      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.
    */

    //  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);

    Vector4 [2d] df = new Vector4[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;

    return df;
  }

}