public class Ramp2dmInit extends Ramp2dm implements InitialPatch2 {

  public final double [1d] dsGet(Point<2> extent) {
    
    double [1d] ds = new double[Ddims];
    
    foreach (pdim in Ddims) ds[pdim] =
      (ends[pdim][+1] - ends[pdim][-1]) / (double) extent[pdim[1]];

    return ds;
  }

  private int nsub = 10;  // number of subintervals
  private double dnsubinv = 1.0 / (double) nsub;
  private double rn2 = dnsubinv * dnsubinv;

  public final Quantities2 [2d]
    qInitial(RectDomain<2> Dlocal, double [1d] ds) {

    Quantities2 [2d] q = new Quantities2[Dlocal];
    double [1d] subds = new double[Ddims];
    foreach (pdim in Ddims) subds[pdim] = ds[pdim] * dnsubinv;
    foreach (p in Dlocal) {
      double left = (double) p[1] * ds[1];
      double right = (double) (p[1] + 1) * ds[1];
      double bottom = (double) p[2] * ds[2];
      double top = (double) (p[2] + 1) * ds[2];
      if (right <= shockloc + bottom / slope) {
	q[p] = q1;
      } else if (left >= shockloc + top / slope) {
	q[p] = q0;
      } else {
	int iwts = 0;
	for (int jsub = 1; jsub <= nsub; jsub++) {
	  for (int ksub = 1; ksub <= nsub; ksub++) {
	    double xsub = left + ((double) jsub - 0.5) * subds[1];
	    double ysub = bottom + ((double) ksub - 0.5) * subds[2];
	    if (xsub <= shockloc + ysub / slope) iwts++;
	  }
	}
	double wt = (double) iwts * rn2;
	q[p] = q1 * wt + q0 * (1.0 - wt);
      }
    }
    return q;
  }
}