double f2(double u, double z, double r, double z0, double h)
{
    double Kp, Ep; // Elliptic
    double p, x;
    double sechsqr, zz0;
    x = (r * r + u * u + z * z) / (2.0 * r * u);
    p = x - sqrt(x*x - 1.0);
    zz0 = z / z0;
    sechsqr = 2.0 / (Exp(zz0) + Exp(-zz0))
    sechsqr *= sechsqr;

    // Kp, Ep = ? (alglib)
    // ...
    //

    return (u * sechsqr * (Kp - Ep) / sqrt(u * p));
}

// x = du, y = dz, z = f
double MMC2D(double N, double r, double z0, double h, double XBEGINVALUE, double YBEGINVALUE, double ZBEGINVALUE, double XENDVALUE, double YENDVALUE, double ZENDVALUE)
{

        double du, dz, f;
        double result;

        double r = 0;
        double ratio = 0.0;
        Random rndt= new Random(1243);
        rndt.Next();
        Random rnd = new Random(rndt);

        for(double i = 0; i < N; i++)
        {

                du = XBEGINVALUE + Math.Tan(XENDVALUE*rnd.NextDouble());
                dz = YBEGINVALUE + Math.Tan(YENDVALUE*rnd.NextDouble());
                f  = ZBEGINVALUE + ZENDVALUE*rnd.NextDouble();

                r += (f <= f2(du, dz, r, z0, h)) ? 1.0 : 0.0;

                //uncomment below to enter the matrix!
                //ratio = (double)r/(double)i;
                //std::cout << "Iteration: " << i << "\tArea: " << ratio*2.0 << "\tRatio: " << ratio << "\tx: " << x << "\ty: " << y << std::endl;
        }
        ratio = r/N;
//      cout.precision(12);
//        std::cout << "Ratio: " << ratio << std::endl;
//        std::cout << "\nResult: " << (DIMENSIONS-1)*ratio*(XENDVALUE-XBEGINVALUE)*(YENDVALUE-YBEGINVALUE) << std::endl;
        result = ratio * (XENDVALUE - XBEGINVALUE) * (YENDVALUE - YBEGINVALUE) * (ZENDVALUE - ZBEGINVALUE);

        return result;
}