Двойной интеграл. Монте-Карло

[Astaroth~]

Дан кратный интеграл, см вложения.


Интеграл беру не до бесконечности, а ограничиваю XENDVALUE, YENDVALUE

q0,z0, L, r - константы

Нужна помощь, пожалуйста проверьте мою версию реализации, дайте советы.
Заранее спасибо.

public static double f2(double u, double z, double r, double z0)
{ //f2(du, dz, r, z0)
     double p, x;
      double sechsqr, zz0;
      x = (r * r + u * u + z * z) / (2.0 * r * u);
      p = x - Math.Sqrt(x * x - 1.0);
 
   
       return u * (-(q0* Math.Exp(-u / L) * Math.Pow(Math.Cosh(z / z0), -2.0)) / L) *
             2.0 * (alglib.ellipticintegralk(p) - alglib.ellipticintegrale(p)) / (Math.PI * Math.Sqrt(r * u * p));
 
 }
 
 // x = du, y = dz, z = f
public static double MMC2D(double N, double r, double z0, double XBEGINVALUE, double
   YBEGINVALUE, double ZBEGINVALUE, double XENDVALUE, double YENDVALUE, double ZENDVALUE)
{
        double du, dz, f;
        double result;
 
        double ratio = 0.0;
        Random rndt = new Random(1243);
        Random rnd = new Random(rndt.Next());
 
        for (double i = 0; i < N; i++)
        {
             rnd.NextDouble();
 
             du = XBEGINVALUE + XENDVALUE * rnd.NextDouble();
             dz = YBEGINVALUE + YENDVALUE * rnd.NextDouble();
             f = ZBEGINVALUE + ZENDVALUE * rnd.NextDouble();
 
             r += (f <= f2(du, dz, r, z0)) ? 1.0 : 0.0;
          }
                      
          ratio = r / N;
          result = ratio * (XENDVALUE - XBEGINVALUE) * (YENDVALUE - YBEGINVALUE) *(ZENDVALUE - ZBEGINVALUE);
          return result;
}

Вызов

double result= MMC2D(100000.0, R , h, 0.0, 0.0, 0.0, A, B, C);
 
double result2 = 4.0 * Math.PI * result;