* Second variable
* @return BVN
*/
@Override
public double op(final double x, final double y) {
final TabulatedGaussLegendre gaussLegendreQuad = new TabulatedGaussLegendre(20);
if (Math.abs(correlation) < 0.3) {
gaussLegendreQuad.setOrder(6);
} else if (Math.abs(correlation) < 0.75) {
gaussLegendreQuad.setOrder(12);
}
final double h = -x;
double k = -y;
double hk = h * k;
double bvn = 0.0;
if (Math.abs(correlation) < 0.925) {
if (Math.abs(correlation) > 0) {
final double asr = Math.asin(correlation);
final Eqn3 f = new Eqn3(h, k, asr);
bvn = gaussLegendreQuad.evaluate(f);
bvn *= asr * (0.25 / Math.PI);
}
bvn += cumnorm.op(-h) * cumnorm.op(-k);
} else {
if (correlation < 0) {
k *= -1;
hk *= -1;
}
if (Math.abs(correlation) < 1) {
final double Ass = (1 - correlation) * (1 + correlation);
double a = Math.sqrt(Ass);
final double bs = (h - k) * (h - k);
final double c = (4 - hk) / 8;
final double d = (12 - hk) / 16;
final double asr = -(bs / Ass + hk) / 2;
if (asr > -100) {
bvn = a * Math.exp(asr) * (1 - c * (bs - Ass) * (1 - d * bs / 5) / 3 + c * d * Ass * Ass / 5);
}
if (-hk < 100) {
final double B = Math.sqrt(bs);
bvn -= Math.exp(-hk / 2) * Constants.M_SQRT2PI * cumnorm.op(-B / a) * B
* (1 - c * bs * (1 - d * bs / 5) / 3);
}
a /= 2;
final Eqn6 f = new Eqn6(a, c, d, bs, hk);
bvn += gaussLegendreQuad.evaluate(f);
bvn /= (-2.0 * Math.PI);
}
if (correlation > 0) {
bvn += cumnorm.op(-Math.max(h, k));