class BjerksundStenslandModel {

  protected double[] bivariateNormDiv(double a, double b, boolean posRho) {

    double rho = posRho ? RHO : -RHO;

    double na = NORMAL.getPDF(a);
    double nb = NORMAL.getPDF(b);
    double x1 = (b - rho * a) / RHO_STAR;
    double x2 = (a - rho * b) / RHO_STAR;
    double nx1 = NORMAL.getPDF(x1);
    double nx2 = NORMAL.getPDF(x2);
    double cnx1 = NORMAL.getCDF(x1);
    double cnx2 = NORMAL.getCDF(x2);
    double[] res = new double[5];
    res[0] = na * cnx1; // dB/da
    res[1] = nb * cnx2; // dB/db
    res[2] = -na * (a * cnx1 + rho / RHO_STAR * nx1); // d^2B/da^2
    res[3] = -nb * (b * cnx2 + rho / RHO_STAR * nx2); // d^2B/db^2
    res[4] = na * nx1 / RHO_STAR;
    return res;
  }

}