Численное интегрирование(проверить код) - C#

Формулировка задачи:

Здравствуйте! Помогите пожалуйста разобраться. Bm всегда выводит значение NaN. Где ошибка? И правильно ли вообще построен код?

                    Листинг программы
                
 
                            
                        
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
 
namespace ConsoleApplication1
{
    class Program
    {
        static double V0,alfak,fk,lambda,gm,result,a0,b0,a,ak,bk,en,b,rect,e;
        
        static int i,j;
        static int nn = 100;
        static int ksi = 1;
        static double m;
 
        public static double Psi(double m, double ksi)
        {
           
                double a = 0.0063f;
                double s = 0.2f;
                double Gamast = 2.5e6f;
                double Tau = 0.09f;
                if (ksi <= a * Math.PI / Tau)
                    return 1 / m - (1 - s) * V0;
                else return Gamast / (m * Gamast);
 
        }
 
        public static double P(double x)
        {
            double h = 0.0104f;
            return (2*h*Math.Abs(x))/(Math.Exp(2*Math.Abs(x))-1)*(1-Math.Exp(- Math.Abs(m)));
        }
 
        //подынтегральная функция
        public static double fx(double x)
        {
            int y = 0;
            double delta = 0.0116f;
            double l = 0.2502f;
            double Mu0 = 12.56e-7f;
            double B0=(2 * Mu0 * Math.Sign(m) * Math.Cosh(m * y)) / ((m - alfak) * Math.Sinh(Math.Abs(m) * delta)) * ((Math.Cos(m * l) - 1) + j * Math.Sin(m * l));
            return gm * B0 * fk;
        }
 
        //решение методом прямоугольников
        private void rectangls()
        {
            double s1;
            double s2;
            double h;
            
            int i,n;
 
            n = 100;
            do
            {
                s1=0;
            s2=0;
            h=(b-a)/n;
 
            for (i=1;i<n;i++)
            {
                s1=s1+h*fx(h*i-h/2);
            }
            n=n*2;
            h=(b-a)/n;
            for (i=1;i<n;i++)
            {
                s2=s2+h*fx(h*i-h/2);
            }
            }
            while (Math.Abs(s2-s1)<e);
            rect=s2;
 
        }
 
        //решение трансцендентного уравнения методом ньютона
        private void ras4et()
        {
            double eps = 0.0001;
            for (m = -79.5; m <= 79.5; m++)
            {
                //a1=0.0063*Math.Pi/0.09
                a0 = 0;
                b0 = Math.PI / 2;
                en = Math.Abs(a - b);
                alfak = b - 0.05;
                i = 1;
 
                while (Math.Abs(en) > eps)
                {
                    alfak = alfak - f(alfak, m) / f1(alfak, m);
                    en = Math.Abs(alfak - b);
                    b = alfak;
                }
            } 
        }
 
        //трансцендентное уравнение
        private static double f(double alfak, double m)
        {
            return Math.Tan(alfak) - (Math.Abs(m) / alfak);
        }
 
        //его производная
        private static double f1(double alfak, double m)
        {
            return 1 / (Math.Pow(Math.Cos(alfak), 2)) + Math.Abs(m) / Math.Pow(alfak, 2);
        }
 
        //расчет суммы Ak и Fk
 
        static void Main()
        {
            int y = 0;
            int ksi = 1;
            int Rml = 5;
            double Tau = 0.09f;
            double a1=a*Math.PI/Tau;
            int j = 1;
            result=0;
 
            for (m = -79.5; m <= 79.5; m++)
            {
               gm=P(m)*Psi(m,ksi);
               lambda = (alfak * alfak + m * m) / (Math.Abs(m) * Math.Abs(gm));
               bk = Math.Sign(lambda) * rect;
               ak=(lambda*bk)/(lambda+j*Rml);
               fk=((2*alfak)/(2*alfak*a1+Math.Sin(2*alfak)*a1*Math.Abs(gm))*Math.Cos(alfak)*y);
               result=result+ak*fk;
               Console.WriteLine("Bm = {0}", result);
             }
             Console.ReadKey();
                Console.WriteLine();
        }
    }
}

Решение задачи: «Численное интегрирование(проверить код)»

textual
                            Листинг программы 
                        
 
                            
                        
class Program
    {
        public delegate double Function(double x);
        public static double Rectangular(Function f, double a, double b, int n)
        {
            double sum = 0.0;
            double h = (b - a) / n;
            for (int i = 0; i < n; i++)
            {
                sum += h * f(a + i * h);
            }
            return sum;
        }
 
        public static double Rectangular(double[] y, double a, double b, int n)
        {
            double sum = 0.0;
            double h = (b - a) / n;
            for (int i = 0; i < n; i++)
            {
                sum += h * y[i];
            }
            return sum;
        }
        static double f(double x)
        {
            return Math.Sin(x);
        }
        static double df(double x)
        {
            return Math.Cos(x);
        }
        static void Main(string[] args)
        {
            int n = 1000;
            double result;
            double a = 0.0;
            double b = 1.0;
 
            Console.WriteLine("\nTesting Rectangular method for integral(df(x)) = f(x) where");
            Console.WriteLine("f(x) = sin(x), df(x) = cos(x) and a=0 to b=1.\n");
 
            result = f(b) - f(a);
            Console.WriteLine("Analytic result         = " + result.ToString());
 
            result = Rectangular(df, a, b, n);
            Console.WriteLine("Result using function   = " + result.ToString());
 
            double[] y = new double[n];
            double h = (b - a) / (n - 1);
            for (int i = 0; i < n; i++)
            {
                double x = i * h;
                y[i] = df(x);
            }
            result = Rectangular(y, a, b, n);
            Console.WriteLine("Result using data array = " + result.ToString());
            Console.ReadLine();
        }
    }