class G2ppPiecewiseConstantModel {

  public double[][] gamma(final G2ppPiecewiseConstantParameters g2parameters, double theta0, double theta1) {
    double[] a = g2parameters.getMeanReversion();
    DoubleArrayList[] sigma = g2parameters.getVolatility();
    int indexStart = 1; // Period in which the time startExpiry is; _volatilityTime[i-1] <= startExpiry < _volatilityTime[i];
    while (theta0 > g2parameters.getVolatilityTime()[indexStart]) {
      indexStart++;
    }
    int indexEnd = indexStart; // Period in which the time endExpiry is; _volatilityTime[i-1] <= endExpiry < _volatilityTime[i];
    while (theta1 > g2parameters.getVolatilityTime()[indexEnd]) {
      indexEnd++;
    }
    int sLen = indexEnd - indexStart + 2;
    double[] s = new double[sLen];
    s[0] = theta0;
    System.arraycopy(g2parameters.getVolatilityTime(), indexStart, s, 1, sLen - 2);
    s[sLen - 1] = theta1;

    double[] gammaii = new double[2];
    double gamma12 = 0.0;
    double[][] exp2as = new double[2][sLen];
    double[] expa0a1s = new double[sLen];
    for (int loopindex = 0; loopindex < sLen; loopindex++) {
      for (int loop = 0; loop < 2; loop++) {
        exp2as[loop][loopindex] = Math.exp(2 * a[loop] * s[loopindex]);
      }
      expa0a1s[loopindex] = Math.exp((a[0] + a[1]) * s[loopindex]);
    }
    for (int loopindex = 0; loopindex < sLen - 1; loopindex++) {
      for (int loop = 0; loop < 2; loop++) {
        gammaii[loop] += sigma[loop].get(indexStart - 1 + loopindex) * sigma[loop].get(indexStart - 1 + loopindex) * (exp2as[loop][loopindex + 1] - exp2as[loop][loopindex]);
      }
      gamma12 += sigma[0].get(indexStart - 1 + loopindex) * sigma[1].get(indexStart - 1 + loopindex) * (expa0a1s[loopindex + 1] - expa0a1s[loopindex]);
    }
    double[][] result = new double[2][2];
    result[0][0] = gammaii[0] / (2 * a[0]);
    result[1][1] = gammaii[1] / (2 * a[1]);
    result[1][0] = gamma12 / (a[0] + a[1]);
    result[0][1] = result[1][0];
    return result;
  }

}