Вычисляется последовательность значений модифицированной функции Бесселя I_a+n(x) для заданных значений a и x и для n=0,1,...,N.
Значения функции Бесселя вычисляются с заданной точностью и записываются в массив.

Структура:

Тип:	-	SUBROUTINE
Имена входа для пользователя:	-	DBESIN
Используемые внешние программы:	-	DGAMMA(C322)

Обращение:

CALL DBESIN(X,A,N,ND,BI), где:

X	-	(REAL*8) аргумент x;
A	-	(REAL*8) порядок a первой модифицированной функции Бесселя вычисляемой последовательности, 0 <= a < 1;
N	-	(INREGER) порядок a+N последней функции Бесселя вычисляемой последовательности;
ND	-	(INREGER) требуемое количество значащих десятичных цифр;
BI	-	(REAL*8) одномерный массив для записи результатов: BI(n+1) (n=0,...N) содержит Ia+n(x).

Ограничения:

Предполагается, что 0 <= a < 1, -100 <= N <= 100, x > 0.
Для отрицательных значений порядка ( N < 0 ) параметр a не может быть нулевым
( 0 < a < 1 ).

Ошибки исполнения:

Если какой-либо из параметров не находится в пределах
0 <= a < 1, -100 <= N <= 100, x > 0,
программа печатает сообщение
DBESIN ... UNREASONABLE ARGUMENT OR ORDER
и делает возврат в вызывающую программу.

Литература:

1. W.Gautshi. Bessel-functions of the first kind.
   Comm. ACM 7 (Aug. 1964), 479.
2. Л.М.Панченко и др., Библиотека программ на фортране, т.IV,
   Дубна, 1983 (подробное описание метода).

Пример:

       . . .
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION BI(3)
       X=5.0D0
       A=0.0D0
       CALL DBESIN(X,A,2,17,BI)
       DO 1 K=1,3
    1     WRITE(*,10) X,A+K-1,BI(K)
       . . .

Результат:

               X                    A+N                I(X)
       5.00000000000000000   .00000000000000000 27.23987182360444000
       5.00000000000000000  1.00000000000000000 24.33564214245052000
       5.00000000000000000  2.00000000000000000 17.50561496662423000

---