Locally-Hosted Collector/ numerical methods library for fortran 90

numerical methods library for fortran 90

!====================================================================
!================ ADTNUMF library. Version  20090204 ================
!====================================================================
! Copyright (C) 2007, 2009 Vsevolod Scherbinin
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
! 02110-1301 USA


Module Matrix

!=====================================================
!   Interface for subprograms which solve the linear algebraic equation array by
!   the Gauss method. The subprograms exist for the cases of double precision 
!   and double complex arrays. They may be called by command:
!   Call SlaeGauss(A,X)
!   where A - is the extended array of the equation array; X - the result vector.
!   Array A should be changed while subprogram works.
!=====================================================
Interface SlaeGauss
    Module Procedure SlaeGaussReal, SlaeGaussComp
End Interface
!=====================================================
!   Interface for subprograms which calculate the matrix determinant. The 
!   subprograms exist for the cases of double precision and double complex
!   arrays. They may be called by command:
!   Call Det(A, S)
!   where A - is the square array; S - the determinant.
!   Array A should be saved.
!=====================================================
Interface Det
    Module Procedure DetReal, DetComp
End Interface
!=====================================================
!   Interface for subprograms which find the inverse matrix
!   Call InvMatr(A)
!   where A - is the square array.
!   Array A should be changed while subprogram works.
!=====================================================
Interface InvMatr
    Module Procedure ReciprocationReal, ReciprocationComp
End Interface
!=====================================================
!   Interface for subprograms which find the inverse matrix
!   Call CondRat(A, Cond)
!   where A - is the square array (double precision or double complex, Cond - the
!   array condition number (double precision).
!   Array A should be saved while subprogram works.
!=====================================================
Interface CondRat
    Module Procedure CondNumbReal, CondNumbComp
End Interface

Contains
!=====================================================
!   Set of algebraic array solving subprogram for the Double Precision arguments.
!=====================================================
Subroutine SlaeGaussReal(A, X)
    Implicit None
    Double Precision :: S, R
    Integer :: FCnt, N1, J, SCnt, K1, N, I, Kst
    Double Precision, Dimension (:,:) :: A
    Double Precision, Dimension (:) :: X
    N=Size(X)
    N1=N+1
    !=====================================
    !   Forward trace
    !=====================================
    Do FCnt=1, N
        K1=FCnt+1
        S=A(FCnt,FCnt)
        J=FCnt
        Do SCnt=K1,N
            If (Abs(R).GT.Abs(S)) Then
                S=A(SCnt,FCnt); J=SCnt
            EndIf
        EndDo
        If (S.EQ.0._8) Return 
        If (J.NE.FCnt) Then
            Do SCnt=FCnt,N1
                R=A(FCnt,SCnt); A(FCnt,SCnt)=A(J,SCnt); A(J,SCnt)=R
            EndDo
        EndIf
        A(FCnt,K1:N1)=A(FCnt,K1:N1)/S
        Do SCnt=K1,N
            A(SCnt,K1:N1)= A(SCnt,K1:N1)- A(FCnt,K1:N1)*A(SCnt,FCnt)
        EndDo
    EndDo
    !=====================================
    !   Counter motion
    !=====================================
    X(N)=A(N,N1)
    Do Scnt=N-1,1,-1
        S=A(SCnt,N1)
        Do J=SCnt+1,N
            S=S-A(SCnt,J)*X(J)
        EndDo
        X(SCnt)=S
    EndDo
Return
End Subroutine SlaeGaussReal

!=====================================================
!   Set of algebraic array solving subprogram for the Double Complex arguments.
!=====================================================
Subroutine SlaeGaussComp(A,X)
    Implicit None
    Double Complex::S, R
    Integer::FCnt, N1, J, SCnt,K1,N, I,Kst
    Double Complex, Dimension (:,:) :: A
    Double Complex, Dimension (:) :: X
    N=Size(X)
    N1=N+1
    !=====================================
    !   Forward trace
    !=====================================
    Do FCnt=1, N
        K1=FCnt+1
        S=A(FCnt,FCnt)
        J=FCnt
        Do SCnt=K1,N
            If (Abs(R).GT.Abs(S)) Then
                S=A(SCnt,FCnt); J=SCnt
            EndIf
        EndDo
        If (S.EQ.(0._8,0._8)) Return 
        If (J.NE.FCnt) Then
            Do SCnt=FCnt,N1
                R=A(FCnt,SCnt); A(FCnt,SCnt)=A(J,SCnt); A(J,SCnt)=R
            EndDo
        EndIf
        A(FCnt,K1:N1)=A(FCnt,K1:N1)/S
        Do SCnt=K1,N
            A(SCnt,K1:N1)= A(SCnt,K1:N1)- A(FCnt,K1:N1)*A(SCnt,FCnt)
        EndDo
    EndDo
    !=====================================
    !   Counter motion
    !=====================================
    X(N)=A(N,N1)
    Do Scnt=N-1,1,-1
        S=A(SCnt,N1)
        Do J=SCnt+1,N
            S=S-A(SCnt,J)*X(J)
        EndDo
        X(SCnt)=S
    EndDo
    Return
End Subroutine SlaeGaussComp

!=====================================================
!   Matrix determinant computing for double precision matrixes.
!=====================================================
Subroutine DetReal(A, Det)
    Implicit None
    Double Precision::S, R, Det, P
    Integer :: Rank, I, K, N1, K1, J
    Double Precision, Dimension (:,:) :: A
    Double Precision, Allocatable, Dimension (:,:) :: Tmp
    Rank=Min(Size(A,1),Size(A,2))
    Allocate(Tmp(Rank,Rank))
    Tmp(:,:)=A(1:Rank,:)
    P=1._8
    Do K=1, Rank-1
        K1=K+1
        S=Tmp(K,K)
        J=K
        Do I=K1,Rank
            R=Tmp(I, K)
            If (Abs(R).GT.Abs(S)) Then
                S=R; J=I
            EndIf
        EndDo
        If (S.EQ.0._8) Return
        If (J.NE.K) Then
            Do I=K, Rank
                R=Tmp(K,I); Tmp(K,I)=Tmp(J,I); Tmp(J,I)=R
            EndDo
            P=-P
        EndIf
        Tmp(K,K1:Rank)=Tmp(K,K1:Rank)/S
        Do I=K1, Rank
            Tmp(I,K1:Rank)=Tmp(I,K1:Rank)-Tmp(K,K1:Rank)*Tmp(I,K)
        EndDo
        P=P*S
    EndDo
    S=P*Tmp(Rank,Rank)
    Det=S
    Deallocate(Tmp)
Return
End Subroutine DetReal

!=====================================================
!   Matrix determinant computing for double complex matrixes.
!=====================================================
Subroutine DetComp(A, Det)
    Implicit None
    Double Complex::S, R, Det, P
    Integer :: Rank, I, K, N1, K1, J
    Double Complex, Dimension (:,:) :: A
    Double Complex, Allocatable, Dimension (:,:) :: Tmp
    Rank=Min(Size(A,1),Size(A,2))
    Allocate(Tmp(Rank,Rank))
    Tmp(:,:)=A(1:Rank,:)
    P=1._8
    Do K=1, Rank-1
        K1=K+1
        S=Tmp(K,K)
        J=K
        Do I=K1,Rank
            R=Tmp(I, K)
            If (Abs(R).GT.Abs(S)) Then
                S=R; J=I
            EndIf
        EndDo
        If (S.EQ.(0._8,0._8)) Return
        If (J.NE.K) Then
            Do I=K, Rank
                R=Tmp(K,I); Tmp(K,I)=Tmp(J,I); Tmp(J,I)=R
            EndDo
            P=-P
        EndIf
        Tmp(K,K1:Rank)=Tmp(K,K1:Rank)/S
        Do I=K1, Rank
            Tmp(I,K1:Rank)=Tmp(I,K1:Rank)-Tmp(K,K1:Rank)*Tmp(I,K)
        EndDo
        P=P*S
    EndDo
    S=P*Tmp(Rank,Rank)
    Det=S
    Deallocate(Tmp)
Return
End Subroutine DetComp

!=====================================================
!   Matrix reciprocation for double precision matrixes. Matrix must be square.
!=====================================================
Subroutine ReciprocationReal(A)
    Implicit None
    Double Precision :: S, R, Det, P
    Integer :: Rank, Drank, I, J, K
    Double Precision, Dimension (:,:) :: A
    Double Precision, Allocatable, Dimension (:,:) :: Tmp
    Rank=Size(A,1)
    J=Size(A,2)
    If (Rank.NE.J) Then
        Write(*,*) 'InvMatr procedure error. The matrix is not square. Subprogram aborted.'
        Return
    EndIf
    If (Rank.EQ.1) Then
        A(1,1)=(1._8,0._8)/A(1,1)
        Return
    EndIf
    Drank=Rank*2
    Allocate(Tmp(Rank,Drank))
    !======================================
    !   Temporary array with initial and unity matrix creating
    !======================================
    Do I=1, Rank
        Tmp(I,1:Rank)=A(I,1:Rank)
        Tmp(I,Rank+1:Drank)=(0._8,0._8)
        Tmp(I,Rank+I)=(1.,0._8)
    EndDo
    !======================================
    !   Cycle by rows
    !======================================
    Do K=1, Rank
        S=Tmp(K,K); J=K
        If (K.NE.Rank) Then
            Do I=K+1,Rank
            If (Abs(Tmp(I,K)).GT.Abs(S)) Then
                S=Tmp(I,K); J=I
            EndIf
            EndDo
        EndIf
        If (Abs(S).EQ.(0._8,0._8)) Then
            Write(*,*) 'InvMatr procedure error. The matrix is degenerate (determinant is equal to zero).'
            Return
        EndIf
        If (J.NE.K) Then
            Do I=K, Drank
                R=Tmp(K,I); Tmp(K,I)=Tmp(J,I); Tmp(J,I)=R
            EndDo
        EndIf
        Tmp(K,K+1:Drank)=Tmp(K,K+1:Drank)/S
        If (K.NE.Rank) Then
            Do I=K+1, Rank
                Tmp(I, K+1:Drank)=Tmp(I, K+1:Drank)-Tmp(K, K+1:Drank)*Tmp(I,K)
            EndDo
        EndIf
        EndDo
        Do J=Rank+1, Drank
            Do I=Rank-1,1,-1
                R=Tmp(I,J)
                Do K=I+1,Rank
                    R=R-Tmp(K,J)*Tmp(I,K)
                EndDo
            Tmp(I,J)=R
            EndDo
        EndDo
        Do I=1,Rank
            A(I,1:Rank)=Tmp(I,Rank+1:Drank)
        EndDo
    Deallocate(Tmp)
Return
End Subroutine ReciprocationReal

!=====================================================
!   Matrix reciprocation for double complex matrixes. Matrix must be square.
!=====================================================
Subroutine ReciprocationComp(A)
    Implicit None
    Double Complex::S, R, Det, P
    Integer :: Rank, Drank, I, J, K
    Double Complex, Dimension (:,:) :: A
    Double Complex, Allocatable, Dimension (:,:) :: Tmp
    Rank=Size(A,1)
    J=Size(A,2)
    If (Rank.NE.J) Then
        Write(*,*) 'InvMatr procedure error. The matrix is not square. Subprogram aborted.'
        Return
    EndIf
    If (Rank.EQ.1) Then
        A(1,1)=(1._8,0._8)/A(1,1)
        Return
    EndIf
    Drank=Rank*2
    Allocate(Tmp(Rank,Drank))
    !======================================
    !   Temporary array with initial and unity matrix creating
    !======================================
    Do I=1, Rank
        Tmp(I,1:Rank)=A(I,1:Rank)
        Tmp(I,Rank+1:Drank)=(0._8,0._8)
        Tmp(I,Rank+I)=(1.,0._8)
    EndDo
    !======================================
    !   Cycle by rows
    !======================================
    Do K=1, Rank
        S=Tmp(K,K); J=K
        If (K.NE.Rank) Then
            Do I=K+1,Rank
            If (Abs(Tmp(I,K)).GT.Abs(S)) Then
                S=Tmp(I,K); J=I
            EndIf
            EndDo
        EndIf
        If (Abs(S).EQ.(0._8,0._8)) Then
            Write(*,*) 'InvMatr procedure error. The matrix is degenerate (determinant is equal to zero).'
            Return
        EndIf
        If (J.NE.K) Then
            Do I=K, Drank
                R=Tmp(K,I); Tmp(K,I)=Tmp(J,I); Tmp(J,I)=R
            EndDo
        EndIf
        Tmp(K,K+1:Drank)=Tmp(K,K+1:Drank)/S
        If (K.NE.Rank) Then
            Do I=K+1, Rank
                Tmp(I, K+1:Drank)=Tmp(I, K+1:Drank)-Tmp(K, K+1:Drank)*Tmp(I,K)
            EndDo
        EndIf
        EndDo
        Do J=Rank+1, Drank
            Do I=Rank-1,1,-1
                R=Tmp(I,J)
                Do K=I+1,Rank
                    R=R-Tmp(K,J)*Tmp(I,K)
                EndDo
            Tmp(I,J)=R
            EndDo
        EndDo
        Do I=1,Rank
            A(I,1:Rank)=Tmp(I,Rank+1:Drank)
        EndDo
    Deallocate(Tmp)
Return
End Subroutine ReciprocationComp

!=====================================================
!   Matrix condition number calculating for double precision matrixes. Matrix must be square.
!=====================================================
Subroutine CondNumbReal(A,Cn)
    Implicit None
    Double Precision::S, R, Cn, P, MNcop, MNinv
    Integer :: Rank, Drank, I, J, K
    Double Precision, Dimension (:,:) :: A
    Double Precision, Allocatable, Dimension (:,:) :: mCop, mInv
    Rank=Size(A,1)
    J=Size(A,2)
    If (Rank.NE.J) Then
        Write(*,*) 'The matrix is not square. Subprogram aborted.'
        Return
    EndIf
    Allocate(mCop(Rank,Rank))
    Allocate(mInv(Rank,Rank))
    mCop=A
    mInv=A
    Call InvMatr(mInv)
    MNcop=Sum(Abs(mCop)*Abs(mCop))  ! Initial matrix norm
    MNinv=Sum(Abs(mInv)*Abs(mInv))      !Inverse matrix norm
    Cn=Sqrt(MNcop*MNinv)
    Deallocate(mCop)
    Deallocate(mInv)
End Subroutine CondNumbReal

!=====================================================
!   Matrix condition number calculating for double precision matrixes. Matrix must be square.
!=====================================================
Subroutine CondNumbComp(A,Cn)
    Implicit None
    Double Precision:: Cn, MNcop, MNinv
    Integer :: Rank, J
    Double Complex, Dimension (:,:) :: A
    Double Complex, Allocatable, Dimension (:,:) :: mCop, mInv
    Rank=Size(A,1)
    J=Size(A,2)
    If (Rank.NE.J) Then
        Write(*,*) 'The matrix is not square. Subprogram aborted.'
        Return
    EndIf
    Allocate(mCop(Rank,Rank))
    Allocate(mInv(Rank,Rank))
    mCop=A
    mInv=A
    Call InvMatr(mInv)
    MNcop=Sum(Abs(mCop)*Abs(mCop))  ! Initial matrix norm
    MNinv=Sum(Abs(mInv)*Abs(mInv))      !Inverse matrix norm
    Cn=Sqrt(MNcop*MNinv)
    Deallocate(mCop)
    Deallocate(mInv)
End Subroutine CondNumbComp

End module Matrix
Module Integrals
!   Use Matrix
    Double Precision, Private:: WG16(16), XG16(16), WG40(40), XG40(40), WG96(96), XG96(96), XL16(16), XL40(40)
!=====================================================
!   Interface for subprograms which Gauss integral by 16 nodes
!   Call IntGauss16(A, B, N, I, F1, [F2], Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, F1 and F2 - is the function name, Dim - array dimension (optional).
!=====================================================
Interface IntGauss16
Module Procedure IntGauss16Real, IntGauss16RealN, IntGauss16Comp, IntGauss16CompN, &
    ArrIntGauss16Real, ArrIntGauss16RealN, ArrIntGauss16Comp, ArrIntGauss16CompN, &
    IntGaussProd16Real, IntGaussProd16RealN, IntGaussProd16Comp, IntGaussProd16CompN, &
    ArrIntGaussProd16Real, ArrIntGaussProd16RealN, ArrIntGaussProd16Comp, ArrIntGaussProd16CompN
End Interface IntGauss16


! Interface IntGauss16
! Module Procedure IntGauss16Real, IntGauss16RealN, IntGauss16Comp, IntGauss16CompN, &
!   ArrIntGauss16Real, ArrIntGauss16RealN, ArrIntGauss16Comp, ArrIntGauss16CompN
! End Interface IntGauss16
! Interface IntGaussProd16
! Module Procedure IntGaussProd16Real, IntGaussProd16RealN, IntGaussProd16Comp, IntGaussProd16CompN, &
!   ArrIntGaussProd16Real, ArrIntGaussProd16RealN, ArrIntGaussProd16Comp, ArrIntGaussProd16CompN
! End Interface IntGaussProd16

!=====================================================
!   Interface for subprograms which Gauss integral by 40 nodes
!   Call IntGauss40(A, B, N, I, E, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, F - is the function name, Dim - array dimension (optional).
!=====================================================
Interface IntGauss40
Module Procedure IntGauss40Real, IntGauss40RealN, IntGauss40Comp, IntGauss40CompN, &
    ArrIntGauss40Real, ArrIntGauss40RealN, ArrIntGauss40Comp, ArrIntGauss40CompN
End Interface IntGauss40
Interface IntGauss
Module Procedure IntGauss40Real, IntGauss40RealN, IntGauss40Comp, IntGauss40CompN, &
    ArrIntGauss40Real, ArrIntGauss40RealN, ArrIntGauss40Comp, ArrIntGauss40CompN
End Interface IntGauss
!=====================================================
!   Interface for subprograms which Gauss integral by 96 nodes
!   Call IntGauss96(A, B, N, I, E, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, F - is the function name, Dim - array dimension (optional).
!=====================================================
Interface IntGauss96
Module Procedure IntGauss96Real, IntGauss96RealN, IntGauss96Comp, IntGauss96CompN, &
    ArrIntGauss96Real, ArrIntGauss96RealN, ArrIntGauss96Comp, ArrIntGauss96CompN
End Interface IntGauss96
!=====================================================
!   Interface for subprograms which Gauss-Kronrod integral by 15 nodes
!   Call IntGaussKronrod15(A, B, N, I, E, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, E - is the error upper bound, F - is the function name, Dim - array
!   dimension (optional).
!=====================================================
Interface IntGaussKronrod15
Module Procedure IntGaussKron15Real, IntGaussKron15RealN, IntGaussKron15Comp, IntGaussKron15CompN, &
ArrIntGaussKron15Real, ArrIntGaussKron15RealN, ArrIntGaussKron15Comp, ArrIntGaussKron15CompN
End Interface IntGaussKronrod15
!=====================================================
!   Interface for subprograms which Gauss-Kronrod integral by 15 nodes
!   Call IntGaussKronrod41(A, B, N, I, E, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, E - is the error upper bound, F - is the function name, Dim - array
!   dimension (optional).
!=====================================================
Interface IntGaussKronrod41
Module Procedure IntGaussKron41Real, IntGaussKron41RealN, IntGaussKron41Comp, IntGaussKron41CompN, &
ArrIntGaussKron41Real, ArrIntGaussKron41RealN, ArrIntGaussKron41Comp, ArrIntGaussKron41CompN
End Interface IntGaussKronrod41
!=====================================================
!   Interface for subprograms which Gauss-Kronrod integral by 15 nodes
!   Call IntGaussKronrod61(A, B, N, I, E, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, E - is the error upper bound, F - is the function name, Dim - array
!   dimension (optional).
!=====================================================
Interface IntGaussKronrod61
Module Procedure IntGaussKron61Real, IntGaussKron61RealN, IntGaussKron61Comp, IntGaussKron61CompN, &
ArrIntGaussKron61Real, ArrIntGaussKron61RealN, ArrIntGaussKron61Comp, ArrIntGaussKron61CompN
End Interface IntGaussKronrod61
!=====================================================
!   Interface for subprograms which trapezoidal Filon method for highly oscillation
!   integrals computing
!   Call IntTrapezFilon(A, B, N, I, E, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, E - is the error upper bound, F - is the function name, Dim - array
!   dimension (optional).
!=====================================================
Interface IntTrapezFilon
Module Procedure IntFilonTrap, ArrIntFilonTrap
End Interface IntTrapezFilon
!=====================================================
!   Interface for subprograms which 16-nodes Iserles method for highly oscillation
!   integrals computing
!   Call IntIserles16(A, B, N, I, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, E - is the error upper bound, F - is the function name, Dim - array
!   dimension (optional).
!=====================================================
Interface IntIserles16
Module Procedure IntIserles16Comp, IntIserles16CompN, ArrIntIserles16Comp, ArrIntIserles16CompN
End Interface IntIserles16
!=====================================================
!   Interface for subprograms which 40-nodes Iserles method for highly oscillation
!   integrals computing
!   Call IntIserles40(A, B, N, I, F, Dim)
!   where A - is the low limit, B - is the high limit, N - is the number of discretes (optional),
!   I - is the result, E - is the error upper bound, F - is the function name, Dim - array
!   dimension (optional).
!=====================================================
Interface IntIserles40
Module Procedure IntIserles40Comp, IntIserles40CompN, ArrIntIserles40Comp, ArrIntIserles40CompN
End Interface IntIserles40
!=====================================================
!   Interface for subprograms which make Simpson-type 5-point cubature
!   Call IntFiveCubature(A, B, C, D, NX, NY, I, F, Dim)
!   where A - is the X axis low limit, B - is the X axis  high limit, C - is the Y axis low limit,
!   D - is the Y axis  high limit, NX - is the number of discretes by X, NY - is the number
!   of discretes by X, I - is the result, F - is the function name, Dim - array dimension
!   (optional).
!=====================================================
Interface IntFiveCubature
Module Procedure IntCubature5Real, ArrIntCubature5Real, IntCubature5Comp, ArrIntCubature5Comp
End Interface IntFiveCubature
!=====================================================
!   Interface for subprograms which make Gauss-type 40-node cubature
!   Call IntGaussCubature(A, B, C, D, NX, NY, I, F, Dim)
!   where A - is the X axis low limit, B - is the X axis  high limit, C - is the Y axis low limit,
!   D - is the Y axis  high limit, NX - is the number of discretes by X, NY - is the number
!   of discretes by X, I - is the result, F - is the function name, Dim - array dimension
!   (optional).
!=====================================================
Interface IntGaussCubature
Module Procedure IntGauss40CubatureReal, IntGauss40CubatureRealN, ArrIntGauss40CubatureReal, &
ArrIntGauss40CubatureRealN, IntGauss40CubatureComp, IntGauss40CubatureCompN, &
ArrIntGauss40CubatureComp, ArrIntGauss40CubatureCompN
End Interface IntGaussCubature

!=====================================================
!   Nodes and weights
!=====================================================
    Data WG16 /0.027152459411754094852_8, 0.062253523938647892863_8, 0.095158511682492784810_8, &
        0.124628971255533872052_8, 0.149595988816576732081_8, 0.169156519395002538189_8, &
        0.182603415044923588867_8, 0.189450610455068496285_8, 0.189450610455068496285_8, &
        0.182603415044923588867_8, 0.169156519395002538189_8, 0.149595988816576732081_8, &
        0.124628971255533872052_8, 0.095158511682492784810_8, 0.062253523938647892863_8, &
        0.027152459411754094852_8/
    Data XG16 /0.989400934991649923596_8, 0.944575023073232576078_8, 0.865631202387831743880_8, &
        0.755404408355003033895_8, 0.617876244402643748447_8, 0.458016777657227386342_8, &
        0.281603550779258913230_8, 0.095012509837637440185_8, -0.095012509837637440185_8, &
        -0.281603550779258913230_8, -0.458016777657227386342_8, -0.617876244402643748447_8, &
        -0.755404408355003033895_8, -0.865631202387831743880_8, -0.944575023073232576078_8, &
        -0.989400934991649923596_8/
    Data WG40 /0.004521277098533191258_8, 0.010498284531152813615_8, 0.016421058381907888713_8, &
        0.022245849194166957262_8, 0.027937006980023401098_8, 0.033460195282547847393_8, &
        0.038782167974472017640_8, 0.043870908185673271992_8, 0.048695807635072232061_8, &
        0.053227846983936824355_8, 0.057439769099391551367_8, 0.061306242492928939167_8, &
        0.064804013456601038075_8, 0.067912045815233903826_8, 0.070611647391286779695_8, &
        0.072886582395804059061_8, 0.074723169057968264200_8, 0.076110361900626242372_8, &
        0.077039818164247965588_8, 0.077505947978424811264_8, 0.077505947978424811264_8, &
        0.077039818164247965588_8, 0.076110361900626242372_8, 0.074723169057968264200_8, &
        0.072886582395804059061_8, 0.070611647391286779695_8, 0.067912045815233903826_8, &
        0.064804013456601038075_8, 0.061306242492928939167_8, 0.057439769099391551367_8, &
        0.053227846983936824355_8, 0.048695807635072232061_8, 0.043870908185673271992_8, &
        0.038782167974472017640_8, 0.033460195282547847393_8, 0.027937006980023401098_8, &
        0.022245849194166957262_8, 0.016421058381907888713_8, 0.010498284531152813615_8, &
        0.004521277098533191258_8/
    Data XG40 /0.998237709710559200350_8, 0.990726238699457006453_8, 0.977259949983774262663_8, &
        0.957916819213791655805_8, 0.932812808278676533361_8, 0.902098806968874296728_8, &
        0.865959503212259503821_8, 0.824612230833311663196_8, 0.778305651426519387695_8, &
        0.727318255189927103281_8, 0.671956684614179548379_8, 0.612553889667980237953_8, &
        0.549467125095128202076_8, 0.483075801686178712909_8, 0.413779204371605001525_8, &
        0.341994090825758473007_8, 0.268152185007253681141_8, 0.192697580701371099716_8, &
        0.116084070675255208483_8, 0.038772417506050821933_8, -0.038772417506050821933_8, &
        -0.116084070675255208483_8, -0.192697580701371099716_8, -0.268152185007253681141_8, &
        -0.341994090825758473007_8, -0.413779204371605001525_8, -0.483075801686178712909_8, &
        -0.549467125095128202076_8, -0.612553889667980237953_8, -0.671956684614179548379_8, &
        -0.727318255189927103281_8, -0.778305651426519387695_8, -0.824612230833311663196_8, &
        -0.865959503212259503821_8, -0.902098806968874296728_8, -0.932812808278676533361_8, &
        -0.957916819213791655805_8, -0.977259949983774262663_8, -0.990726238699457006453_8, &
        -0.998237709710559200350_8/
    Data WG96 /0.000796792065552012429_8, 0.001853960788946921732_8, 0.002910731817934946408_8, &
        0.003964554338444686674_8, 0.005014202742927517693_8, 0.006058545504235961683_8, &
        0.007096470791153865269_8, 0.008126876925698759217_8, 0.009148671230783386633_8, &
        0.010160770535008415758_8, 0.011162102099838498591_8, 0.012151604671088319635_8, &
        0.013128229566961572637_8, 0.014090941772314860916_8, 0.015038721026994938006_8, &
        0.015970562902562291381_8, 0.016885479864245172450_8, 0.017782502316045260838_8, &
        0.018660679627411467385_8, 0.019519081140145022410_8, 0.020356797154333324595_8, &
        0.021172939892191298988_8, 0.021966644438744349195_8, 0.022737069658329374001_8, &
        0.023483399085926219842_8, 0.024204841792364691282_8, 0.024900633222483610288_8, &
        0.025570036005349361499_8, 0.026212340735672413913_8, 0.026826866725591762198_8, &
        0.027412962726029242823_8, 0.027970007616848334440_8, 0.028497411065085385646_8, &
        0.028994614150555236543_8, 0.029461089958167905970_8, 0.029896344136328385984_8, &
        0.030299915420827593794_8, 0.030671376123669149014_8, 0.031010332586313837423_8, &
        0.031316425596861355813_8, 0.031589330770727168558_8, 0.031828758894411006535_8, &
        0.032034456231992663218_8, 0.032206204794030250669_8, 0.032343822568575928429_8, &
        0.032447163714064269364_8, 0.032516118713868835987_8, 0.032550614492363166242_8, &
        0.032550614492363166242_8, 0.032516118713868835987_8, 0.032447163714064269364_8, &
        0.032343822568575928429_8, 0.032206204794030250669_8, 0.032034456231992663218_8, &
        0.031828758894411006535_8, 0.031589330770727168558_8, 0.031316425596861355813_8, &
        0.031010332586313837423_8, 0.030671376123669149014_8, 0.030299915420827593794_8, &
        0.029896344136328385984_8, 0.029461089958167905970_8, 0.028994614150555236543_8, &
        0.028497411065085385646_8, 0.027970007616848334440_8, 0.027412962726029242823_8, &
        0.026826866725591762198_8, 0.026212340735672413913_8, 0.025570036005349361499_8, &
        0.024900633222483610288_8, 0.024204841792364691282_8, 0.023483399085926219842_8, &
        0.022737069658329374001_8, 0.021966644438744349195_8, 0.021172939892191298988_8, &
        0.020356797154333324595_8, 0.019519081140145022410_8, 0.018660679627411467385_8, &
        0.017782502316045260838_8, 0.016885479864245172450_8, 0.015970562902562291381_8, &
        0.015038721026994938006_8, 0.014090941772314860916_8, 0.013128229566961572637_8, &
        0.012151604671088319635_8, 0.011162102099838498591_8, 0.010160770535008415758_8, &
        0.009148671230783386633_8, 0.008126876925698759217_8, 0.007096470791153865269_8, &
        0.006058545504235961683_8, 0.005014202742927517693_8, 0.003964554338444686674_8, &
        0.002910731817934946408_8, 0.001853960788946921732_8, 0.000796792065552012429_8/
    Data XG96 /0.999689503883230766828_8, 0.998364375863181677724_8, 0.995981842987209290650_8, &
        0.992543900323762624572_8, 0.988054126329623799481_8, 0.982517263563014677447_8, &
        0.975939174585136466453_8, 0.968326828463264212174_8, 0.959688291448742539300_8, &
        0.950032717784437635756_8, 0.939370339752755216932_8, 0.927712456722308690965_8, &
        0.915071423120898074206_8, 0.901460635315852341319_8, 0.886894517402420416057_8, &
        0.871388505909296502874_8, 0.854959033434601455463_8, 0.837623511228187121494_8, &
        0.819400310737931675539_8, 0.800308744139140817229_8, 0.780369043867433217604_8, &
        0.759602341176647498703_8, 0.738030643744400132851_8, 0.715676812348967626225_8, &
        0.692564536642171561344_8, 0.668718310043916153953_8, 0.644163403784967106798_8, &
        0.618925840125468570386_8, 0.593032364777572080684_8, 0.566510418561397168404_8, &
        0.539388108324357436227_8, 0.511694177154667673586_8, 0.483457973920596359768_8, &
        0.454709422167743008636_8, 0.425478988407300545365_8, 0.395797649828908603285_8, &
        0.365696861472313635031_8, 0.335208522892625422616_8, 0.304364944354496353024_8, &
        0.273198812591049141487_8, 0.241743156163840012328_8, 0.210031310460567203603_8, &
        0.178096882367618602759_8, 0.145973714654896941989_8, 0.113695850110665920911_8, &
        0.081297495464425558994_8, 0.048812985136049731112_8, 0.016276744849602969579_8, &
        -0.016276744849602969579_8, -0.048812985136049731112_8, -0.081297495464425558994_8, &
        -0.113695850110665920911_8, -0.145973714654896941989_8, -0.178096882367618602759_8, &
        -0.210031310460567203603_8, -0.241743156163840012328_8, -0.273198812591049141487_8, &
        -0.304364944354496353024_8, -0.335208522892625422616_8, -0.365696861472313635031_8, &
        -0.395797649828908603285_8, -0.425478988407300545365_8, -0.454709422167743008636_8, &
        -0.483457973920596359768_8, -0.511694177154667673586_8, -0.539388108324357436227_8, &
        -0.566510418561397168404_8, -0.593032364777572080684_8, -0.618925840125468570386_8, &
        -0.644163403784967106798_8, -0.668718310043916153953_8, -0.692564536642171561344_8, &
        -0.715676812348967626225_8, -0.738030643744400132851_8, -0.759602341176647498703_8, &
        -0.780369043867433217604_8, -0.800308744139140817229_8, -0.819400310737931675539_8, &
        -0.837623511228187121494_8, -0.854959033434601455463_8, -0.871388505909296502874_8, &
        -0.886894517402420416057_8, -0.901460635315852341319_8, -0.915071423120898074206_8, &
        -0.927712456722308690965_8, -0.939370339752755216932_8, -0.950032717784437635756_8, &
        -0.959688291448742539300_8, -0.968326828463264212174_8, -0.975939174585136466453_8, &
        -0.982517263563014677447_8, -0.988054126329623799481_8, -0.992543900323762624572_8, &
        -0.995981842987209290650_8, -0.998364375863181677724_8, -0.999689503883230766828_8/
    Data XL16 /1._8, 0.98478402313511_8, 0.94960026654674_8, 0.89600414593091_8, 0.82619435144125_8, &
        0.74302971094357_8, 0.64991523445038_8, 0.55066313676097_8, 0.44933686323903_8, &
        0.35008476554962_8, 0.25697028905643_8, 0.17380564855875_8, 0.10399585406909_8, &
        0.05039973345326_8, 0.01521597686489_8, 0._8/
    Data XL40 /1._8, 0.99764896462217_8, 0.99213314035875_8, 0.98350503824399_8, 0.97181988247180_8, &
        0.95715166984510_8, 0.93959317173967_8, 0.91925541138905_8, 0.89626697630078_8, &
        0.87077320957369_8, 0.84293529254216_8, 0.81292922637763_8, 0.78094471964736_8, &
        0.74718398906263_8, 0.71186048107778_8, 0.67519752245709_8, 0.63742690835716_8, &
        0.59878743685946_8, 0.55952339922249_8, 0.51988303540109_8, 0.48011696459891_8, &
        0.44047660077751_8, 0.40121256314054_8, 0.36257309164284_8, 0.32480247754291_8, &
        0.28813951892222_8, 0.25281601093737_8, 0.21905528035264_8, 0.18707077362237_8, &
        0.15706470745784_8, 0.12922679042631_8, 0.10373302369922_8, 0.08074458861095_8, &
        0.06040682826033_8, 0.04284833015490_8, 0.02818011752820_8, 0.01649496175601_8, &
        0.00786685964125_8, 0.00235103537783_8, 0._8 /
Contains
!=====================================================
!   Gauss method of 16 order integration for double precision functions.
!=====================================================
Recursive Subroutine IntGauss16Real(Alim, Blim, PResult, FuncS)
    Implicit None
    Double Precision::  Alim, Blim, Plength
    Double Precision:: Y, PResult, FuncS
    Integer:: Kstep, Lstep
    External FuncS
    PResult=0._8
    Do Kstep=1,16
        Y=0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        PResult=PResult+WG16(Kstep)*FuncS(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine IntGauss16Real

Recursive Subroutine IntGauss16RealN(Llim, Hlim, Razb, LResult, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Precision:: PResult, LResult, FuncN
    Integer:: Kstep, Lstep, Razb
    External FuncN
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGauss16Real(Alim, Blim, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGauss16Real(Alim, Blim, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntGauss16RealN

!=====================================================
!   Gauss method of 16 order integration for double precision array functions.
!=====================================================
Recursive Subroutine ArrIntGauss16Real(Alim, Blim, PResult, FuncS, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength
    Double Precision:: Y
    Integer:: Kstep, Lstep, Razm
    Double Precision, Dimension(Razm):: PResult
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    PResult=0._8
    Do Kstep=1,16
        Y=0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        PResult=PResult+WG16(Kstep)*FuncS(Y, Razm)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine ArrIntGauss16Real

Recursive Subroutine ArrIntGauss16RealN(Llim, Hlim, Razb, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Precision, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGauss16Real(Alim, Blim, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGauss16Real(Alim, Blim, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGauss16RealN

!=====================================================
!   Gauss method of 16 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntGauss16Comp(Alim, Blim, PResult, FuncS)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Double Complex:: PResult, FuncS
    Integer:: Kstep, Lstep, Razb
    External FuncS
    PResult=(0._8,0._8)
    Do Kstep=1,16
        Y=0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        PResult=PResult+WG16(Kstep)*FuncS(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine IntGauss16Comp

Recursive Subroutine IntGauss16CompN(Llim, Hlim, Razb, LResult,FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Complex:: PResult, LResult, FuncN
    Integer:: Kstep, Lstep, Razb
    External FuncN
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGauss16Comp(Alim, Blim, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGauss16Comp(Alim, Blim, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntGauss16CompN

!=====================================================
!   Gauss method of 16 order integration for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntGauss16Comp(Alim, Blim, PResult, FuncS, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Integer:: Kstep, Lstep, Razm
    Double Complex, Dimension(Razm):: PResult, Tmp
    Interface
        Function FuncS(Z, N)
            Integer:: N
            Double Precision:: Z
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    PResult(1:Razm)=(0._8,0._8)
    Do Kstep=1,16
        Y = 0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        Tmp = FuncS(Y, Razm)
        PResult = PResult + WG16(Kstep)*Tmp
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine ArrIntGauss16Comp

Recursive Subroutine ArrIntGauss16CompN(Llim, Hlim, Razb, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Dcmplx(Real(Razb),0._8)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGauss16Comp(Alim, Blim, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGauss16Comp(Alim, Blim, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGauss16CompN

!=====================================================
!   Gauss method of 16 order integration for double precision function product.
!   Интегрирование методом Гаусса 16 порядка произведения двух функций.
!=====================================================
Recursive Subroutine IntGaussProd16Real(Alim, Blim, PResult, FuncS1, FuncS2)
    Implicit None
    Double Precision::  Alim, Blim, Plength
    Double Precision:: Y, PResult, FuncS1, FuncS2
    Integer:: Kstep, Lstep
    External FuncS1, FuncS2
    PResult=0._8
    Do Kstep=1,16
        Y=0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        PResult=PResult+WG16(Kstep)*FuncS1(Y)*FuncS2(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine IntGaussProd16Real

Recursive Subroutine IntGaussProd16RealN(Llim, Hlim, Razb, LResult, FuncN1, FuncN2)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Precision:: PResult, LResult, FuncN1, FuncN2
    Integer:: Kstep, Lstep, Razb
    External FuncN1, FuncN2
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussProd16Real(Alim, Blim, PResult, FuncN1, FuncN2)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussProd16Real(Alim, Blim, PResult, FuncN1, FuncN2)
    LResult=LResult+PResult
Return
End Subroutine IntGaussProd16RealN

!=====================================================
!   Gauss method of 16 order integration for double precision array functions products.
!=====================================================
Recursive Subroutine ArrIntGaussProd16Real(Alim, Blim, PResult, FuncS1, FuncS2, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength
    Double Precision:: Y
    Integer:: Kstep, Lstep, Razm
    Double Precision, Dimension(Razm):: PResult
    Interface
        Function FuncS1(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS1
        End Function FuncS1
        Function FuncS2(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS2
        End Function FuncS2
    End Interface
    PResult=0._8
    Do Kstep=1,16
        Y=0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        PResult=PResult+WG16(Kstep)*FuncS1(Y, Razm)*FuncS2(Y, Razm)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine ArrIntGaussProd16Real

Recursive Subroutine ArrIntGaussProd16RealN(Llim, Hlim, Razb, LResult, FuncN1, FuncN2, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Precision, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN1(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN1
        End Function FuncN1
        Function FuncN2(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN2
        End Function FuncN2
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussProd16Real(Alim, Blim, PResult, FuncN1, FuncN2, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussProd16Real(Alim, Blim, PResult, FuncN1, FuncN2, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGaussProd16RealN

!=====================================================
!   Gauss method of 16 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntGaussProd16Comp(Alim, Blim, PResult, FuncS1, FuncS2)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Double Complex:: PResult, FuncS1, FuncS2
    Integer:: Kstep, Lstep, Razb
    External FuncS1, FuncS2
    PResult=(0._8,0._8)
    Do Kstep=1,16
        Y=0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        PResult=PResult+WG16(Kstep)*FuncS1(Y)*FuncS2(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine IntGaussProd16Comp

Recursive Subroutine IntGaussProd16CompN(Llim, Hlim, Razb, LResult, FuncN1, FuncN2)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Complex:: PResult, LResult, FuncN1, FuncN2
    Integer:: Kstep, Lstep, Razb
    External FuncN1, FuncN2
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussProd16Comp(Alim, Blim, PResult, FuncN1, FuncN2)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussProd16Comp(Alim, Blim, PResult, FuncN1, FuncN2)
    LResult=LResult+PResult
Return
End Subroutine IntGaussProd16CompN

!=====================================================
!   Gauss method of 16 order integration for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntGaussProd16Comp(Alim, Blim, PResult, FuncS1, FuncS2, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Integer:: Kstep, Lstep, Razm
    Double Complex, Dimension(Razm):: PResult, Tmp
    Interface
        Function FuncS1(Z, N)
            Integer:: N
            Double Precision:: Z
            Double Complex, Dimension(N):: FuncS1
        End Function FuncS1
        Function FuncS2(Z, N)
            Integer:: N
            Double Precision:: Z
            Double Complex, Dimension(N):: FuncS2
        End Function FuncS2
    End Interface
    PResult(1:Razm)=(0._8,0._8)
    Do Kstep=1,16
        Y = 0.5_8*((Blim-Alim)*XG16(Kstep)+Blim+Alim)
        Tmp = FuncS1(Y, Razm)*FuncS2(Y, Razm)
        PResult = PResult + WG16(Kstep)*Tmp
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine ArrIntGaussProd16Comp

Recursive Subroutine ArrIntGaussProd16CompN(Llim, Hlim, Razb, LResult, FuncN1, FuncN2, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN1(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN1
        End Function FuncN1
        Function FuncN2(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN2
        End Function FuncN2
    End Interface
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Dcmplx(Real(Razb),0._8)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussProd16Comp(Alim, Blim, PResult, FuncN1, FuncN2, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussProd16Comp(Alim, Blim, PResult, FuncN1, FuncN2, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGaussProd16CompN


!=====================================================
!   Gauss method of 40 order integration for double precision functions.
!=====================================================
Recursive Subroutine IntGauss40Real(Alim, Blim, PResult, FuncS)
    Implicit None
    Double Precision:: Alim, Blim, Plength
    Double Precision:: Y, PResult, FuncS
    Integer:: Kstep, Lstep
    External FuncS
    PResult=0._8
    Do Kstep=1,40
        Y=0.5_8*((Blim-Alim)*XG40(Kstep)+Blim+Alim)
        PResult=PResult+WG40(Kstep)*FuncS(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine IntGauss40Real

Recursive Subroutine IntGauss40RealN(Llim, Hlim, Razb, LResult, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Precision:: Y, PResult, LResult, FuncN
    Integer:: Kstep, Lstep, Razb
    External FuncN
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGauss40Real(Alim, Blim, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGauss40Real(Alim, Blim, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntGauss40RealN

!=====================================================
!   Gauss method of 40 order integration for double precision array functions.
!=====================================================
Recursive Subroutine ArrIntGauss40Real(Alim, Blim, PResult, FuncS, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength
    Double Precision:: Y
    Integer:: Kstep, Lstep, Razm
    Double Precision, Dimension(Razm):: PResult
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    PResult=0._8
    Do Kstep=1,40
        Y=0.5_8*((Blim-Alim)*XG40(Kstep)+Blim+Alim)
        PResult=PResult+WG40(Kstep)*FuncS(Y, Razm)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine ArrIntGauss40Real

Recursive Subroutine ArrIntGauss40RealN(Llim, Hlim, Razb, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Precision, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGauss40Real(Alim, Blim, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGauss40Real(Alim, Blim, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGauss40RealN

!=====================================================
!   Gauss method of 40 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntGauss40Comp(Alim, Blim, PResult, FuncS)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Double Complex:: PResult, FuncS
    Integer:: Kstep, Lstep
    External FuncS
    PResult=(0._8,0._8)
    Do Kstep=1,40
        Y=0.5_8*((Blim-Alim)*XG40(Kstep)+Blim+Alim)
        PResult=PResult+WG40(Kstep)*FuncS(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine IntGauss40Comp

Recursive Subroutine IntGauss40CompN(Llim, Hlim, Razb,LResult,FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Complex:: PResult, LResult, FuncN
    Integer:: Kstep, Lstep, Razb
    External FuncN
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGauss40Comp(Alim, Blim, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGauss16Comp(Alim, Blim, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntGauss40CompN

!=====================================================
!   Gauss method of 40 order integration for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntGauss40Comp(Alim, Blim, PResult, FuncS, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Integer:: Kstep, Lstep, Razm
    Double Complex, Dimension(Razm):: PResult
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    PResult=(0._8,0._8)
    Do Kstep=1,40
        Y=0.5_8*((Blim-Alim)*XG40(Kstep)+Blim+Alim)
        PResult=PResult+WG40(Kstep)*FuncS(Y, Razm)
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine ArrIntGauss40Comp

Recursive Subroutine ArrIntGauss40CompN(Llim, Hlim, Razb, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Dcmplx(Real(Razb),0._8)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGauss40Comp(Alim, Blim, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGauss40Comp(Alim, Blim, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGauss40CompN

!=====================================================
!   Gauss method of 96 order integration for double precision functions.
!=====================================================
Recursive Subroutine IntGauss96Real(Alim, Blim, PResult, FuncS)
    Implicit None
    Double Precision:: Alim, Blim, Plength
    Double Precision:: Y, PResult, FuncS
    Integer:: Kstep, Lstep
    External FuncS
    PResult=0._8
    Do Kstep=1,96
        Y=0.5_8*((Blim-Alim)*XG96(Kstep)+Blim+Alim)
        PResult=PResult+WG96(Kstep)*FuncS(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine IntGauss96Real

Recursive Subroutine IntGauss96RealN(Llim, Hlim, Razb, LResult, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Precision:: Y, PResult, LResult, FuncN
    Integer:: Kstep, Lstep, Razb
    External FuncN
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGauss96Real(Alim, Blim, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGauss96Real(Alim, Blim, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntGauss96RealN

!=====================================================
!   Gauss method of 96 order integration for double precision array functions.
!=====================================================
Recursive Subroutine ArrIntGauss96Real(Alim, Blim, PResult, FuncS, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength
    Double Precision:: Y
    Integer:: Kstep, Lstep, Razm
    Double Precision, Dimension(Razm):: PResult
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    PResult=0._8
    Do Kstep=1,96
        Y=0.5_8*((Blim-Alim)*XG96(Kstep)+Blim+Alim)
        PResult=PResult+WG96(Kstep)*FuncS(Y, Razm)
    EndDo
    PResult=PResult*(Blim-Alim)*0.5_8
Return
End Subroutine ArrIntGauss96Real

Recursive Subroutine ArrIntGauss96RealN(Llim, Hlim, Razb, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Precision, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGauss96Real(Alim, Blim, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGauss96Real(Alim, Blim, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGauss96RealN

!=====================================================
!   Gauss method of 96 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntGauss96Comp(Alim, Blim, PResult, FuncS)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Double Complex:: PResult, FuncS
    Integer:: Kstep, Lstep
    External FuncS
    PResult=(0._8,0._8)
    Do Kstep=1,96
        Y=0.5_8*((Blim-Alim)*XG96(Kstep)+Blim+Alim)
        PResult=PResult+WG96(Kstep)*FuncS(Y)
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine IntGauss96Comp

Recursive Subroutine IntGauss96CompN(Llim, Hlim, Razb,LResult, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Double Complex:: PResult, LResult, FuncN
    Integer:: Kstep, Lstep, Razb
    External FuncN
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGauss96Comp(Alim, Blim, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGauss16Comp(Alim, Blim, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntGauss96CompN

!=====================================================
!   Gauss method of 96 order integration for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntGauss96Comp(Alim, Blim, PResult, FuncS, Razm)
    Implicit None
    Double Precision:: Alim, Blim, Plength, Y
    Integer:: Kstep, Lstep, Razm
    Double Complex, Dimension(Razm):: PResult
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    PResult=(0._8,0._8)
    Do Kstep=1,96
        Y=0.5_8*((Blim-Alim)*XG96(Kstep)+Blim+Alim)
        PResult=PResult+WG96(Kstep)*FuncS(Y, Razm)
    EndDo
    PResult=PResult*(Blim-Alim)*(0.5_8,0._8)
Return
End Subroutine ArrIntGauss96Comp

Recursive Subroutine ArrIntGauss96CompN(Llim, Hlim, Razb, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Kstep, Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=(0._8,0._8)
    Plength=(Hlim-Llim)/Dcmplx(Real(Razb),0._8)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGauss96Comp(Alim, Blim, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGauss96Comp(Alim, Blim, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntGauss96CompN

!=====================================================
!   Gauss-Kronrod method of 15 order integration for double precision functions.
!=====================================================
Recursive Subroutine IntGaussKron15Real(Alim, Blim, IK, Error, FuncS)
    Implicit None
    Double Precision:: WG(7), WK(15), X(15), Alim, Blim, Error
    Double Precision:: Y, IG, IK, FuncS, C1, C2
    Integer:: Kstep, Lstep, I
    External FuncS
    Data WG(1:2) /0.129484966168869693270611432679082_8, 0.279705391489276667901467771423780_8/
    Data WG(3:4) /0.381830050505118944950369775488975_8, 0.417959183673469387755102040816327_8/
    Data WG(5:6) /0.381830050505118944950369775488975_8, 0.279705391489276667901467771423780_8/
    Data WG(7) /0.129484966168869693270611432679082_8/
    Data WK(1:2) /0.022935322010529224963732008058970_8, 0.063092092629978553290700663189204_8/
    Data WK(3:4) /0.104790010322250183839876322541518_8, 0.140653259715525918745189590510238_8/
    Data WK(5:6) /0.169004726639267902826583426598550_8, 0.190350578064785409913256402421014_8/
    Data WK(7:8) /0.204432940075298892414161999234649_8, 0.209482141084727828012999174891714_8/
    Data WK(9:10) /0.204432940075298892414161999234649_8, 0.190350578064785409913256402421014_8/
    Data WK(11:12) /0.169004726639267902826583426598550_8, 0.140653259715525918745189590510238_8/
    Data WK(13:14) /0.104790010322250183839876322541518_8, 0.063092092629978553290700663189204_8/
    Data WK(15) /0.022935322010529224963732008058970_8/
    Data X(1:2) /0.991455371120812639206854697526329_8, 0.949107912342758524526189684047851_8/
    Data X(3:4) /0.864864423359769072789712788640926_8, 0.741531185599394439863864773280788_8/
    Data X(5:6) /0.586087235467691130294144838258730_8, 0.405845151377397166906606412076961_8/
    Data X(7:8) /0.207784955007898467600689403773245_8, 0.000000000000000000000000000000000_8/
    Data X(9:10) /-0.207784955007898467600689403773245_8, -0.405845151377397166906606412076961_8/
    Data X(11:12) /-0.586087235467691130294144838258730_8,-0.741531185599394439863864773280788_8/
    Data X(13:14) /-0.864864423359769072789712788640926_8, -0.949107912342758524526189684047851_8/
    Data X(15) /-0.991455371120812639206854697526329_8/
    IG=0._8
    IK=0._8
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,15
        Y=C1*X(Kstep)+C2
        Y=FuncS(Y)
        IK=IK+WK(Kstep)*Y
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Y*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine IntGaussKron15Real

Recursive Subroutine IntGaussKron15RealN(Llim, Hlim, Razb, LResult, Err, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Err, Perr
    Double Precision:: PResult, LResult, FuncN
    Integer:: Lstep, Razb
    External FuncN
    LResult=0._8
    Err=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussKron15Real(Alim, Blim, PResult, Perr, FuncN)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussKron15Real(Alim, Blim, PResult, Perr, FuncN)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine IntGaussKron15RealN

!=====================================================
!   Gauss-Kronrod method of 15 order integration for double precision array functions.
!=====================================================
Recursive Subroutine ArrIntGaussKron15Real(Alim, Blim, IK, Error, FuncS, Razm)
    Implicit None
    Double Precision:: WG(7), WK(15), X(15), Alim, Blim
    Double Precision:: Y, C1, C2
    Integer:: Kstep, Lstep, I, Razm
    Double Precision, Dimension(Razm):: Tmp, IG, IK, Error
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Data WG(1:2) /0.129484966168869693270611432679082_8, 0.279705391489276667901467771423780_8/
    Data WG(3:4) /0.381830050505118944950369775488975_8, 0.417959183673469387755102040816327_8/
    Data WG(5:6) /0.381830050505118944950369775488975_8, 0.279705391489276667901467771423780_8/
    Data WG(7) /0.129484966168869693270611432679082_8/
    Data WK(1:2) /0.022935322010529224963732008058970_8, 0.063092092629978553290700663189204_8/
    Data WK(3:4) /0.104790010322250183839876322541518_8, 0.140653259715525918745189590510238_8/
    Data WK(5:6) /0.169004726639267902826583426598550_8, 0.190350578064785409913256402421014_8/
    Data WK(7:8) /0.204432940075298892414161999234649_8, 0.209482141084727828012999174891714_8/
    Data WK(9:10) /0.204432940075298892414161999234649_8, 0.190350578064785409913256402421014_8/
    Data WK(11:12) /0.169004726639267902826583426598550_8, 0.140653259715525918745189590510238_8/
    Data WK(13:14) /0.104790010322250183839876322541518_8, 0.063092092629978553290700663189204_8/
    Data WK(15) /0.022935322010529224963732008058970_8/
    Data X(1:2) /0.991455371120812639206854697526329_8, 0.949107912342758524526189684047851_8/
    Data X(3:4) /0.864864423359769072789712788640926_8, 0.741531185599394439863864773280788_8/
    Data X(5:6) /0.586087235467691130294144838258730_8, 0.405845151377397166906606412076961_8/
    Data X(7:8) /0.207784955007898467600689403773245_8, 0.000000000000000000000000000000000_8/
    Data X(9:10) /-0.207784955007898467600689403773245_8, -0.405845151377397166906606412076961_8/
    Data X(11:12) /-0.586087235467691130294144838258730_8,-0.741531185599394439863864773280788_8/
    Data X(13:14) /-0.864864423359769072789712788640926_8, -0.949107912342758524526189684047851_8/
    Data X(15) /-0.991455371120812639206854697526329_8/
    IG=0._8
    IK=0._8
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,15
        Y=C1*X(Kstep)+C2
        Tmp=FuncS(Y, Razm)
        IK=IK+WK(Kstep)*Tmp
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Tmp*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine ArrIntGaussKron15Real

Recursive Subroutine ArrIntGaussKron15RealN(Llim, Hlim, Razb, LResult, Err, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Lstep, Razb, Razm
    Double Precision, Dimension(Razm):: PResult, LResult, Err, Perr
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Err=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussKron15Real(Alim, Blim, PResult, Perr, FuncN, Razm)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussKron15Real(Alim, Blim, PResult, Perr, FuncN, Razm)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine ArrIntGaussKron15RealN

!=====================================================
!   Gauss-Kronrod method of 15 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntGaussKron15Comp(Alim, Blim, IK, Error, FuncS)
    Implicit None
    Double Precision:: WG(7), WK(15), X(15), Alim, Blim, Error, Y
    Double Complex:: IG, IK, FuncS, C1, C2
    Integer:: Kstep, Lstep, I
    External FuncS
    Data WG(1:2) /0.129484966168869693270611432679082_8, 0.279705391489276667901467771423780_8/
    Data WG(3:4) /0.381830050505118944950369775488975_8, 0.417959183673469387755102040816327_8/
    Data WG(5:6) /0.381830050505118944950369775488975_8, 0.279705391489276667901467771423780_8/
    Data WG(7) /0.129484966168869693270611432679082_8/
    Data WK(1:2) /0.022935322010529224963732008058970_8, 0.063092092629978553290700663189204_8/
    Data WK(3:4) /0.104790010322250183839876322541518_8, 0.140653259715525918745189590510238_8/
    Data WK(5:6) /0.169004726639267902826583426598550_8, 0.190350578064785409913256402421014_8/
    Data WK(7:8) /0.204432940075298892414161999234649_8, 0.209482141084727828012999174891714_8/
    Data WK(9:10) /0.204432940075298892414161999234649_8, 0.190350578064785409913256402421014_8/
    Data WK(11:12) /0.169004726639267902826583426598550_8, 0.140653259715525918745189590510238_8/
    Data WK(13:14) /0.104790010322250183839876322541518_8, 0.063092092629978553290700663189204_8/
    Data WK(15) /0.022935322010529224963732008058970_8/
    Data X(1:2) /0.991455371120812639206854697526329_8, 0.949107912342758524526189684047851_8/
    Data X(3:4) /0.864864423359769072789712788640926_8, 0.741531185599394439863864773280788_8/
    Data X(5:6) /0.586087235467691130294144838258730_8, 0.405845151377397166906606412076961_8/
    Data X(7:8) /0.207784955007898467600689403773245_8, 0.000000000000000000000000000000000_8/
    Data X(9:10) /-0.207784955007898467600689403773245_8, -0.405845151377397166906606412076961_8/
    Data X(11:12) /-0.586087235467691130294144838258730_8,-0.741531185599394439863864773280788_8/
    Data X(13:14) /-0.864864423359769072789712788640926_8, -0.949107912342758524526189684047851_8/
    Data X(15) /-0.991455371120812639206854697526329_8/
    IG=0._8
    IK=0._8
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,15
        Y=C1*X(Kstep)+C2
        Y=FuncS(Y)
        IK=IK+WK(Kstep)*Y
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Y*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine IntGaussKron15Comp

Recursive Subroutine IntGaussKron15CompN(Llim, Hlim, Razb, LResult, Err, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Err, Perr
    Double Complex:: PResult, LResult, FuncN
    Integer:: Lstep, Razb
    External FuncN
    LResult=0._8
    Err=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussKron15Comp(Alim, Blim, PResult, Perr, FuncN)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussKron15Comp(Alim, Blim, PResult, Perr, FuncN)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine IntGaussKron15CompN

!=====================================================
!   Gauss-Kronrod method of 15 order integration for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntGaussKron15Comp(Alim, Blim, IK, Error, FuncS, Razm)
    Implicit None
    Double Precision:: WG(7), WK(15), X(15), Alim, Blim, Y, C1, C2
    Integer:: Kstep, Lstep, I, Razm
    Double Complex, Dimension(Razm):: Tmp, IG, IK
    Double Precision, Dimension(Razm):: Error
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Data WG(1:2) /0.129484966168869693270611432679082_8, 0.279705391489276667901467771423780_8/
    Data WG(3:4) /0.381830050505118944950369775488975_8, 0.417959183673469387755102040816327_8/
    Data WG(5:6) /0.381830050505118944950369775488975_8, 0.279705391489276667901467771423780_8/
    Data WG(7) /0.129484966168869693270611432679082_8/
    Data WK(1:2) /0.022935322010529224963732008058970_8, 0.063092092629978553290700663189204_8/
    Data WK(3:4) /0.104790010322250183839876322541518_8, 0.140653259715525918745189590510238_8/
    Data WK(5:6) /0.169004726639267902826583426598550_8, 0.190350578064785409913256402421014_8/
    Data WK(7:8) /0.204432940075298892414161999234649_8, 0.209482141084727828012999174891714_8/
    Data WK(9:10) /0.204432940075298892414161999234649_8, 0.190350578064785409913256402421014_8/
    Data WK(11:12) /0.169004726639267902826583426598550_8, 0.140653259715525918745189590510238_8/
    Data WK(13:14) /0.104790010322250183839876322541518_8, 0.063092092629978553290700663189204_8/
    Data WK(15) /0.022935322010529224963732008058970_8/
    Data X(1:2) /0.991455371120812639206854697526329_8, 0.949107912342758524526189684047851_8/
    Data X(3:4) /0.864864423359769072789712788640926_8, 0.741531185599394439863864773280788_8/
    Data X(5:6) /0.586087235467691130294144838258730_8, 0.405845151377397166906606412076961_8/
    Data X(7:8) /0.207784955007898467600689403773245_8, 0.000000000000000000000000000000000_8/
    Data X(9:10) /-0.207784955007898467600689403773245_8, -0.405845151377397166906606412076961_8/
    Data X(11:12) /-0.586087235467691130294144838258730_8,-0.741531185599394439863864773280788_8/
    Data X(13:14) /-0.864864423359769072789712788640926_8, -0.949107912342758524526189684047851_8/
    Data X(15) /-0.991455371120812639206854697526329_8/
    IG=(0._8,0._8)
    IK=(0._8,0._8)
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,15
        Y=C1*X(Kstep)+C2
        Tmp=FuncS(Y, Razm)
        IK=IK+WK(Kstep)*Tmp
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Tmp*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine ArrIntGaussKron15Comp

Recursive Subroutine ArrIntGaussKron15CompN(Llim, Hlim, Razb, LResult, Err, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Double Precision, Dimension(Razm):: Err, Perr
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=(0._8,0._8)
    Err=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussKron15Comp(Alim, Blim, PResult, Perr, FuncN, Razm)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussKron15Comp(Alim, Blim, PResult, Perr, FuncN, Razm)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine ArrIntGaussKron15CompN

!=====================================================
!   Gauss-Kronrod method of 41 order integration for double precision functions.
!=====================================================
Recursive Subroutine IntGaussKron41Real(Alim, Blim, IK, Error, FuncS)
    Implicit None
    Double Precision:: WG(20), WK(41), X(41), Alim, Blim, Error, Y, C1, C2
    Double Precision:: IG, IK, FuncS
    Integer:: Kstep, Lstep, I
    External FuncS
    Data WG(1:2) /0.017614007139152118311861962351853_8, 0.040601429800386941331039952274932_8/
    Data WG(3:4) /0.062672048334109063569506535187042_8, 0.083276741576704748724758143222046_8/
    Data WG(5:6) /0.101930119817240435036750135480350_8, 0.118194531961518417312377377711382_8/
    Data WG(7:8) /0.131688638449176626898494499748163_8, 0.142096109318382051329298325067165_8/
    Data WG(9:10) /0.149172986472603746787828737001969_8, 0.152753387130725850698084331955098_8/
    Data WG(11:12) /0.152753387130725850698084331955098_8, 0.149172986472603746787828737001969_8/
    Data WG(13:14) /0.142096109318382051329298325067165_8, 0.131688638449176626898494499748163_8/
    Data WG(15:16) /0.118194531961518417312377377711382_8, 0.101930119817240435036750135480350_8/
    Data WG(17:18) /0.083276741576704748724758143222046_8, 0.062672048334109063569506535187042_8/
    Data WG(19:20) /0.040601429800386941331039952274932_8, 0.017614007139152118311861962351853_8/
    Data WK(1:2) /0.003073583718520531501218293246031_8, 0.008600269855642942198661787950102_8/
    Data WK(3:4) /0.014626169256971252983787960308868_8, 0.020388373461266523598010231432755_8/
    Data WK(5:6) /0.025882133604951158834505067096153_8, 0.031287306777032798958543119323801_8/
    Data WK(7:8) /0.036600169758200798030557240707211_8, 0.041668873327973686263788305936895_8/
    Data WK(9:10) /0.046434821867497674720231880926108_8, 0.050944573923728691932707670050345_8/
    Data WK(11:12) /0.055195105348285994744832372419777_8, 0.059111400880639572374967220648594_8/
    Data WK(13:14) /0.062653237554781168025870122174255_8, 0.065834597133618422111563556969398_8/
    Data WK(15:16) /0.068648672928521619345623411885368_8, 0.071054423553444068305790361723210_8/
    Data WK(17:18) /0.073030690332786667495189417658913_8, 0.074582875400499188986581418362488_8/
    Data WK(19:20) /0.075704497684556674659542775376617_8, 0.076377867672080736705502835038061_8/
    Data WK(21) /0.076600711917999656445049901530102_8/
    Data WK(22:23) /0.076377867672080736705502835038061_8, 0.075704497684556674659542775376617_8/
    Data WK(24:25) /0.074582875400499188986581418362488_8, 0.073030690332786667495189417658913_8/
    Data WK(26:27) /0.071054423553444068305790361723210_8, 0.068648672928521619345623411885368_8/
    Data WK(28:29) /0.065834597133618422111563556969398_8, 0.062653237554781168025870122174255_8/
    Data WK(30:31) /0.059111400880639572374967220648594_8, 0.055195105348285994744832372419777_8/
    Data WK(32:33) /0.050944573923728691932707670050345_8, 0.046434821867497674720231880926108_8/
    Data WK(34:35) /0.041668873327973686263788305936895_8, 0.036600169758200798030557240707211_8/
    Data WK(36:37) /0.031287306777032798958543119323801_8, 0.025882133604951158834505067096153_8/
    Data WK(38:39) /0.020388373461266523598010231432755_8, 0.014626169256971252983787960308868_8/
    Data WK(40:41) /0.008600269855642942198661787950102_8, 0.003073583718520531501218293246031_8/
    Data X(1:2) /0.998859031588277663838315576545863_8, 0.993128599185094924786122388471320_8/
    Data X(3:4) /0.981507877450250259193342994720217_8, 0.963971927277913791267666131197277_8/
    Data X(5:6) /0.940822633831754753519982722212443_8, 0.912234428251325905867752441203298_8/
    Data X(7:8) /0.878276811252281976077442995113078_8, 0.839116971822218823394529061701521_8/
    Data X(9:10) /0.795041428837551198350638833272788_8, 0.746331906460150792614305070355642_8/
    Data X(11:12) /0.693237656334751384805490711845932_8, 0.636053680726515025452836696226286_8/
    Data X(13:14) /0.575140446819710315342946036586425_8, 0.510867001950827098004364050955251_8/
    Data X(15:16) /0.443593175238725103199992213492640_8, 0.373706088715419560672548177024927_8/
    Data X(17:18) /0.301627868114913004320555356858592_8, 0.227785851141645078080496195368575_8/
    Data X(19:20) /0.152605465240922675505220241022678_8, 0.076526521133497333754640409398838_8/
    Data X(21) /0._8/
    Data X(22:23) /-0.076526521133497333754640409398838_8, -0.152605465240922675505220241022678_8/
    Data X(24:25) /-0.227785851141645078080496195368575_8, -0.301627868114913004320555356858592_8/
    Data X(26:27) /-0.373706088715419560672548177024927_8, -0.443593175238725103199992213492640_8/
    Data X(28:29) /-0.510867001950827098004364050955251_8, -0.575140446819710315342946036586425_8/
    Data X(30:31) /-0.636053680726515025452836696226286_8, -0.693237656334751384805490711845932_8/
    Data X(32:33) /-0.746331906460150792614305070355642_8, -0.795041428837551198350638833272788_8/
    Data X(34:35) /-0.839116971822218823394529061701521_8, -0.878276811252281976077442995113078_8/
    Data X(36:37) /-0.912234428251325905867752441203298_8, -0.940822633831754753519982722212443_8/
    Data X(38:39) /-0.963971927277913791267666131197277_8, -0.981507877450250259193342994720217_8/
    Data X(40:41) /-0.993128599185094924786122388471320_8, -0.998859031588277663838315576545863_8/
    IG=0._8
    IK=0._8
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,41
        Y=C1*X(Kstep)+C2
        Y=FuncS(Y)
        IK=IK+WK(Kstep)*Y
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Y*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine IntGaussKron41Real

Recursive Subroutine IntGaussKron41RealN(Llim, Hlim, Razb, LResult, Err, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Err, Perr
    Double Precision:: PResult, LResult, FuncN
    Integer:: Lstep, Razb
    External FuncN
    LResult=0._8
    Err=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussKron41Real(Alim, Blim, PResult, Perr, FuncN)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussKron41Real(Alim, Blim, PResult, Perr, FuncN)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine IntGaussKron41RealN

!=====================================================
!   Gauss-Kronrod method of 41 order integration for double precision array functions.
!=====================================================
Recursive Subroutine ArrIntGaussKron41Real(Alim, Blim, IK, Error, FuncS, Razm)
    Implicit None
    Double Precision:: WG(20), WK(41), X(41), Alim, Blim
    Double Precision:: Y, C1, C2
    Integer:: Kstep, Lstep, I, Razm
    Double Precision, Dimension(Razm):: Tmp, IG, IK, Error
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Data WG(1:2) /0.017614007139152118311861962351853_8, 0.040601429800386941331039952274932_8/
    Data WG(3:4) /0.062672048334109063569506535187042_8, 0.083276741576704748724758143222046_8/
    Data WG(5:6) /0.101930119817240435036750135480350_8, 0.118194531961518417312377377711382_8/
    Data WG(7:8) /0.131688638449176626898494499748163_8, 0.142096109318382051329298325067165_8/
    Data WG(9:10) /0.149172986472603746787828737001969_8, 0.152753387130725850698084331955098_8/
    Data WG(11:12) /0.152753387130725850698084331955098_8, 0.149172986472603746787828737001969_8/
    Data WG(13:14) /0.142096109318382051329298325067165_8, 0.131688638449176626898494499748163_8/
    Data WG(15:16) /0.118194531961518417312377377711382_8, 0.101930119817240435036750135480350_8/
    Data WG(17:18) /0.083276741576704748724758143222046_8, 0.062672048334109063569506535187042_8/
    Data WG(19:20) /0.040601429800386941331039952274932_8, 0.017614007139152118311861962351853_8/
    Data WK(1:2) /0.003073583718520531501218293246031_8, 0.008600269855642942198661787950102_8/
    Data WK(3:4) /0.014626169256971252983787960308868_8, 0.020388373461266523598010231432755_8/
    Data WK(5:6) /0.025882133604951158834505067096153_8, 0.031287306777032798958543119323801_8/
    Data WK(7:8) /0.036600169758200798030557240707211_8, 0.041668873327973686263788305936895_8/
    Data WK(9:10) /0.046434821867497674720231880926108_8, 0.050944573923728691932707670050345_8/
    Data WK(11:12) /0.055195105348285994744832372419777_8, 0.059111400880639572374967220648594_8/
    Data WK(13:14) /0.062653237554781168025870122174255_8, 0.065834597133618422111563556969398_8/
    Data WK(15:16) /0.068648672928521619345623411885368_8, 0.071054423553444068305790361723210_8/
    Data WK(17:18) /0.073030690332786667495189417658913_8, 0.074582875400499188986581418362488_8/
    Data WK(19:20) /0.075704497684556674659542775376617_8, 0.076377867672080736705502835038061_8/
    Data WK(21) /0.076600711917999656445049901530102_8/
    Data WK(22:23) /0.076377867672080736705502835038061_8, 0.075704497684556674659542775376617_8/
    Data WK(24:25) /0.074582875400499188986581418362488_8, 0.073030690332786667495189417658913_8/
    Data WK(26:27) /0.071054423553444068305790361723210_8, 0.068648672928521619345623411885368_8/
    Data WK(28:29) /0.065834597133618422111563556969398_8, 0.062653237554781168025870122174255_8/
    Data WK(30:31) /0.059111400880639572374967220648594_8, 0.055195105348285994744832372419777_8/
    Data WK(32:33) /0.050944573923728691932707670050345_8, 0.046434821867497674720231880926108_8/
    Data WK(34:35) /0.041668873327973686263788305936895_8, 0.036600169758200798030557240707211_8/
    Data WK(36:37) /0.031287306777032798958543119323801_8, 0.025882133604951158834505067096153_8/
    Data WK(38:39) /0.020388373461266523598010231432755_8, 0.014626169256971252983787960308868_8/
    Data WK(40:41) /0.008600269855642942198661787950102_8, 0.003073583718520531501218293246031_8/
    Data X(1:2) /0.998859031588277663838315576545863_8, 0.993128599185094924786122388471320_8/
    Data X(3:4) /0.981507877450250259193342994720217_8, 0.963971927277913791267666131197277_8/
    Data X(5:6) /0.940822633831754753519982722212443_8, 0.912234428251325905867752441203298_8/
    Data X(7:8) /0.878276811252281976077442995113078_8, 0.839116971822218823394529061701521_8/
    Data X(9:10) /0.795041428837551198350638833272788_8, 0.746331906460150792614305070355642_8/
    Data X(11:12) /0.693237656334751384805490711845932_8, 0.636053680726515025452836696226286_8/
    Data X(13:14) /0.575140446819710315342946036586425_8, 0.510867001950827098004364050955251_8/
    Data X(15:16) /0.443593175238725103199992213492640_8, 0.373706088715419560672548177024927_8/
    Data X(17:18) /0.301627868114913004320555356858592_8, 0.227785851141645078080496195368575_8/
    Data X(19:20) /0.152605465240922675505220241022678_8, 0.076526521133497333754640409398838_8/
    Data X(21) /0._8/
    Data X(22:23) /-0.076526521133497333754640409398838_8, -0.152605465240922675505220241022678_8/
    Data X(24:25) /-0.227785851141645078080496195368575_8, -0.301627868114913004320555356858592_8/
    Data X(26:27) /-0.373706088715419560672548177024927_8, -0.443593175238725103199992213492640_8/
    Data X(28:29) /-0.510867001950827098004364050955251_8, -0.575140446819710315342946036586425_8/
    Data X(30:31) /-0.636053680726515025452836696226286_8, -0.693237656334751384805490711845932_8/
    Data X(32:33) /-0.746331906460150792614305070355642_8, -0.795041428837551198350638833272788_8/
    Data X(34:35) /-0.839116971822218823394529061701521_8, -0.878276811252281976077442995113078_8/
    Data X(36:37) /-0.912234428251325905867752441203298_8, -0.940822633831754753519982722212443_8/
    Data X(38:39) /-0.963971927277913791267666131197277_8, -0.981507877450250259193342994720217_8/
    Data X(40:41) /-0.993128599185094924786122388471320_8, -0.998859031588277663838315576545863_8/
    IG=0._8
    IK=0._8
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,41
        Y=C1*X(Kstep)+C2
        Tmp=FuncS(Y, Razm)
        IK=IK+WK(Kstep)*Tmp
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Tmp*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine ArrIntGaussKron41Real

Recursive Subroutine ArrIntGaussKron41RealN(Llim, Hlim, Razb, LResult, Err, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Lstep, Razb, Razm
    Double Precision, Dimension(Razm):: PResult, LResult, Err, Perr
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Err=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussKron41Real(Alim, Blim, PResult, Perr, FuncN, Razm)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussKron41Real(Alim, Blim, PResult, Perr, FuncN, Razm)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine ArrIntGaussKron41RealN

!=====================================================
!   Gauss-Kronrod method of 41 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntGaussKron41Comp(Alim, Blim, IK, Error, FuncS)
    Implicit None
    Double Precision:: WG(20), WK(41), X(41), Alim, Blim, Error, Y, C1, C2
    Double Complex:: IG, IK, FuncS
    Integer:: Kstep, Lstep, I
    External FuncS
    Data WG(1:2) /0.017614007139152118311861962351853_8, 0.040601429800386941331039952274932_8/
    Data WG(3:4) /0.062672048334109063569506535187042_8, 0.083276741576704748724758143222046_8/
    Data WG(5:6) /0.101930119817240435036750135480350_8, 0.118194531961518417312377377711382_8/
    Data WG(7:8) /0.131688638449176626898494499748163_8, 0.142096109318382051329298325067165_8/
    Data WG(9:10) /0.149172986472603746787828737001969_8, 0.152753387130725850698084331955098_8/
    Data WG(11:12) /0.152753387130725850698084331955098_8, 0.149172986472603746787828737001969_8/
    Data WG(13:14) /0.142096109318382051329298325067165_8, 0.131688638449176626898494499748163_8/
    Data WG(15:16) /0.118194531961518417312377377711382_8, 0.101930119817240435036750135480350_8/
    Data WG(17:18) /0.083276741576704748724758143222046_8, 0.062672048334109063569506535187042_8/
    Data WG(19:20) /0.040601429800386941331039952274932_8, 0.017614007139152118311861962351853_8/
    Data WK(1:2) /0.003073583718520531501218293246031_8, 0.008600269855642942198661787950102_8/
    Data WK(3:4) /0.014626169256971252983787960308868_8, 0.020388373461266523598010231432755_8/
    Data WK(5:6) /0.025882133604951158834505067096153_8, 0.031287306777032798958543119323801_8/
    Data WK(7:8) /0.036600169758200798030557240707211_8, 0.041668873327973686263788305936895_8/
    Data WK(9:10) /0.046434821867497674720231880926108_8, 0.050944573923728691932707670050345_8/
    Data WK(11:12) /0.055195105348285994744832372419777_8, 0.059111400880639572374967220648594_8/
    Data WK(13:14) /0.062653237554781168025870122174255_8, 0.065834597133618422111563556969398_8/
    Data WK(15:16) /0.068648672928521619345623411885368_8, 0.071054423553444068305790361723210_8/
    Data WK(17:18) /0.073030690332786667495189417658913_8, 0.074582875400499188986581418362488_8/
    Data WK(19:20) /0.075704497684556674659542775376617_8, 0.076377867672080736705502835038061_8/
    Data WK(21) /0.076600711917999656445049901530102_8/
    Data WK(22:23) /0.076377867672080736705502835038061_8, 0.075704497684556674659542775376617_8/
    Data WK(24:25) /0.074582875400499188986581418362488_8, 0.073030690332786667495189417658913_8/
    Data WK(26:27) /0.071054423553444068305790361723210_8, 0.068648672928521619345623411885368_8/
    Data WK(28:29) /0.065834597133618422111563556969398_8, 0.062653237554781168025870122174255_8/
    Data WK(30:31) /0.059111400880639572374967220648594_8, 0.055195105348285994744832372419777_8/
    Data WK(32:33) /0.050944573923728691932707670050345_8, 0.046434821867497674720231880926108_8/
    Data WK(34:35) /0.041668873327973686263788305936895_8, 0.036600169758200798030557240707211_8/
    Data WK(36:37) /0.031287306777032798958543119323801_8, 0.025882133604951158834505067096153_8/
    Data WK(38:39) /0.020388373461266523598010231432755_8, 0.014626169256971252983787960308868_8/
    Data WK(40:41) /0.008600269855642942198661787950102_8, 0.003073583718520531501218293246031_8/
    Data X(1:2) /0.998859031588277663838315576545863_8, 0.993128599185094924786122388471320_8/
    Data X(3:4) /0.981507877450250259193342994720217_8, 0.963971927277913791267666131197277_8/
    Data X(5:6) /0.940822633831754753519982722212443_8, 0.912234428251325905867752441203298_8/
    Data X(7:8) /0.878276811252281976077442995113078_8, 0.839116971822218823394529061701521_8/
    Data X(9:10) /0.795041428837551198350638833272788_8, 0.746331906460150792614305070355642_8/
    Data X(11:12) /0.693237656334751384805490711845932_8, 0.636053680726515025452836696226286_8/
    Data X(13:14) /0.575140446819710315342946036586425_8, 0.510867001950827098004364050955251_8/
    Data X(15:16) /0.443593175238725103199992213492640_8, 0.373706088715419560672548177024927_8/
    Data X(17:18) /0.301627868114913004320555356858592_8, 0.227785851141645078080496195368575_8/
    Data X(19:20) /0.152605465240922675505220241022678_8, 0.076526521133497333754640409398838_8/
    Data X(21) /0._8/
    Data X(22:23) /-0.076526521133497333754640409398838_8, -0.152605465240922675505220241022678_8/
    Data X(24:25) /-0.227785851141645078080496195368575_8, -0.301627868114913004320555356858592_8/
    Data X(26:27) /-0.373706088715419560672548177024927_8, -0.443593175238725103199992213492640_8/
    Data X(28:29) /-0.510867001950827098004364050955251_8, -0.575140446819710315342946036586425_8/
    Data X(30:31) /-0.636053680726515025452836696226286_8, -0.693237656334751384805490711845932_8/
    Data X(32:33) /-0.746331906460150792614305070355642_8, -0.795041428837551198350638833272788_8/
    Data X(34:35) /-0.839116971822218823394529061701521_8, -0.878276811252281976077442995113078_8/
    Data X(36:37) /-0.912234428251325905867752441203298_8, -0.940822633831754753519982722212443_8/
    Data X(38:39) /-0.963971927277913791267666131197277_8, -0.981507877450250259193342994720217_8/
    Data X(40:41) /-0.993128599185094924786122388471320_8, -0.998859031588277663838315576545863_8/
    IG=(0._8,0._8)
    IK=(0._8,0._8)
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,41
        Y=C1*X(Kstep)+C2
        Y=FuncS(Y)
        IK=IK+WK(Kstep)*Y
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Y*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine IntGaussKron41Comp

Recursive Subroutine IntGaussKron41CompN(Llim, Hlim, Razb, LResult, Err, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Err, Perr
    Double Complex:: PResult, LResult, FuncN
    Integer:: Lstep, Razb
    External FuncN
    LResult=(0._8,0._8)
    Err=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussKron41Comp(Alim, Blim, PResult, Perr, FuncN)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussKron41Comp(Alim, Blim, PResult, Perr, FuncN)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine IntGaussKron41CompN

!=====================================================
!   Gauss-Kronrod method of 41 order integration for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntGaussKron41Comp(Alim, Blim, IK, Error, FuncS, Razm)
    Implicit None
    Double Precision:: WG(20), WK(41), X(41), Alim, Blim, Y, C1, C2
    Integer:: Kstep, Lstep, I, Razm
    Double Complex, Dimension(Razm):: Tmp, IG, IK
    Double Precision, Dimension(Razm):: Error
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Data WG(1:2) /0.017614007139152118311861962351853_8, 0.040601429800386941331039952274932_8/
    Data WG(3:4) /0.062672048334109063569506535187042_8, 0.083276741576704748724758143222046_8/
    Data WG(5:6) /0.101930119817240435036750135480350_8, 0.118194531961518417312377377711382_8/
    Data WG(7:8) /0.131688638449176626898494499748163_8, 0.142096109318382051329298325067165_8/
    Data WG(9:10) /0.149172986472603746787828737001969_8, 0.152753387130725850698084331955098_8/
    Data WG(11:12) /0.152753387130725850698084331955098_8, 0.149172986472603746787828737001969_8/
    Data WG(13:14) /0.142096109318382051329298325067165_8, 0.131688638449176626898494499748163_8/
    Data WG(15:16) /0.118194531961518417312377377711382_8, 0.101930119817240435036750135480350_8/
    Data WG(17:18) /0.083276741576704748724758143222046_8, 0.062672048334109063569506535187042_8/
    Data WG(19:20) /0.040601429800386941331039952274932_8, 0.017614007139152118311861962351853_8/
    Data WK(1:2) /0.003073583718520531501218293246031_8, 0.008600269855642942198661787950102_8/
    Data WK(3:4) /0.014626169256971252983787960308868_8, 0.020388373461266523598010231432755_8/
    Data WK(5:6) /0.025882133604951158834505067096153_8, 0.031287306777032798958543119323801_8/
    Data WK(7:8) /0.036600169758200798030557240707211_8, 0.041668873327973686263788305936895_8/
    Data WK(9:10) /0.046434821867497674720231880926108_8, 0.050944573923728691932707670050345_8/
    Data WK(11:12) /0.055195105348285994744832372419777_8, 0.059111400880639572374967220648594_8/
    Data WK(13:14) /0.062653237554781168025870122174255_8, 0.065834597133618422111563556969398_8/
    Data WK(15:16) /0.068648672928521619345623411885368_8, 0.071054423553444068305790361723210_8/
    Data WK(17:18) /0.073030690332786667495189417658913_8, 0.074582875400499188986581418362488_8/
    Data WK(19:20) /0.075704497684556674659542775376617_8, 0.076377867672080736705502835038061_8/
    Data WK(21) /0.076600711917999656445049901530102_8/
    Data WK(22:23) /0.076377867672080736705502835038061_8, 0.075704497684556674659542775376617_8/
    Data WK(24:25) /0.074582875400499188986581418362488_8, 0.073030690332786667495189417658913_8/
    Data WK(26:27) /0.071054423553444068305790361723210_8, 0.068648672928521619345623411885368_8/
    Data WK(28:29) /0.065834597133618422111563556969398_8, 0.062653237554781168025870122174255_8/
    Data WK(30:31) /0.059111400880639572374967220648594_8, 0.055195105348285994744832372419777_8/
    Data WK(32:33) /0.050944573923728691932707670050345_8, 0.046434821867497674720231880926108_8/
    Data WK(34:35) /0.041668873327973686263788305936895_8, 0.036600169758200798030557240707211_8/
    Data WK(36:37) /0.031287306777032798958543119323801_8, 0.025882133604951158834505067096153_8/
    Data WK(38:39) /0.020388373461266523598010231432755_8, 0.014626169256971252983787960308868_8/
    Data WK(40:41) /0.008600269855642942198661787950102_8, 0.003073583718520531501218293246031_8/
    Data X(1:2) /0.998859031588277663838315576545863_8, 0.993128599185094924786122388471320_8/
    Data X(3:4) /0.981507877450250259193342994720217_8, 0.963971927277913791267666131197277_8/
    Data X(5:6) /0.940822633831754753519982722212443_8, 0.912234428251325905867752441203298_8/
    Data X(7:8) /0.878276811252281976077442995113078_8, 0.839116971822218823394529061701521_8/
    Data X(9:10) /0.795041428837551198350638833272788_8, 0.746331906460150792614305070355642_8/
    Data X(11:12) /0.693237656334751384805490711845932_8, 0.636053680726515025452836696226286_8/
    Data X(13:14) /0.575140446819710315342946036586425_8, 0.510867001950827098004364050955251_8/
    Data X(15:16) /0.443593175238725103199992213492640_8, 0.373706088715419560672548177024927_8/
    Data X(17:18) /0.301627868114913004320555356858592_8, 0.227785851141645078080496195368575_8/
    Data X(19:20) /0.152605465240922675505220241022678_8, 0.076526521133497333754640409398838_8/
    Data X(21) /0._8/
    Data X(22:23) /-0.076526521133497333754640409398838_8, -0.152605465240922675505220241022678_8/
    Data X(24:25) /-0.227785851141645078080496195368575_8, -0.301627868114913004320555356858592_8/
    Data X(26:27) /-0.373706088715419560672548177024927_8, -0.443593175238725103199992213492640_8/
    Data X(28:29) /-0.510867001950827098004364050955251_8, -0.575140446819710315342946036586425_8/
    Data X(30:31) /-0.636053680726515025452836696226286_8, -0.693237656334751384805490711845932_8/
    Data X(32:33) /-0.746331906460150792614305070355642_8, -0.795041428837551198350638833272788_8/
    Data X(34:35) /-0.839116971822218823394529061701521_8, -0.878276811252281976077442995113078_8/
    Data X(36:37) /-0.912234428251325905867752441203298_8, -0.940822633831754753519982722212443_8/
    Data X(38:39) /-0.963971927277913791267666131197277_8, -0.981507877450250259193342994720217_8/
    Data X(40:41) /-0.993128599185094924786122388471320_8, -0.998859031588277663838315576545863_8/
    IG=(0._8,0._8)
    IK=(0._8,0._8)
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,41
        Y=C1*X(Kstep)+C2
        Tmp=FuncS(Y, Razm)
        IK=IK+WK(Kstep)*Tmp
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Tmp*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine ArrIntGaussKron41Comp

Recursive Subroutine ArrIntGaussKron41CompN(Llim, Hlim, Razb, LResult, Err, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Double Precision, Dimension(Razm):: Err, Perr
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=(0._8,0._8)
    Err=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussKron41Comp(Alim, Blim, PResult, Perr, FuncN, Razm)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussKron41Comp(Alim, Blim, PResult, Perr, FuncN, Razm)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine ArrIntGaussKron41CompN

!=====================================================
!   Gauss-Kronrod method of 61 order integration for double precision functions.
!=====================================================
Recursive Subroutine IntGaussKron61Real(Alim, Blim, IK, Error, FuncS)
    Implicit None
    Double Precision:: WG(30), WK(61), X(61), Alim, Blim, Error
    Double Precision:: Y, IG, IK, FuncS, C1, C2
    Integer:: Kstep, Lstep, I
    External FuncS
    Data WG(1:2) /0.007968192496166605615465883474674_8, 0.018466468311090959142302131912047_8/
    Data WG(3:4) /0.028784707883323369349719179611292_8, 0.038799192569627049596801936446348_8/
    Data WG(5:6) /0.048402672830594052902938140422808_8, 0.057493156217619066481721689402056_8/
    Data WG(7:8) /0.065974229882180495128128515115962_8, 0.073755974737705206268243850022191_8/
    Data WG(9:10) /0.080755895229420215354694938460530_8, 0.086899787201082979802387530715126_8/
    Data WG(11:12) /0.092122522237786128717632707087619_8, 0.096368737174644259639468626351810_8/
    Data WG(13:14) /0.099593420586795267062780282103569_8, 0.101762389748405504596428952168554_8/
    Data WG(15:16) /0.102852652893558840341285636705415_8, 0.102852652893558840341285636705415_8/
    Data WG(17:18) /0.101762389748405504596428952168554_8, 0.099593420586795267062780282103569_8/
    Data WG(19:20) /0.096368737174644259639468626351810_8, 0.092122522237786128717632707087619_8/
    Data WG(21:22) /0.086899787201082979802387530715126_8, 0.080755895229420215354694938460530_8/
    Data WG(23:24) /0.073755974737705206268243850022191_8, 0.065974229882180495128128515115962_8/
    Data WG(25:26) /0.057493156217619066481721689402056_8, 0.048402672830594052902938140422808_8/
    Data WG(27:28) /0.038799192569627049596801936446348_8, 0.028784707883323369349719179611292_8/
    Data WG(29:30) /0.018466468311090959142302131912047_8, 0.007968192496166605615465883474674_8/   
    Data WK(1:2) /0.001389013698677007624551591226760_8, 0.003890461127099884051267201844516_8/
    Data WK(3:4) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(5:6) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(7:8) /0.016920889189053272627572289420322_8, 0.019414141193942381173408951050128_8/
    Data WK(9:10) /0.021828035821609192297167485738339_8, 0.024191162078080601365686370725232_8/
    Data WK(11:12) /0.026509954882333101610601709335075_8, 0.028754048765041292843978785354334_8/
    Data WK(13:14) /0.030907257562387762472884252943092_8, 0.032981447057483726031814191016854_8/
    Data WK(15:16) /0.034979338028060024137499670731468_8, 0.036882364651821229223911065617136_8/
    Data WK(17:18) /0.038678945624727592950348651532281_8, 0.040374538951535959111995279752468_8/
    Data WK(19:20) /0.041969810215164246147147541285970_8, 0.043452539701356069316831728117073_8/
    Data WK(21:22) /0.044814800133162663192355551616723_8, 0.046059238271006988116271735559374_8/
    Data WK(23:24) /0.047185546569299153945261478181099_8, 0.048185861757087129140779492298305_8/
    Data WK(25:26) /0.049055434555029778887528165367238_8, 0.049795683427074206357811569379942_8/
    Data WK(27:28) /0.050405921402782346840893085653585_8, 0.050881795898749606492297473049805_8/
    Data WK(29:30) /0.051221547849258772170656282604944_8, 0.051426128537459025933862879215781_8/
        Data WK(31) /0.051494729429451567558340433647099_8/
    Data WK(32:33) /0.051426128537459025933862879215781_8, 0.051221547849258772170656282604944_8/
    Data WK(34:35) /0.050881795898749606492297473049805_8, 0.050405921402782346840893085653585_8/
    Data WK(36:37) /0.049795683427074206357811569379942_8, 0.049055434555029778887528165367238_8/
    Data WK(38:39) /0.048185861757087129140779492298305_8, 0.047185546569299153945261478181099_8/
    Data WK(40:41) /0.046059238271006988116271735559374_8, 0.044814800133162663192355551616723_8/
    Data WK(42:43) /0.043452539701356069316831728117073_8, 0.041969810215164246147147541285970_8/
    Data WK(44:45) /0.040374538951535959111995279752468_8, 0.038678945624727592950348651532281_8/
    Data WK(46:47) /0.036882364651821229223911065617136_8, 0.034979338028060024137499670731468_8/
    Data WK(48:49) /0.032981447057483726031814191016854_8, 0.030907257562387762472884252943092_8/
    Data WK(50:51) /0.028754048765041292843978785354334_8, 0.026509954882333101610601709335075_8/
    Data WK(52:53) /0.024191162078080601365686370725232_8, 0.021828035821609192297167485738339_8/
    Data WK(54:55) /0.019414141193942381173408951050128_8, 0.016920889189053272627572289420322_8/
    Data WK(56:57) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(58:59) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(60:61) /0.003890461127099884051267201844516_8, 0.001389013698677007624551591226760_8/   
    Data X(1:2) /0.999484410050490637571325895705811_8, 0.996893484074649540271630050918695_8/
    Data X(3:4) /0.991630996870404594858628366109486_8, 0.983668123279747209970032581605663_8/
    Data X(5:6) /0.973116322501126268374693868423707_8, 0.960021864968307512216871025581798_8/
    Data X(7:8) /0.944374444748559979415831324037439_8, 0.926200047429274325879324277080474_8/
    Data X(9:10) /0.905573307699907798546522558925958_8, 0.882560535792052681543116462530226_8/
    Data X(11:12) /0.857205233546061098958658510658944_8, 0.829565762382768397442898119732502_8/
    Data X(13:14) /0.799727835821839083013668942322683_8, 0.767777432104826194917977340974503_8/
    Data X(15:16) /0.733790062453226804726171131369528_8, 0.697850494793315796932292388026640_8/
    Data X(17:18) /0.660061064126626961370053668149271_8, 0.620526182989242861140477556431189_8/
    Data X(19:20) /0.579345235826361691756024932172540_8, 0.536624148142019899264169793311073_8/
    Data X(21:22) /0.492480467861778574993693061207709_8, 0.447033769538089176780609900322854_8/
    Data X(23:24) /0.400401254830394392535476211542661_8, 0.352704725530878113471037207089374_8/
    Data X(25:26) /0.304073202273625077372677107199257_8, 0.254636926167889846439805129817805_8/
    Data X(27:28) /0.204525116682309891438957671002025_8, 0.153869913608583546963794672743256_8/
    Data X(29:30) /0.102806937966737030147096751318001_8, 0.051471842555317695833025213166723_8/
    Data X(31) /0._8/
    Data X(32:33) /-0.051471842555317695833025213166723_8, -0.102806937966737030147096751318001_8/
    Data X(34:35) /-0.153869913608583546963794672743256_8, -0.204525116682309891438957671002025_8/
    Data X(36:37) /-0.254636926167889846439805129817805_8, -0.304073202273625077372677107199257_8/
    Data X(38:39) /-0.352704725530878113471037207089374_8, -0.400401254830394392535476211542661_8/
    Data X(40:41) /-0.447033769538089176780609900322854_8, -0.492480467861778574993693061207709_8/
    Data X(42:43) /-0.536624148142019899264169793311073_8, -0.579345235826361691756024932172540_8/
    Data X(44:45) /-0.620526182989242861140477556431189_8, -0.660061064126626961370053668149271_8/
    Data X(46:47) /-0.697850494793315796932292388026640_8, -0.733790062453226804726171131369528_8/
    Data X(48:49) /-0.767777432104826194917977340974503_8, -0.799727835821839083013668942322683_8/
    Data X(50:51) /-0.829565762382768397442898119732502_8, -0.857205233546061098958658510658944_8/
    Data X(52:53) /-0.882560535792052681543116462530226_8, -0.905573307699907798546522558925958_8/
    Data X(54:55) /-0.926200047429274325879324277080474_8, -0.944374444748559979415831324037439_8/
    Data X(56:57) /-0.960021864968307512216871025581798_8, -0.973116322501126268374693868423707_8/
    Data X(58:59) /-0.983668123279747209970032581605663_8, -0.991630996870404594858628366109486_8/
    Data X(60:61) /-0.996893484074649540271630050918695_8, -0.999484410050490637571325895705811_8/
    IG=0._8
    IK=0._8
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,61
        Y=C1*X(Kstep)+C2
        Y=FuncS(Y)
        IK=IK+WK(Kstep)*Y
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Y*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine IntGaussKron61Real

Recursive Subroutine IntGaussKron61RealN(Llim, Hlim, Razb, LResult, Err, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Err, Perr
    Double Precision:: PResult, LResult, FuncN
    Integer:: Lstep, Razb
    External FuncN
    LResult=0._8
    Err=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussKron61Real(Alim, Blim, PResult, Perr, FuncN)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussKron61Real(Alim, Blim, PResult, Perr, FuncN)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine IntGaussKron61RealN

!=====================================================
!   Gauss-Kronrod method of 61 order integration for double precision array functions.
!=====================================================
Recursive Subroutine ArrIntGaussKron61Real(Alim, Blim, IK, Error, FuncS, Razm)
    Implicit None
    Double Precision:: WG(30), WK(61), X(61), Alim, Blim
    Double Precision:: Y, C1, C2
    Integer:: Kstep, Lstep, I, Razm
    Double Precision, Dimension(Razm):: Tmp, IG, IK, Error
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Data WG(1:2) /0.007968192496166605615465883474674_8, 0.018466468311090959142302131912047_8/
    Data WG(3:4) /0.028784707883323369349719179611292_8, 0.038799192569627049596801936446348_8/
    Data WG(5:6) /0.048402672830594052902938140422808_8, 0.057493156217619066481721689402056_8/
    Data WG(7:8) /0.065974229882180495128128515115962_8, 0.073755974737705206268243850022191_8/
    Data WG(9:10) /0.080755895229420215354694938460530_8, 0.086899787201082979802387530715126_8/
    Data WG(11:12) /0.092122522237786128717632707087619_8, 0.096368737174644259639468626351810_8/
    Data WG(13:14) /0.099593420586795267062780282103569_8, 0.101762389748405504596428952168554_8/
    Data WG(15:16) /0.102852652893558840341285636705415_8, 0.102852652893558840341285636705415_8/
    Data WG(17:18) /0.101762389748405504596428952168554_8, 0.099593420586795267062780282103569_8/
    Data WG(19:20) /0.096368737174644259639468626351810_8, 0.092122522237786128717632707087619_8/
    Data WG(21:22) /0.086899787201082979802387530715126_8, 0.080755895229420215354694938460530_8/
    Data WG(23:24) /0.073755974737705206268243850022191_8, 0.065974229882180495128128515115962_8/
    Data WG(25:26) /0.057493156217619066481721689402056_8, 0.048402672830594052902938140422808_8/
    Data WG(27:28) /0.038799192569627049596801936446348_8, 0.028784707883323369349719179611292_8/
    Data WG(29:30) /0.018466468311090959142302131912047_8, 0.007968192496166605615465883474674_8/   
    Data WK(1:2) /0.001389013698677007624551591226760_8, 0.003890461127099884051267201844516_8/
    Data WK(3:4) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(5:6) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(7:8) /0.016920889189053272627572289420322_8, 0.019414141193942381173408951050128_8/
    Data WK(9:10) /0.021828035821609192297167485738339_8, 0.024191162078080601365686370725232_8/
    Data WK(11:12) /0.026509954882333101610601709335075_8, 0.028754048765041292843978785354334_8/
    Data WK(13:14) /0.030907257562387762472884252943092_8, 0.032981447057483726031814191016854_8/
    Data WK(15:16) /0.034979338028060024137499670731468_8, 0.036882364651821229223911065617136_8/
    Data WK(17:18) /0.038678945624727592950348651532281_8, 0.040374538951535959111995279752468_8/
    Data WK(19:20) /0.041969810215164246147147541285970_8, 0.043452539701356069316831728117073_8/
    Data WK(21:22) /0.044814800133162663192355551616723_8, 0.046059238271006988116271735559374_8/
    Data WK(23:24) /0.047185546569299153945261478181099_8, 0.048185861757087129140779492298305_8/
    Data WK(25:26) /0.049055434555029778887528165367238_8, 0.049795683427074206357811569379942_8/
    Data WK(27:28) /0.050405921402782346840893085653585_8, 0.050881795898749606492297473049805_8/
    Data WK(29:30) /0.051221547849258772170656282604944_8, 0.051426128537459025933862879215781_8/
        Data WK(31) /0.051494729429451567558340433647099_8/
    Data WK(32:33) /0.051426128537459025933862879215781_8, 0.051221547849258772170656282604944_8/
    Data WK(34:35) /0.050881795898749606492297473049805_8, 0.050405921402782346840893085653585_8/
    Data WK(36:37) /0.049795683427074206357811569379942_8, 0.049055434555029778887528165367238_8/
    Data WK(38:39) /0.048185861757087129140779492298305_8, 0.047185546569299153945261478181099_8/
    Data WK(40:41) /0.046059238271006988116271735559374_8, 0.044814800133162663192355551616723_8/
    Data WK(42:43) /0.043452539701356069316831728117073_8, 0.041969810215164246147147541285970_8/
    Data WK(44:45) /0.040374538951535959111995279752468_8, 0.038678945624727592950348651532281_8/
    Data WK(46:47) /0.036882364651821229223911065617136_8, 0.034979338028060024137499670731468_8/
    Data WK(48:49) /0.032981447057483726031814191016854_8, 0.030907257562387762472884252943092_8/
    Data WK(50:51) /0.028754048765041292843978785354334_8, 0.026509954882333101610601709335075_8/
    Data WK(52:53) /0.024191162078080601365686370725232_8, 0.021828035821609192297167485738339_8/
    Data WK(54:55) /0.019414141193942381173408951050128_8, 0.016920889189053272627572289420322_8/
    Data WK(56:57) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(58:59) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(60:61) /0.003890461127099884051267201844516_8, 0.001389013698677007624551591226760_8/   
    Data X(1:2) /0.999484410050490637571325895705811_8, 0.996893484074649540271630050918695_8/
    Data X(3:4) /0.991630996870404594858628366109486_8, 0.983668123279747209970032581605663_8/
    Data X(5:6) /0.973116322501126268374693868423707_8, 0.960021864968307512216871025581798_8/
    Data X(7:8) /0.944374444748559979415831324037439_8, 0.926200047429274325879324277080474_8/
    Data X(9:10) /0.905573307699907798546522558925958_8, 0.882560535792052681543116462530226_8/
    Data X(11:12) /0.857205233546061098958658510658944_8, 0.829565762382768397442898119732502_8/
    Data X(13:14) /0.799727835821839083013668942322683_8, 0.767777432104826194917977340974503_8/
    Data X(15:16) /0.733790062453226804726171131369528_8, 0.697850494793315796932292388026640_8/
    Data X(17:18) /0.660061064126626961370053668149271_8, 0.620526182989242861140477556431189_8/
    Data X(19:20) /0.579345235826361691756024932172540_8, 0.536624148142019899264169793311073_8/
    Data X(21:22) /0.492480467861778574993693061207709_8, 0.447033769538089176780609900322854_8/
    Data X(23:24) /0.400401254830394392535476211542661_8, 0.352704725530878113471037207089374_8/
    Data X(25:26) /0.304073202273625077372677107199257_8, 0.254636926167889846439805129817805_8/
    Data X(27:28) /0.204525116682309891438957671002025_8, 0.153869913608583546963794672743256_8/
    Data X(29:30) /0.102806937966737030147096751318001_8, 0.051471842555317695833025213166723_8/
    Data X(31) /0._8/
    Data X(32:33) /-0.051471842555317695833025213166723_8, -0.102806937966737030147096751318001_8/
    Data X(34:35) /-0.153869913608583546963794672743256_8, -0.204525116682309891438957671002025_8/
    Data X(36:37) /-0.254636926167889846439805129817805_8, -0.304073202273625077372677107199257_8/
    Data X(38:39) /-0.352704725530878113471037207089374_8, -0.400401254830394392535476211542661_8/
    Data X(40:41) /-0.447033769538089176780609900322854_8, -0.492480467861778574993693061207709_8/
    Data X(42:43) /-0.536624148142019899264169793311073_8, -0.579345235826361691756024932172540_8/
    Data X(44:45) /-0.620526182989242861140477556431189_8, -0.660061064126626961370053668149271_8/
    Data X(46:47) /-0.697850494793315796932292388026640_8, -0.733790062453226804726171131369528_8/
    Data X(48:49) /-0.767777432104826194917977340974503_8, -0.799727835821839083013668942322683_8/
    Data X(50:51) /-0.829565762382768397442898119732502_8, -0.857205233546061098958658510658944_8/
    Data X(52:53) /-0.882560535792052681543116462530226_8, -0.905573307699907798546522558925958_8/
    Data X(54:55) /-0.926200047429274325879324277080474_8, -0.944374444748559979415831324037439_8/
    Data X(56:57) /-0.960021864968307512216871025581798_8, -0.973116322501126268374693868423707_8/
    Data X(58:59) /-0.983668123279747209970032581605663_8, -0.991630996870404594858628366109486_8/
    Data X(60:61) /-0.996893484074649540271630050918695_8, -0.999484410050490637571325895705811_8/
    IG=0._8
    IK=0._8
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,61
        Y=C1*X(Kstep)+C2
        Tmp=FuncS(Y, Razm)
        IK=IK+WK(Kstep)*Tmp
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Tmp*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine ArrIntGaussKron61Real

Recursive Subroutine ArrIntGaussKron61RealN(Llim, Hlim, Razb, LResult, Err, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Lstep, Razb, Razm
    Double Precision, Dimension(Razm):: PResult, LResult, Err, Perr
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Precision, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Err=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussKron61Real(Alim, Blim, PResult, Perr, FuncN, Razm)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussKron61Real(Alim, Blim, PResult, Perr, FuncN, Razm)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine ArrIntGaussKron61RealN

!=====================================================
!   Gauss-Kronrod method of 61 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntGaussKron61Comp(Alim, Blim, IK, Error, FuncS)
    Implicit None
    Double Precision:: WG(30), WK(61), X(61), Alim, Blim, Error, Y, C1, C2
    Double Complex:: IG, IK, FuncS
    Integer:: Kstep, Lstep, I
    External FuncS
    Data WG(1:2) /0.007968192496166605615465883474674_8, 0.018466468311090959142302131912047_8/
    Data WG(3:4) /0.028784707883323369349719179611292_8, 0.038799192569627049596801936446348_8/
    Data WG(5:6) /0.048402672830594052902938140422808_8, 0.057493156217619066481721689402056_8/
    Data WG(7:8) /0.065974229882180495128128515115962_8, 0.073755974737705206268243850022191_8/
    Data WG(9:10) /0.080755895229420215354694938460530_8, 0.086899787201082979802387530715126_8/
    Data WG(11:12) /0.092122522237786128717632707087619_8, 0.096368737174644259639468626351810_8/
    Data WG(13:14) /0.099593420586795267062780282103569_8, 0.101762389748405504596428952168554_8/
    Data WG(15:16) /0.102852652893558840341285636705415_8, 0.102852652893558840341285636705415_8/
    Data WG(17:18) /0.101762389748405504596428952168554_8, 0.099593420586795267062780282103569_8/
    Data WG(19:20) /0.096368737174644259639468626351810_8, 0.092122522237786128717632707087619_8/
    Data WG(21:22) /0.086899787201082979802387530715126_8, 0.080755895229420215354694938460530_8/
    Data WG(23:24) /0.073755974737705206268243850022191_8, 0.065974229882180495128128515115962_8/
    Data WG(25:26) /0.057493156217619066481721689402056_8, 0.048402672830594052902938140422808_8/
    Data WG(27:28) /0.038799192569627049596801936446348_8, 0.028784707883323369349719179611292_8/
    Data WG(29:30) /0.018466468311090959142302131912047_8, 0.007968192496166605615465883474674_8/   
    Data WK(1:2) /0.001389013698677007624551591226760_8, 0.003890461127099884051267201844516_8/
    Data WK(3:4) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(5:6) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(7:8) /0.016920889189053272627572289420322_8, 0.019414141193942381173408951050128_8/
    Data WK(9:10) /0.021828035821609192297167485738339_8, 0.024191162078080601365686370725232_8/
    Data WK(11:12) /0.026509954882333101610601709335075_8, 0.028754048765041292843978785354334_8/
    Data WK(13:14) /0.030907257562387762472884252943092_8, 0.032981447057483726031814191016854_8/
    Data WK(15:16) /0.034979338028060024137499670731468_8, 0.036882364651821229223911065617136_8/
    Data WK(17:18) /0.038678945624727592950348651532281_8, 0.040374538951535959111995279752468_8/
    Data WK(19:20) /0.041969810215164246147147541285970_8, 0.043452539701356069316831728117073_8/
    Data WK(21:22) /0.044814800133162663192355551616723_8, 0.046059238271006988116271735559374_8/
    Data WK(23:24) /0.047185546569299153945261478181099_8, 0.048185861757087129140779492298305_8/
    Data WK(25:26) /0.049055434555029778887528165367238_8, 0.049795683427074206357811569379942_8/
    Data WK(27:28) /0.050405921402782346840893085653585_8, 0.050881795898749606492297473049805_8/
    Data WK(29:30) /0.051221547849258772170656282604944_8, 0.051426128537459025933862879215781_8/
        Data WK(31) /0.051494729429451567558340433647099_8/
    Data WK(32:33) /0.051426128537459025933862879215781_8, 0.051221547849258772170656282604944_8/
    Data WK(34:35) /0.050881795898749606492297473049805_8, 0.050405921402782346840893085653585_8/
    Data WK(36:37) /0.049795683427074206357811569379942_8, 0.049055434555029778887528165367238_8/
    Data WK(38:39) /0.048185861757087129140779492298305_8, 0.047185546569299153945261478181099_8/
    Data WK(40:41) /0.046059238271006988116271735559374_8, 0.044814800133162663192355551616723_8/
    Data WK(42:43) /0.043452539701356069316831728117073_8, 0.041969810215164246147147541285970_8/
    Data WK(44:45) /0.040374538951535959111995279752468_8, 0.038678945624727592950348651532281_8/
    Data WK(46:47) /0.036882364651821229223911065617136_8, 0.034979338028060024137499670731468_8/
    Data WK(48:49) /0.032981447057483726031814191016854_8, 0.030907257562387762472884252943092_8/
    Data WK(50:51) /0.028754048765041292843978785354334_8, 0.026509954882333101610601709335075_8/
    Data WK(52:53) /0.024191162078080601365686370725232_8, 0.021828035821609192297167485738339_8/
    Data WK(54:55) /0.019414141193942381173408951050128_8, 0.016920889189053272627572289420322_8/
    Data WK(56:57) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(58:59) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(60:61) /0.003890461127099884051267201844516_8, 0.001389013698677007624551591226760_8/   
    Data X(1:2) /0.999484410050490637571325895705811_8, 0.996893484074649540271630050918695_8/
    Data X(3:4) /0.991630996870404594858628366109486_8, 0.983668123279747209970032581605663_8/
    Data X(5:6) /0.973116322501126268374693868423707_8, 0.960021864968307512216871025581798_8/
    Data X(7:8) /0.944374444748559979415831324037439_8, 0.926200047429274325879324277080474_8/
    Data X(9:10) /0.905573307699907798546522558925958_8, 0.882560535792052681543116462530226_8/
    Data X(11:12) /0.857205233546061098958658510658944_8, 0.829565762382768397442898119732502_8/
    Data X(13:14) /0.799727835821839083013668942322683_8, 0.767777432104826194917977340974503_8/
    Data X(15:16) /0.733790062453226804726171131369528_8, 0.697850494793315796932292388026640_8/
    Data X(17:18) /0.660061064126626961370053668149271_8, 0.620526182989242861140477556431189_8/
    Data X(19:20) /0.579345235826361691756024932172540_8, 0.536624148142019899264169793311073_8/
    Data X(21:22) /0.492480467861778574993693061207709_8, 0.447033769538089176780609900322854_8/
    Data X(23:24) /0.400401254830394392535476211542661_8, 0.352704725530878113471037207089374_8/
    Data X(25:26) /0.304073202273625077372677107199257_8, 0.254636926167889846439805129817805_8/
    Data X(27:28) /0.204525116682309891438957671002025_8, 0.153869913608583546963794672743256_8/
    Data X(29:30) /0.102806937966737030147096751318001_8, 0.051471842555317695833025213166723_8/
    Data X(31) /0._8/
    Data X(32:33) /-0.051471842555317695833025213166723_8, -0.102806937966737030147096751318001_8/
    Data X(34:35) /-0.153869913608583546963794672743256_8, -0.204525116682309891438957671002025_8/
    Data X(36:37) /-0.254636926167889846439805129817805_8, -0.304073202273625077372677107199257_8/
    Data X(38:39) /-0.352704725530878113471037207089374_8, -0.400401254830394392535476211542661_8/
    Data X(40:41) /-0.447033769538089176780609900322854_8, -0.492480467861778574993693061207709_8/
    Data X(42:43) /-0.536624148142019899264169793311073_8, -0.579345235826361691756024932172540_8/
    Data X(44:45) /-0.620526182989242861140477556431189_8, -0.660061064126626961370053668149271_8/
    Data X(46:47) /-0.697850494793315796932292388026640_8, -0.733790062453226804726171131369528_8/
    Data X(48:49) /-0.767777432104826194917977340974503_8, -0.799727835821839083013668942322683_8/
    Data X(50:51) /-0.829565762382768397442898119732502_8, -0.857205233546061098958658510658944_8/
    Data X(52:53) /-0.882560535792052681543116462530226_8, -0.905573307699907798546522558925958_8/
    Data X(54:55) /-0.926200047429274325879324277080474_8, -0.944374444748559979415831324037439_8/
    Data X(56:57) /-0.960021864968307512216871025581798_8, -0.973116322501126268374693868423707_8/
    Data X(58:59) /-0.983668123279747209970032581605663_8, -0.991630996870404594858628366109486_8/
    Data X(60:61) /-0.996893484074649540271630050918695_8, -0.999484410050490637571325895705811_8/
    IG=(0._8,0._8)
    IK=(0._8,0._8)
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,61
        Y=C1*X(Kstep)+C2
        Y=FuncS(Y)
        IK=IK+WK(Kstep)*Y
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Y*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine IntGaussKron61Comp

Recursive Subroutine IntGaussKron61CompN(Llim, Hlim, Razb, LResult, Err, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Err, Perr
    Double Complex:: PResult, LResult, FuncN
    Integer:: Lstep, Razb
    External FuncN
    LResult=(0._8,0._8)
    Err=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntGaussKron61Comp(Alim, Blim, PResult, Perr, FuncN)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntGaussKron61Comp(Alim, Blim, PResult, Perr, FuncN)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine IntGaussKron61CompN

!=====================================================
!   Gauss-Kronrod method of 61 order integration for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntGaussKron61Comp(Alim, Blim, IK, Error, FuncS, Razm)
    Implicit None
    Double Precision:: WG(30), WK(61), X(61), Alim, Blim, Y, C1, C2
    Integer:: Kstep, Lstep, I, Razm
    Double Complex, Dimension(Razm):: Tmp, IG, IK
    Double Precision, Dimension(Razm):: Error
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Data WG(1:2) /0.007968192496166605615465883474674_8, 0.018466468311090959142302131912047_8/
    Data WG(3:4) /0.028784707883323369349719179611292_8, 0.038799192569627049596801936446348_8/
    Data WG(5:6) /0.048402672830594052902938140422808_8, 0.057493156217619066481721689402056_8/
    Data WG(7:8) /0.065974229882180495128128515115962_8, 0.073755974737705206268243850022191_8/
    Data WG(9:10) /0.080755895229420215354694938460530_8, 0.086899787201082979802387530715126_8/
    Data WG(11:12) /0.092122522237786128717632707087619_8, 0.096368737174644259639468626351810_8/
    Data WG(13:14) /0.099593420586795267062780282103569_8, 0.101762389748405504596428952168554_8/
    Data WG(15:16) /0.102852652893558840341285636705415_8, 0.102852652893558840341285636705415_8/
    Data WG(17:18) /0.101762389748405504596428952168554_8, 0.099593420586795267062780282103569_8/
    Data WG(19:20) /0.096368737174644259639468626351810_8, 0.092122522237786128717632707087619_8/
    Data WG(21:22) /0.086899787201082979802387530715126_8, 0.080755895229420215354694938460530_8/
    Data WG(23:24) /0.073755974737705206268243850022191_8, 0.065974229882180495128128515115962_8/
    Data WG(25:26) /0.057493156217619066481721689402056_8, 0.048402672830594052902938140422808_8/
    Data WG(27:28) /0.038799192569627049596801936446348_8, 0.028784707883323369349719179611292_8/
    Data WG(29:30) /0.018466468311090959142302131912047_8, 0.007968192496166605615465883474674_8/   
    Data WK(1:2) /0.001389013698677007624551591226760_8, 0.003890461127099884051267201844516_8/
    Data WK(3:4) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(5:6) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(7:8) /0.016920889189053272627572289420322_8, 0.019414141193942381173408951050128_8/
    Data WK(9:10) /0.021828035821609192297167485738339_8, 0.024191162078080601365686370725232_8/
    Data WK(11:12) /0.026509954882333101610601709335075_8, 0.028754048765041292843978785354334_8/
    Data WK(13:14) /0.030907257562387762472884252943092_8, 0.032981447057483726031814191016854_8/
    Data WK(15:16) /0.034979338028060024137499670731468_8, 0.036882364651821229223911065617136_8/
    Data WK(17:18) /0.038678945624727592950348651532281_8, 0.040374538951535959111995279752468_8/
    Data WK(19:20) /0.041969810215164246147147541285970_8, 0.043452539701356069316831728117073_8/
    Data WK(21:22) /0.044814800133162663192355551616723_8, 0.046059238271006988116271735559374_8/
    Data WK(23:24) /0.047185546569299153945261478181099_8, 0.048185861757087129140779492298305_8/
    Data WK(25:26) /0.049055434555029778887528165367238_8, 0.049795683427074206357811569379942_8/
    Data WK(27:28) /0.050405921402782346840893085653585_8, 0.050881795898749606492297473049805_8/
    Data WK(29:30) /0.051221547849258772170656282604944_8, 0.051426128537459025933862879215781_8/
        Data WK(31) /0.051494729429451567558340433647099_8/
    Data WK(32:33) /0.051426128537459025933862879215781_8, 0.051221547849258772170656282604944_8/
    Data WK(34:35) /0.050881795898749606492297473049805_8, 0.050405921402782346840893085653585_8/
    Data WK(36:37) /0.049795683427074206357811569379942_8, 0.049055434555029778887528165367238_8/
    Data WK(38:39) /0.048185861757087129140779492298305_8, 0.047185546569299153945261478181099_8/
    Data WK(40:41) /0.046059238271006988116271735559374_8, 0.044814800133162663192355551616723_8/
    Data WK(42:43) /0.043452539701356069316831728117073_8, 0.041969810215164246147147541285970_8/
    Data WK(44:45) /0.040374538951535959111995279752468_8, 0.038678945624727592950348651532281_8/
    Data WK(46:47) /0.036882364651821229223911065617136_8, 0.034979338028060024137499670731468_8/
    Data WK(48:49) /0.032981447057483726031814191016854_8, 0.030907257562387762472884252943092_8/
    Data WK(50:51) /0.028754048765041292843978785354334_8, 0.026509954882333101610601709335075_8/
    Data WK(52:53) /0.024191162078080601365686370725232_8, 0.021828035821609192297167485738339_8/
    Data WK(54:55) /0.019414141193942381173408951050128_8, 0.016920889189053272627572289420322_8/
    Data WK(56:57) /0.011823015253496341742232898853251_8, 0.014369729507045804812451432443580_8/
    Data WK(58:59) /0.006630703915931292173319826369750_8, 0.009273279659517763428441146892024_8/
    Data WK(60:61) /0.003890461127099884051267201844516_8, 0.001389013698677007624551591226760_8/   
    Data X(1:2) /0.999484410050490637571325895705811_8, 0.996893484074649540271630050918695_8/
    Data X(3:4) /0.991630996870404594858628366109486_8, 0.983668123279747209970032581605663_8/
    Data X(5:6) /0.973116322501126268374693868423707_8, 0.960021864968307512216871025581798_8/
    Data X(7:8) /0.944374444748559979415831324037439_8, 0.926200047429274325879324277080474_8/
    Data X(9:10) /0.905573307699907798546522558925958_8, 0.882560535792052681543116462530226_8/
    Data X(11:12) /0.857205233546061098958658510658944_8, 0.829565762382768397442898119732502_8/
    Data X(13:14) /0.799727835821839083013668942322683_8, 0.767777432104826194917977340974503_8/
    Data X(15:16) /0.733790062453226804726171131369528_8, 0.697850494793315796932292388026640_8/
    Data X(17:18) /0.660061064126626961370053668149271_8, 0.620526182989242861140477556431189_8/
    Data X(19:20) /0.579345235826361691756024932172540_8, 0.536624148142019899264169793311073_8/
    Data X(21:22) /0.492480467861778574993693061207709_8, 0.447033769538089176780609900322854_8/
    Data X(23:24) /0.400401254830394392535476211542661_8, 0.352704725530878113471037207089374_8/
    Data X(25:26) /0.304073202273625077372677107199257_8, 0.254636926167889846439805129817805_8/
    Data X(27:28) /0.204525116682309891438957671002025_8, 0.153869913608583546963794672743256_8/
    Data X(29:30) /0.102806937966737030147096751318001_8, 0.051471842555317695833025213166723_8/
    Data X(31) /0._8/
    Data X(32:33) /-0.051471842555317695833025213166723_8, -0.102806937966737030147096751318001_8/
    Data X(34:35) /-0.153869913608583546963794672743256_8, -0.204525116682309891438957671002025_8/
    Data X(36:37) /-0.254636926167889846439805129817805_8, -0.304073202273625077372677107199257_8/
    Data X(38:39) /-0.352704725530878113471037207089374_8, -0.400401254830394392535476211542661_8/
    Data X(40:41) /-0.447033769538089176780609900322854_8, -0.492480467861778574993693061207709_8/
    Data X(42:43) /-0.536624148142019899264169793311073_8, -0.579345235826361691756024932172540_8/
    Data X(44:45) /-0.620526182989242861140477556431189_8, -0.660061064126626961370053668149271_8/
    Data X(46:47) /-0.697850494793315796932292388026640_8, -0.733790062453226804726171131369528_8/
    Data X(48:49) /-0.767777432104826194917977340974503_8, -0.799727835821839083013668942322683_8/
    Data X(50:51) /-0.829565762382768397442898119732502_8, -0.857205233546061098958658510658944_8/
    Data X(52:53) /-0.882560535792052681543116462530226_8, -0.905573307699907798546522558925958_8/
    Data X(54:55) /-0.926200047429274325879324277080474_8, -0.944374444748559979415831324037439_8/
    Data X(56:57) /-0.960021864968307512216871025581798_8, -0.973116322501126268374693868423707_8/
    Data X(58:59) /-0.983668123279747209970032581605663_8, -0.991630996870404594858628366109486_8/
    Data X(60:61) /-0.996893484074649540271630050918695_8, -0.999484410050490637571325895705811_8/
    IG=(0._8,0._8)
    IK=(0._8,0._8)
    C1=0.5_8*(Blim-Alim)
    C2=0.5_8*(Blim+Alim)
    Do Kstep=1,61
        Y=C1*X(Kstep)+C2
        Tmp=FuncS(Y, Razm)
        IK=IK+WK(Kstep)*Tmp
        If (Mod(Kstep,2).EQ.0) Then
            IG=IG+Tmp*WG(Kstep/2)
        EndIf
    EndDo
    IG=IG*C1
    IK=IK*C1
    Error=Abs(IG-IK)
Return
End Subroutine ArrIntGaussKron61Comp

Recursive Subroutine ArrIntGaussKron61CompN(Llim, Hlim, Razb, LResult, Err, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength
    Integer:: Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Double Precision, Dimension(Razm):: Err, Perr
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=(0._8,0._8)
    Err=(0._8,0._8)
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntGaussKron61Comp(Alim, Blim, PResult, Perr, FuncN, Razm)
        LResult=LResult+PResult
        Err=Err+Perr
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntGaussKron61Comp(Alim, Blim, PResult, Perr, FuncN, Razm)
    LResult=LResult+PResult
    Err=Err+Perr
Return
End Subroutine ArrIntGaussKron61CompN

!=====================================================
!   Trapezoidal Filon method.
!=====================================================
Recursive Subroutine IntFilonTrap(Llim, Hlim, Frac, W, Res, Func)
    Implicit None
    Double Precision:: Llim, Hlim, Step, Xcurr, W, TmpR, TmpRR
    Double Complex:: Res, Tmp, I=(0._8,1._8)
    Integer:: Frac, Cnt
    Interface
        Double Complex Function Func(X)
            Double Precision:: X
        End Function Func
    End Interface
    Res=(0._8, 0._8)
    Step = (Hlim - Llim)/Dble(Frac)
    Xcurr = Llim
    Do Cnt=1, Frac-1
        Xcurr = Xcurr + Step
        Res = Res + Func(Xcurr)
    EndDo
    TmpR = 0.5_8*W*Step
    TmpRR = sin(TmpR)/TmpR
    Res = Res*TmpRR*TmpRR*Step
    Tmp = ((1._8,0._8)/(I*W) + ((1._8,0._8) - Exp(-I*W*Step))/(W*W*Step))
    Res = Res + Tmp*Func(Hlim)
    Tmp = ((1._8,0._8)/(I*W) + (Exp(I*W*Step) - (1._8,0._8))/(W*W*Step))
    Res = Res + Tmp*Func(Llim)
    Return
End Subroutine IntFilonTrap

Recursive Subroutine ArrIntFilonTrap(Llim, Hlim, Frac, W, Res, FuncS, Razm)
    Implicit None
    Double Precision:: Llim, Hlim, Step, Xcurr, W, TmpR, TmpRR
    Integer:: Frac, Cnt, Razm
    Double Complex:: Tmp, I=(0._8,1._8)
    Double Complex, Dimension(Razm):: Res
    Interface
        Function FuncS(Z, N)
            Integer:: N
            Double Precision:: Z
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Res=(0._8, 0._8)
    Step = (Hlim - Llim)/Dble(Frac)
    Xcurr = Llim
    Do Cnt=1, Frac-1
        Xcurr = Xcurr + Step
        Res = Res + FuncS(Xcurr, Razm)
    EndDo
    TmpR = 0.5_8*W*Step
    TmpRR = sin(TmpR)/TmpR
    Res = Res*TmpRR*TmpRR*Step
    Tmp = ((1._8,0._8)/(I*W) + ((1._8,0._8) - Exp(-I*W*Step))/(W*W*Step))
    Res = Res + Tmp*FuncS(Hlim, Razm)
    Tmp = ((1._8,0._8)/(I*W) + (Exp(I*W*Step) - (1._8,0._8))/(W*W*Step))
    Res = Res + Tmp*FuncS(Llim, Razm)
    Return
End Subroutine ArrIntFilonTrap

!=====================================================
!   Iserles method of 16 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntIserles16Comp(Llim, Hlim, Woo, ResL, Func)
    Use Matrix
    Implicit None
    Double Precision:: Llim, Hlim, Wo, Woo, YL
    Double Complex:: WL(16), BJL(16,17), Iw, Eiw, ResL
    Integer:: Kstep, Cnt
    Interface
        Double Complex Function Func(X)
            Double Precision:: X
        End Function Func
    End Interface
    Wo=Woo*(Hlim-Llim)
!   Weight matrix filling
!   Right parts calculation
    Iw=(0._8,1._8)*Wo
    Eiw=Exp(Iw)
    Iw=(1._8,0._8)/Iw
    BJL(1,17) = (Eiw-(1._8, 0._8))*Iw
    Do Kstep=2,16
        BJL(Kstep,17)=(Eiw-(Kstep-1)*BJL(Kstep-1,17))*Iw
    EndDo
!   Left parts filling  
    Do Cnt=1,16
        BJL(1,Cnt)=(1._8,0._8)
    EndDo
    Do Kstep=2,16
        Do Cnt=1,16
            BJL(Kstep,Cnt)=BJL(Kstep-1,Cnt)*XL16(Cnt)
        EndDo
    EndDo
!   Weights resolving
    Call SlaeGauss(BJL,WL)
!   Quadrature usage
    ResL=(0._8, 0._8)
    Do Kstep=1, 16
        YL=((Hlim-Llim)*XL16(Kstep)+Llim)
        ResL=ResL+WL(Kstep)*Func(YL)*Exp(-(0._8,1._8)*Woo*YL)
    EndDo
!   Result normalizing
    ResL=ResL*(Hlim-Llim)*Exp((0._8,1._8)*Woo*Llim)
Return
End Subroutine IntIserles16Comp

Recursive Subroutine IntIserles16CompN(Llim, Hlim, Razb, Woo, LResult, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Woo
    Double Complex:: PResult, LResult
    Integer:: Kstep, Lstep, Razb
    Interface
        Double Complex Function FuncN(X)
            Double Precision:: X
        End Function FuncN
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntIserles16Comp(Alim, Blim, Woo, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntIserles16Comp(Alim, Blim, Woo, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntIserles16CompN

!=====================================================
!   Iserles method of 16 order integration for double complex array functions.
!=====================================================

Recursive Subroutine ArrIntIserles16Comp(Llim, Hlim, Woo, ResL, FuncS, Razm)
    Use Matrix
    Implicit None
    Double Precision:: Llim, Hlim, Wo, Woo, YL
    Double Complex:: WL(16), BJL(16,17), Iw, Eiw
    Integer:: Kstep, Cnt, Razm
    Double Complex, Dimension(Razm):: ResL
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Wo=Woo*(Hlim-Llim)
!   Weight matrix filling
!   Right parts calculation
    Iw=(0._8,1._8)*Wo
    Eiw=Exp(Iw)
    Iw=(1._8,0._8)/Iw
    BJL(1,17) = (Eiw-(1._8, 0._8))*Iw
    Do Kstep=2,16
        BJL(Kstep,17)=(Eiw-(Kstep-1)*BJL(Kstep-1,17))*Iw
    EndDo
!   Left parts filling  
    Do Cnt=1,16
        BJL(1,Cnt)=(1._8,0._8)
    EndDo
    Do Kstep=2,16
        Do Cnt=1,16
            BJL(Kstep,Cnt)=BJL(Kstep-1,Cnt)*XL16(Cnt)
        EndDo
    EndDo
!   Weights resolving
    Call SlaeGauss(BJL,WL)
!   Quadrature usage
    ResL=(0._8, 0._8)
    Do Kstep=1, 16
        YL=((Hlim-Llim)*XL16(Kstep)+Llim)
        ResL=ResL+WL(Kstep)*FuncS(YL, Razm)*Exp(-(0._8,1._8)*Woo*YL)
    EndDo
!   Result normalizing
    ResL=ResL*(Hlim-Llim)*Exp((0._8,1._8)*Woo*Llim)
Return
End Subroutine ArrIntIserles16Comp

Recursive Subroutine ArrIntIserles16CompN(Llim, Hlim, Razb, Woo, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Woo
    Integer:: Kstep, Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntIserles16Comp(Alim, Blim, Woo, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntIserles16Comp(Alim, Blim, Woo, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntIserles16CompN

!=====================================================
!   Iserles method of 40 order integration for double complex functions.
!=====================================================
Recursive Subroutine IntIserles40Comp(Llim, Hlim, Woo, ResL, Func)
    Use Matrix
    Implicit None
    Double Precision:: Llim, Hlim, Wo, Woo, YL
    Double Complex:: WL(40), BJL(40,41), Iw, Eiw, ResL
    Integer:: Kstep, Cnt
    Interface
        Double Complex Function Func(X)
            Double Precision:: X
        End Function Func
    End Interface
    Wo=Woo*(Hlim-Llim)
!   Weight matrix filling
!   Right parts calculation
    Iw=(0._8,1._8)*Wo
    Eiw=Exp(Iw)
    Iw=(1._8,0._8)/Iw
    BJL(1,41) = (Eiw-(1._8, 0._8))*Iw
    Do Kstep=2,40
        BJL(Kstep,41)=(Eiw-(Kstep-1)*BJL(Kstep-1,41))*Iw
    EndDo
!   Left parts filling  
    Do Cnt=1,40
        BJL(1,Cnt)=(1._8,0._8)
    EndDo
    Do Kstep=2,40
        Do Cnt=1,40
            BJL(Kstep,Cnt)=BJL(Kstep-1,Cnt)*XL40(Cnt)
        EndDo
    EndDo
!   Weights resolving
    Call SlaeGauss(BJL,WL)
!   Quadrature usage
    ResL=(0._8, 0._8)
    Do Kstep=1, 40
        YL=((Hlim-Llim)*XL40(Kstep)+Llim)
        ResL=ResL+WL(Kstep)*Func(YL)*Exp(-(0._8,1._8)*Woo*YL)
    EndDo
!   Result normalizing
    ResL=ResL*(Hlim-Llim)*Exp((0._8,1._8)*Woo*Llim)
Return
End Subroutine IntIserles40Comp

Recursive Subroutine IntIserles40CompN(Llim, Hlim, Razb, Woo, LResult, FuncN)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Woo
    Double Complex:: PResult, LResult
    Integer:: Kstep, Lstep, Razb
    Interface
        Double Complex Function FuncN(X)
            Double Precision:: X
        End Function FuncN
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call IntIserles40Comp(Alim, Blim, Woo, PResult, FuncN)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call IntIserles40Comp(Alim, Blim, Woo, PResult, FuncN)
    LResult=LResult+PResult
Return
End Subroutine IntIserles40CompN

!=====================================================
!   Iserles method of 40 order integration for double complex array functions.
!=====================================================

Recursive Subroutine ArrIntIserles40Comp(Llim, Hlim, Woo, ResL, FuncS, Razm)
    Use Matrix
    Implicit None
    Double Precision:: Llim, Hlim, Wo, Woo, YL
    Double Complex:: WL(40), BJL(40,41), Iw, Eiw
    Integer:: Kstep, Cnt, Razm
    Double Complex, Dimension(Razm):: ResL
    Interface
        Function FuncS(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncS
        End Function FuncS
    End Interface
    Wo=Woo*(Hlim-Llim)
!   Weight matrix filling
!   Right parts calculation
    Iw=(0._8,1._8)*Wo
    Eiw=Exp(Iw)
    Iw=(1._8,0._8)/Iw
    BJL(1,41) = (Eiw-(1._8, 0._8))*Iw
    Do Kstep=2,40
        BJL(Kstep,41)=(Eiw-(Kstep-1)*BJL(Kstep-1,41))*Iw
    EndDo
!   Left parts filling  
    Do Cnt=1,40
        BJL(1,Cnt)=(1._8,0._8)
    EndDo
    Do Kstep=2,40
        Do Cnt=1,40
            BJL(Kstep,Cnt)=BJL(Kstep-1,Cnt)*XL40(Cnt)
        EndDo
    EndDo
!   Weights resolving
    Call SlaeGauss(BJL,WL)
!   Quadrature usage
    ResL=(0._8, 0._8)
    Do Kstep=1, 40
        YL=((Hlim-Llim)*XL40(Kstep)+Llim)
        ResL=ResL+WL(Kstep)*FuncS(YL, Razm)*Exp(-(0._8,1._8)*Woo*YL)
    EndDo
!   Result normalizing
    ResL=ResL*(Hlim-Llim)*Exp((0._8,1._8)*Woo*Llim)
Return
End Subroutine ArrIntIserles40Comp

Recursive Subroutine ArrIntIserles40CompN(Llim, Hlim, Razb, Woo, LResult, FuncN, Razm)
    Implicit None
    Double Precision:: Hlim, Llim, Alim, Blim, Plength, Woo
    Integer:: Kstep, Lstep, Razb, Razm
    Double Complex, Dimension(Razm):: PResult, LResult
    Interface
        Function FuncN(X, N)
            Integer:: N
            Double Precision:: X
            Double Complex, Dimension(N):: FuncN
        End Function FuncN
    End Interface
    LResult=0._8
    Plength=(Hlim-Llim)/Real(Razb)
    Alim=Llim
    Do LStep=1,Razb-1
        Blim=Alim+Plength
        Call ArrIntIserles40Comp(Alim, Blim, Woo, PResult, FuncN, Razm)
        LResult=LResult+PResult
        Alim=Blim
    EndDo
    Blim=Hlim
    Call ArrIntIserles40Comp(Alim, Blim, Woo, PResult, FuncN, Razm)
    LResult=LResult+PResult
Return
End Subroutine ArrIntIserles40CompN

!=====================================================
!   Simpson-type 5-point two-dimensional cubature for double precision functions.
!=====================================================
Recursive Subroutine IntCubature5Real(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res1, Func)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, StepX, StepY
    Double Precision:: Res1
    Integer:: FracX, FracY, CntX, CntY
    Interface
        Double Precision Function Func(X, Y)
            Double Precision:: X, Y
        End Function Func
    End Interface
    Res1 = 0._8
    StepX = (HlimX-LlimX)/Dble(FracX)
    StepY = (HlimY-LlimY)/Dble(FracY)
    Xcurr = LlimX + StepX*0.5_8
    Do CntX=1, 2*FracX-1
        Ycurr = LlimY + StepY*0.5_8
        Do CntY=1, 2*FracY-1
            If (Mod(CntX,2).EQ.0.AND.Mod(CntY,2).EQ.0) Then
            Else
                Res1 = Res1 + Func(Xcurr, Ycurr)
            EndIf
            Ycurr = Ycurr + StepY*0.5_8
        EndDo
        Xcurr = Xcurr + StepX*0.5_8
    EndDo
    Res1 = 2._8*Res1
    Xcurr = LlimX + StepX*0.5_8
    Ycurr = LlimY + StepY*0.5_8
    Do CntX=1, FracX
        Res1= Res1 + Func(Xcurr, LlimY)+Func(Xcurr, HlimY)
        Xcurr = Xcurr + StepX
    EndDo
    Do CntY=1, FracY
        Res1=Res1+Func(LlimY, Ycurr)+Func(HlimY, Ycurr)
        Ycurr = Ycurr + StepY
    EndDo
    Res1 = Res1*StepX*StepY*0.16666666666666666666666666666666666667_8
Return
End Subroutine IntCubature5Real

!=====================================================
!   Simpson-type 5-point two-dimensional cubature for double precision array functions.
!=====================================================
Recursive Subroutine ArrIntCubature5Real(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res1, Func, N)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, StepX, StepY
    Double Precision:: Res1(N)
    Integer:: FracX, FracY, CntX, CntY, N
    Interface
        Function Func(X, Y, N)
            Double Precision:: X, Y
            Integer:: N
            Double Precision:: Func(N)
        End Function Func
    End Interface
    Res1 = 0._8
    StepX = (HlimX-LlimX)/Dble(FracX)
    StepY = (HlimY-LlimY)/Dble(FracY)
    Xcurr = LlimX + StepX*0.5_8
    Do CntX=1, 2*FracX-1
        Ycurr = LlimY + StepY*0.5_8
        Do CntY=1, 2*FracY-1
            If (Mod(CntX,2).EQ.0.AND.Mod(CntY,2).EQ.0) Then
            Else
                Res1 = Res1 + Func(Xcurr, Ycurr, N)
            EndIf
            Ycurr = Ycurr + StepY*0.5_8
        EndDo
        Xcurr = Xcurr + StepX*0.5_8
    EndDo
    Res1 = 2._8*Res1
    Xcurr = LlimX + StepX*0.5_8
    Ycurr = LlimY + StepY*0.5_8
    Do CntX=1, FracX
        Res1= Res1 + Func(Xcurr, LlimY, N)+Func(Xcurr, HlimY, N)
        Xcurr = Xcurr + StepX
    EndDo
    Do CntY=1, FracY
        Res1=Res1+Func(LlimY, Ycurr, N)+Func(HlimY, Ycurr, N)
        Ycurr = Ycurr + StepY
    EndDo
    Res1 = Res1*StepX*StepY*0.16666666666666666666666666666666666667_8
Return
End Subroutine ArrIntCubature5Real

!=====================================================
!   Simpson-type 5-point two-dimensional cubature for double precision functions.
!=====================================================
Recursive Subroutine IntCubature5Comp(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res1, Func)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, StepX, StepY
    Double complex:: Res1
    Integer:: FracX, FracY, CntX, CntY
    Interface
        Double Complex Function Func(X, Y)
            Double Precision:: X, Y
        End Function Func
    End Interface
    Res1 = (0._8, 0._8)
    StepX = (HlimX-LlimX)/Dble(FracX)
    StepY = (HlimY-LlimY)/Dble(FracY)
    Xcurr = LlimX + StepX*0.5_8
    Do CntX=1, 2*FracX-1
        Ycurr = LlimY + StepY*0.5_8
        Do CntY=1, 2*FracY-1
            If (Mod(CntX,2).EQ.0.AND.Mod(CntY,2).EQ.0) Then
            Else
                Res1 = Res1 + Func(Xcurr, Ycurr)
            EndIf
            Ycurr = Ycurr + StepY*0.5_8
        EndDo
        Xcurr = Xcurr + StepX*0.5_8
    EndDo
    Res1 = 2._8*Res1
    Xcurr = LlimX + StepX*0.5_8
    Ycurr = LlimY + StepY*0.5_8
    Do CntX=1, FracX
        Res1= Res1 + Func(Xcurr, LlimY)+Func(Xcurr, HlimY)
        Xcurr = Xcurr + StepX
    EndDo
    Do CntY=1, FracY
        Res1=Res1+Func(LlimY, Ycurr)+Func(HlimY, Ycurr)
        Ycurr = Ycurr + StepY
    EndDo
    Res1 = Res1*StepX*StepY*(0.16666666666666666666666666666666666667_8, 0._8)
Return
End Subroutine IntCubature5Comp
!=====================================================
!   Simpson-type 5-point two-dimensional cubature for double complex array functions.
!=====================================================
Recursive Subroutine ArrIntCubature5Comp(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res1, Func, N)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, StepX, StepY
    Double Complex:: Res1(N)
    Integer:: FracX, FracY, CntX, CntY, N
    Interface
        Function Func(X, Y, N)
            Double Precision:: X, Y
            Integer:: N
            Double Complex:: Func(N)
        End Function Func
    End Interface
    Res1 = (0._8, 0._8)
    StepX = (HlimX-LlimX)/Dble(FracX)
    StepY = (HlimY-LlimY)/Dble(FracY)
    Xcurr = LlimX + StepX*0.5_8
    Do CntX=1, 2*FracX-1
        Ycurr = LlimY + StepY*0.5_8
        Do CntY=1, 2*FracY-1
            If (Mod(CntX,2).EQ.0.AND.Mod(CntY,2).EQ.0) Then
            Else
                Res1 = Res1 + Func(Xcurr, Ycurr, N)
            EndIf
            Ycurr = Ycurr + StepY*0.5_8
        EndDo
        Xcurr = Xcurr + StepX*0.5_8
    EndDo
    Res1 = 2._8*Res1
    Xcurr = LlimX + StepX*0.5_8
    Ycurr = LlimY + StepY*0.5_8
    Do CntX=1, FracX
        Res1= Res1 + Func(Xcurr, LlimY, N)+Func(Xcurr, HlimY, N)
        Xcurr = Xcurr + StepX
    EndDo
    Do CntY=1, FracY
        Res1=Res1+Func(LlimY, Ycurr, N)+Func(HlimY, Ycurr, N)
        Ycurr = Ycurr + StepY
    EndDo
    Res1 = Res1*StepX*StepY*0.16666666666666666666666666666666666667_8
Return
End Subroutine ArrIntCubature5Comp

!=====================================================
!   Gauss-type 40-node two-dimensional cubature for double precision functions.
!=====================================================
Subroutine IntGauss40CubatureReal(LlimX, HlimX, LlimY, HlimY, Res, Func)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr
    Double Precision:: Res
    Integer:: KsX, KsY
    Interface
        Double Precision Function Func(X, Y)
            Double Precision:: X, Y
        End Function Func
    End Interface
    Res=0._8
    Do KsX=1,40
        Xcurr=0.5_8*((HlimX-LlimX)*XG40(KsX)+HlimX+LlimX)
        Do KsY=1,40
            Ycurr=0.5_8*((HlimY-LlimY)*XG40(KsY)+HlimY+LlimY)
            Res=Res + WG40(KsX)*WG40(KsY)*Func(Xcurr, Ycurr)
        EndDo
    EndDo
    Res=Res*(HlimX-LlimX)*(HlimY-LlimY)*0.25_8
Return
End Subroutine IntGauss40CubatureReal

Subroutine IntGauss40CubatureRealN(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res, Func)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, BlimX, AlimX, BlimY, AlimY, StepX, StepY
    Double Precision:: Res, Pres
    Integer:: FracX, FracY, CntX, CntY
    Interface
        Double Precision Function Func(X, Y)
            Double Precision:: X, Y
        End Function Func
    End Interface
    Res=0._8
    StepX=(HlimX-LlimX)/Real(FracX)
    StepY=(HlimY-LlimY)/Real(FracY)
    AlimX=LlimX
    BlimX=AlimX+StepX
    Do CntX=1, FracX
        AlimY=LlimY
        BlimY=AlimY+StepY
        Do CntY=1, FracY
            Call IntGauss40CubatureReal(AlimX, BlimX, AlimY, BlimY, Pres, Func)
            Res = Res + Pres
            AlimY=BlimY
            BlimY=BlimY+StepY
        EndDo
        AlimX=BlimX
        BlimX=BlimX+StepX
    EndDo
Return
End Subroutine IntGauss40CubatureRealN

!=====================================================
!   Gauss-type 40-node two-dimensional cubature for double precision arra functions.
!=====================================================
Subroutine ArrIntGauss40CubatureReal(LlimX, HlimX, LlimY, HlimY, Res, Func, N)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr
    Integer:: N, KsX, KsY
    Double Precision:: Res(N)
    Interface
        Function Func(X, Y, N)
            Double Precision:: X, Y
            Integer:: N
            Double Precision:: Func(N)
        End Function Func
    End Interface
    Res=0._8
    Do KsX=1,40
        Xcurr=0.5_8*((HlimX-LlimX)*XG40(KsX)+HlimX+LlimX)
        Do KsY=1,40
            Ycurr=0.5_8*((HlimY-LlimY)*XG40(KsY)+HlimY+LlimY)
            Res=Res + WG40(KsX)*WG40(KsY)*Func(Xcurr, Ycurr, N)
        EndDo
    EndDo
    Res=Res*(HlimX-LlimX)*(HlimY-LlimY)*0.25_8
Return
End Subroutine ArrIntGauss40CubatureReal

Subroutine ArrIntGauss40CubatureRealN(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res, Func, N)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, BlimX, AlimX, BlimY, AlimY, StepX, StepY
    Integer:: FracX, FracY, CntX, CntY, N
    Double Precision:: Res(N), Pres(N)
    Interface
        Function Func(X, Y, N)
            Double Precision:: X, Y
            Integer:: N
            Double Precision:: Func(N)
        End Function Func
    End Interface
    Res=0._8
    StepX=(HlimX-LlimX)/Real(FracX)
    StepY=(HlimY-LlimY)/Real(FracY)
    AlimX=LlimX
    BlimX=AlimX+StepX
    Do CntX=1, FracX
        AlimY=LlimY
        BlimY=AlimY+StepY
        Do CntY=1, FracY
            Call ArrIntGauss40CubatureReal(AlimX, BlimX, AlimY, BlimY, Pres, Func, N)
            Res = Res + Pres
            AlimY=BlimY
            BlimY=BlimY+StepY
        EndDo
        AlimX=BlimX
        BlimX=BlimX+StepX
    EndDo
Return
End Subroutine ArrIntGauss40CubatureRealN

!=====================================================
!   Gauss-type 40-node two-dimensional cubature for double complex functions.
!=====================================================
Subroutine IntGauss40CubatureComp(LlimX, HlimX, LlimY, HlimY, Res, Func)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr
    Double Complex:: Res
    Integer:: KsX, KsY
    Interface
        Double Complex Function Func(X, Y)
            Double Precision:: X, Y
        End Function Func
    End Interface
    Res=(0._8, 0._8)
    Do KsX=1,40
        Xcurr=0.5_8*((HlimX-LlimX)*XG40(KsX)+HlimX+LlimX)
        Do KsY=1,40
            Ycurr=0.5_8*((HlimY-LlimY)*XG40(KsY)+HlimY+LlimY)
            Res=Res + WG40(KsX)*WG40(KsY)*Func(Xcurr, Ycurr)
        EndDo
    EndDo
    Res=Res*(HlimX-LlimX)*(HlimY-LlimY)*(0.25_8, 0._8)
Return
End Subroutine IntGauss40CubatureComp

Subroutine IntGauss40CubatureCompN(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res, Func)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, BlimX, AlimX, BlimY, AlimY, StepX, StepY
    Double Complex:: Res, Pres
    Integer:: FracX, FracY, CntX, CntY
    Interface
        Double Complex Function Func(X, Y)
            Double Precision:: X, Y
        End Function Func
    End Interface
    Res=(0._8, 0._8)
    StepX=(HlimX-LlimX)/Real(FracX)
    StepY=(HlimY-LlimY)/Real(FracY)
    AlimX=LlimX
    BlimX=AlimX+StepX
    Do CntX=1, FracX
        AlimY=LlimY
        BlimY=AlimY+StepY
        Do CntY=1, FracY
            Call IntGauss40CubatureComp(AlimX, BlimX, AlimY, BlimY, Pres, Func)
            Res = Res + Pres
            AlimY=BlimY
            BlimY=BlimY+StepY
        EndDo
        AlimX=BlimX
        BlimX=BlimX+StepX
    EndDo
Return
End Subroutine IntGauss40CubatureCompN

!=====================================================
!   Gauss-type 40-node two-dimensional cubature for double precision arra functions.
!=====================================================
Subroutine ArrIntGauss40CubatureComp(LlimX, HlimX, LlimY, HlimY, Res, Func, N)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr
    Integer:: N, KsX, KsY
    Double Complex:: Res(N)
    Interface
        Function Func(X, Y, N)
            Double Precision:: X, Y
            Integer:: N
            Double Complex:: Func(N)
        End Function Func
    End Interface
    Res=(0._8, 0._8)
    Do KsX=1,40
        Xcurr=0.5_8*((HlimX-LlimX)*XG40(KsX)+HlimX+LlimX)
        Do KsY=1,40
            Ycurr=0.5_8*((HlimY-LlimY)*XG40(KsY)+HlimY+LlimY)
            Res=Res + WG40(KsX)*WG40(KsY)*Func(Xcurr, Ycurr, N)
        EndDo
    EndDo
    Res=Res*(HlimX-LlimX)*(HlimY-LlimY)*(0.25_8, 0._8)
Return
End Subroutine ArrIntGauss40CubatureComp

Subroutine ArrIntGauss40CubatureCompN(LlimX, HlimX, LlimY, HlimY, FracX, FracY, Res, Func, N)
    Implicit None
    Double Precision:: LlimX, HlimX, LlimY, HlimY, Xcurr, Ycurr, BlimX, AlimX, BlimY, AlimY, StepX, StepY
    Integer:: FracX, FracY, CntX, CntY, N
    Double Complex:: Res(N), Pres(N)
    Interface
        Function Func(X, Y, N)
            Double Precision:: X, Y
            Integer:: N
            Double Complex:: Func(N)
        End Function Func
    End Interface
    Res=(0._8, 0._8)
    StepX=(HlimX-LlimX)/Real(FracX)
    StepY=(HlimY-LlimY)/Real(FracY)
    AlimX=LlimX
    BlimX=AlimX+StepX
    Do CntX=1, FracX
        AlimY=LlimY
        BlimY=AlimY+StepY
        Do CntY=1, FracY
            Call ArrIntGauss40CubatureComp(AlimX, BlimX, AlimY, BlimY, Pres, Func, N)
            Res = Res + Pres
            AlimY=BlimY
            BlimY=BlimY+StepY
        EndDo
        AlimX=BlimX
        BlimX=BlimX+StepX
    EndDo
Return
End Subroutine ArrIntGauss40CubatureCompN

End module Integrals
Module Specfun

!=====================================================
!   Interface for subprograms which calculate the Bessel functions of first kind
!   X=BesJ(o, Z)
!   where o - is the function order (integer or double precision);
!   Z - is the argument (double precision or double complex).
!=====================================================
Interface BesJ
    Module Procedure JioReal, JroReal, JioComp, JroComp
End Interface BesJ
!=====================================================
!   Interface for subprograms which calculate the Bessel functions of second kind
!   (Newman function)
!   X=BesN(o, Z)
!   where o - is the function order (integer or double precision);
!   Z - is the argument (double precision or double complex).
!=====================================================
Interface BesN
    Module Procedure NioReal, NroReal, NioComp, NroComp
End Interface BesN
!=====================================================
!   Interface for subprograms which calculate the Hankel functions of first kind
!   X=BesH1(o, Z)
!   where o - is the function order (integer or double precision);
!   Z - is the argument (double precision or double complex).
!=====================================================
Interface BesH1
    Module Procedure H1ioReal, H1ioComp
End Interface BesH1
!=====================================================
!   Interface for subprograms which calculate the Hankel functions of first kind
!   X=BesH2(o, Z)
!   where o - is the function order (integer or double precision);
!   Z - is the argument (double precision or double complex).
!=====================================================
Interface BesH2
    Module Procedure H2ioReal, H2ioComp
End Interface BesH2
!=====================================================
!   Interface for subprograms which calculate the factorial
!   K=Factorial(N)
!   where N - is the integer argument.
!=====================================================
Interface Factorial
    Module Procedure Factorial
End Interface Factorial
!=====================================================
!   Interface for subprograms which calculate the gamma-function
!   K=Gamma(N)
!   where N - is the double precision argument.
!=====================================================
Interface Gamma
    Module Procedure GammaReal, GammaComp
End Interface Gamma
!=====================================================
!   Interface for subprograms which calculate the digamma-function
!   K=Digamma(N)
!   where N - is the integer argument.
!=====================================================
Interface Digamma
    Module Procedure DigammaInt
End Interface Digamma

!=====================================================
!   Auxillary interfaces.
!=====================================================


!=====================================================
!   Module procedures.
!=====================================================
Contains
!=====================================================
!   The factorial
!=====================================================
Integer Function Factorial(N)
    Implicit None
    Integer N, K
    Factorial=1
    If (N.NE.0) Then
        Do K=N,1,-1
            Factorial=Factorial*K
        EndDo
    EndIf
End Function Factorial

!=====================================================
!   Gamma function for double precision arguments.
!=====================================================
Double Precision Function GammaReal(X)
    Implicit None
    Double Precision:: X, Y, Z, Tmp
    Integer:: I
    Y = X
    Y = Y+20._8
    Z = Y
    !Stirling formula usage:
    Y=1._8/Y
    Tmp=1._8+Y*(0.0833333333333333_8+Y*(0.00347222222222222_8-Y*(0.002681327160493827_8+Y*0.0002294720936213992_8)))
    Tmp=Tmp*Exp(-Z)*(Z**(Z-0.5_8))*2.50662827463100045_8
    !Back Recursion
    Do I=1,20
        Tmp=Tmp/(Z-I)
    EndDo
    GammaReal=Tmp
End Function GammaReal

Double Complex Function GammaComp(Z)
    Implicit None
    Double Complex:: Z
    GammaComp=(0._8,0._8)
End Function GammaComp

!=====================================================
!   Digamma function for integer arguments.
!=====================================================
Double Precision Function DigammaInt(N)
    Implicit None
    Integer:: N, K
    DigammaInt=-0.577215664901532860606512_8
    Do K=1,N-1
        DigammaInt=DigammaInt+1._8/Real(K)
    EndDo
End Function DigammaInt

!=====================================================
!   Bessel functions of integer order and double precision argument
!=====================================================
Double Precision Function JioReal(p,x)
    Implicit None
    Integer:: p, k, p2
    Double Precision:: Ck, Tmp, X, xp2
    Double Complex:: cCk, cTmp, cHio
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=0.5_8*x
        Ck=1._8/Gamma(p+1._8)
        JioReal = Ck
        Do k=1,35
            Ck = -Ck*xp2*xp2/(k*(k+p))
            JioReal = JioReal+Ck
        EndDo
        JioReal=JioReal*(xp2**p)
    Else
        Tmp=Dble(p)
        JioReal=Real(AsyHank1(Tmp,Dcmplx(x,0._8)))
    EndIf
Return
End Function JioReal

!=====================================================
!   Bessel functions of real order and double precision argument
!=====================================================
Double Precision Function JroReal(p,x)
    Implicit None
    Integer:: k
    Double Precision:: Ck, Tmp, X, xp2,  p, p2
    Double Complex:: cCk, cTmp, cHio
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=0.5_8*x
        Ck=1._8/Gamma(p+1._8)
        JroReal = Ck
        Do k=1,35
            Ck = -Ck*xp2*xp2/(k*(k+p))
            JroReal = JroReal+Ck
        EndDo
        JroReal=JroReal*(xp2**p)
    Else
        JroReal=Real(AsyHank1(p,Dcmplx(x,0._8)))
    EndIf
Return
End Function JroReal

!=====================================================
!   Bessel functions of integer order and double complex argument
!=====================================================
Double Complex Function JioComp(p,x)
    Implicit None
    Integer:: p, k, p2
    Double Complex:: Ck, Tmp, X, xp2
    Double Complex:: cCk, cTmp, cHio
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=(0.5_8,0._8)*x
        Ck=1._8/Gamma(p+1._8)
        JioComp = Ck
        Do k=1,35
            Ck = -Ck*xp2*xp2/(k*(k+p))
            JioComp = JioComp+Ck
        EndDo
        JioComp=JioComp*(xp2**p)
    Else
        Tmp=Dble(p)
        JioComp=(0.5_8,0._8)*(AsyHank1(Dble(p),x)+AsyHank2(Dble(p),x))
    EndIf
Return
End Function JioComp

!=====================================================
!   Bessel functions of real order and double complex argument
!=====================================================
Double Complex Function JroComp(p,x)
    Implicit None
    Integer:: k
    Double Precision:: p, p2
    Double Complex:: Ck, Tmp, X, xp2
    Double Complex:: cCk, cTmp, cHio
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=0.5_8*x
        Ck=1._8/Gamma(p+1._8)
        JroComp = Ck
        Do k=1,35
            Ck = -Ck*xp2*xp2/(k*(k+p))
            JroComp = JroComp+Ck
        EndDo
        JroComp=Jrocomp*(xp2**p)
    Else
        JroComp=0.5_8*(AsyHank1(p,x)+AsyHank2(p,x))
    EndIf
Return
End Function JroComp

!=====================================================
!   Newman functions of integer order and double precision argument
!=====================================================
Double Precision Function NioReal(p,x)
    Implicit none
    Integer:: p, k,p2
    Double Precision:: x, xp2, xp2p, Ak, Bk, Dk, Den
    Double Precision:: Tmp=0._8
    Double Complex:: cCk, cTmp, cHio
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
        Double Precision Function DigammaInt(N)
            Integer:: N
        End function DigammaInt
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=0.5_8*x
        Ak=1._8/Gamma(p+1._8)
        p2=p+1
        Tmp=Digamma(p2)
        Bk=2._8*Log(xp2)+0.577215664901532860606512_8-Tmp
        NioReal = Ak*Bk
        Do k=1, 35
            Den=1._8/(k*(k+p))
            Ak=-Ak*xp2*xp2*Den
            Bk=Bk-(p+k+k)*Den
            NioReal=NioReal+Ak*Bk
        EndDo
        xp2p=xp2**p
        NioReal=NioReal*xp2p
        Tmp=0._8
        If (p.GE.1) Then
            Dk=Gamma(Dble(p))
            Tmp=Dk
            Do k=1, p-1
                Dk=Dk*xp2*xp2/(k*(p-k))
                Tmp=Tmp+Dk
            EndDo
            Tmp=Tmp/xp2p
        EndIf
        NioReal=0.31830988618379067153776752674503_8*(NioReal-Tmp)
    Else
        Tmp=Dble(p)
        NioReal=Aimag(AsyHank1(Tmp,Dcmplx(x,0._8)))
    EndIf
End Function NioReal

!=====================================================
!   Newman functions of integer order and double precision argument
!=====================================================
Double Precision Function NroReal(p,x)
    Implicit none
    Double Precision:: x, p
    Interface
        Double Precision Function JioReal(N, Y)
            Integer:: N
            Double Precision:: Y
        End function JioReal
        Double Precision Function NioReal(N,Y)
            Integer:: N
            Double Precision:: Y
        End function NioReal
    End Interface
    If (Abs(p-Anint(p)).LE.0.0000005_8) Then
        NroReal=BesN(Int(Anint(p)), x)
    Else
        Write(*,*) 'Not worked!'
        NroReal=0._8
    EndIf
End Function NroReal

!=====================================================
!   Newman functions of integer order and double complex argument
!=====================================================
Double Complex Function NioComp(p,x)
    Implicit none
    Integer:: p, k,p2
    Double Complex:: x, xp2, xp2p, Ak, Bk, Dk, Den
    Double Precision:: Tmp=0._8
    Double Complex:: cCk, cTmp, cHio
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
        Double Precision Function DigammaInt(N)
            Integer:: N
        End function DigammaInt
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=(0.5_8,0._8)*x
        Ak=(1._8,0._8)/Gamma(p+1._8)
        p2=p+1
        Tmp=Digamma(p2)
        Bk=2._8*Log(xp2)+(0.577215664901532860606512_8,0._8)-Dcmplx(Tmp,0._8)
        NioComp = Ak*Bk
        Do k=1, 35
            Den=(1._8,0._8)/(k*(k+p))
            Ak=-Ak*xp2*xp2*Den
            Bk=Bk-(p+k+k)*Den
            NioComp=NioComp+Ak*Bk
        EndDo
        xp2p=xp2**p
        NioComp=NioComp*xp2p
        Tmp=0._8
        If (p.GE.1) Then
            Dk=Gamma(Dble(p))
            Tmp=Dk
            Do k=1, p-1
                Dk=Dk*xp2*xp2/(k*(p-k))
                Tmp=Tmp+Dk
            EndDo
            Tmp=Tmp/xp2p
        EndIf
        NioComp=0.31830988618379067153776752674503_8*(NioComp-Tmp)
    Else
        Tmp=Dble(p)
        NioComp=(0._8,-0.5_8)*(AsyHank1(Tmp,x)-AsyHank2(Tmp,x))
    EndIf
End Function NioComp

!=====================================================
!   Newman functions of integer order and double precision argument
!=====================================================
Double Complex Function NroComp(p,x)
    Implicit none
    Double Complex:: x
    Double Precision:: Tmp, p, p2
    Interface
        Double Precision Function JioComp(N, Y)
            Double Precision:: N
            Double Complex:: Y
        End function JioComp
        Double Precision Function NioComp(N,Y)
            Double Precision:: N
            Double Complex:: Y
        End function NioComp
    End Interface
    If (Abs(p-Anint(p)).LE.0.0000005_8) Then
        NroComp=BesN(Int(Anint(p)), x)
    Else
        Write(*,*) 'Not worked!'
        NroComp=(0._8,0._8)
    EndIf
End Function NroComp

!=====================================================
!   Hankel functions first kind, integer order and double precision argument
!=====================================================
Double Complex Function H1ioReal(p,x)
    Implicit none
    Integer:: p, k, p2
    Double Precision:: x, xp2, xp2p, Tmp=0._8
    Double Complex:: cCk, cTmp, cHio, Ak, Bk, Dk, Den
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
        Double Precision Function DigammaInt(N)
            Integer:: N
        End function DigammaInt
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=(0.5_8,0._8)*x
        Ak=(1._8,0._8)/Gamma(p+1._8)
        p2=p+1
        Tmp=Digamma(p2)
        Bk=(2._8,0._8)*Log(xp2)+(0.577215664901532860606512_8,0._8)-Dcmplx(Tmp,0._8)-(0._8,3.1415926535897932384626433832795_8)
        H1ioReal = Ak*Bk
        Do k=1, 35
            Den=(1._8,0._8)/(k*(k+p))
            Ak=-Ak*xp2*xp2*Den
            Bk=Bk-(p+k+k)*Den
            H1ioReal=H1ioReal+Ak*Bk
        EndDo
        xp2p=xp2**p
        H1ioReal=H1ioReal*xp2p
        cTmp=(0._8,0._8)
        If (p.GE.1) Then
            Dk=Dcmplx(Gamma(Dble(p)),0._8)
            cTmp=Dk
            Do k=1, p-1
                Dk=Dk*xp2*xp2/(k*(p-k))
                cTmp=cTmp+Dk
            EndDo
            cTmp=cTmp/xp2p
        EndIf
        H1ioReal=(0._8,0.31830988618379067153776752674503_8)*(H1ioReal-cTmp)
    Else
        H1ioReal=AsyHank1(Dble(p),Dcmplx(x,0._8))
    EndIf
End Function H1ioReal

!=====================================================
!   Hankel functions first kind, integer order and double complex argument
!=====================================================
Double Complex Function H1ioComp(p, x)
    Implicit none
    Integer:: p, k, p2
    Double Precision:: Tmp=0._8
    Double Complex:: cCk, cTmp, cHio, Ak, Bk, Dk, Den, x, xp2, xp2p
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
        Double Precision Function DigammaInt(N)
            Integer:: N
        End function DigammaInt
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=(0.5_8,0._8)*x
        Ak=(1._8,0._8)/Gamma(p+1._8)
        p2=p+1
        Tmp=Digamma(p2)
        Bk=(2._8,0._8)*Log(xp2)+(0.577215664901532860606512_8,0._8)-Dcmplx(Tmp,0._8)-(0._8,3.1415926535897932384626433832795_8)
        H1ioComp = Ak*Bk
        Do k=1, 35
            Den=(1._8,0._8)/(k*(k+p))
            Ak=-Ak*xp2*xp2*Den
            Bk=Bk-(p+k+k)*Den
            H1ioComp=H1ioComp+Ak*Bk
        EndDo
        xp2p=xp2**p
        H1ioComp=H1ioComp*xp2p
        cTmp=(0._8,0._8)
        If (p.GE.1) Then
            Dk=Dcmplx(Gamma(Dble(p)),0._8)
            cTmp=Dk
            Do k=1, p-1
                Dk=Dk*xp2*xp2/(k*(p-k))
                cTmp=cTmp+Dk
            EndDo
            cTmp=cTmp/xp2p
        EndIf
        H1ioComp=(0._8,0.31830988618379067153776752674503_8)*(H1ioComp-cTmp)
    Else
        Tmp=Dble(p)
        H1ioComp=AsyHank1(Tmp,x)
    EndIf
End Function H1ioComp

!=====================================================
!   Hankel functions second kind, integer order and double precision argument
!=====================================================
Double Complex Function H2ioReal(p,x)
    Implicit none
    Integer:: p
    Double Precision:: x
    Double Complex:: cTmp
    Interface
        Double Complex Function H1ioReal(O, Y)
            Integer:: O
            Double Precision:: Y
        End function H1ioReal
    End Interface
    cTmp=BesH1(p, x)
    H2ioReal=Conjg(cTmp)
End Function H2ioReal

!=====================================================
!   Hankel functions second kind, integer order and double complexn argument
!=====================================================
Double Complex Function H2ioComp(p,x)
    Implicit none
    Integer:: p, k, p2
    Double Precision:: Tmp=0._8
    Double Complex:: cCk, cTmp, cHio, Ak, Bk, Dk, Den, x, xp2, xp2p
    Interface
        Double Precision Function GammaReal(Y)
            Double Precision:: Y
        End function GammaReal
        Double Precision Function DigammaInt(N)
            Integer:: N
        End function DigammaInt
    End Interface
    If (Abs(x).LE.20._8) Then
        xp2=(0.5_8,0._8)*x
        Ak=(1._8,0._8)/Gamma(p+1._8)
        p2=p+1
        Tmp=Digamma(p2)
        Bk=(2._8,0._8)*Log(xp2)+(0.577215664901532860606512_8,0._8)-Dcmplx(Tmp,0._8)+(0._8,3.1415926535897932384626433832795_8)
        H2ioComp = Ak*Bk
        Do k=1, 35
            Den=(1._8,0._8)/(k*(k+p))
            Ak=-Ak*xp2*xp2*Den
            Bk=Bk-(p+k+k)*Den
            H2ioComp=H2ioComp+Ak*Bk
        EndDo
        xp2p=xp2**p
        H2ioComp=H2ioComp*xp2p
        cTmp=(0._8,0._8)
        If (p.GE.1) Then
            Dk=Dcmplx(Gamma(Dble(p)),0._8)
            cTmp=Dk
            Do k=1, p-1
                Dk=Dk*xp2*xp2/(k*(p-k))
                cTmp=cTmp+Dk
            EndDo
            cTmp=cTmp/xp2p
        EndIf
        H2ioComp=-(0._8,0.31830988618379067153776752674503_8)*(H2ioComp-cTmp)
    Else
        Tmp=Dble(p)
        H2ioComp=AsyHank2(Tmp,x)
    EndIf
End Function H2ioComp

!=====================================================
!   Hankel functions of first kind asymptotics
!=====================================================
Double Complex Function AsyHank1(p,x)
    Implicit None
    Integer:: k
    Double Precision:: p, p2, Tmp=0._8
    Double Complex:: cCk, cTmp, cHio, x
    If (Aimag(x).LT.0._8) Then
        AsyHank1=conjg(AsyHank2(p,x))
    Else
        cCk=(1._8,0._8)
        cHio=cCk
        p2=p*p
        Do k=1,60
            cTmp=Dcmplx(k-0.5_8,0._8)
            cCk=cCk*(0._8,0.5_8)*(p2-cTmp*cTmp)/(k*x)
            cHio=cHio+cCk
        EndDo
        cTmp=Exp((0._8,1._8)*(x-0.78539816339744830961566084581988_8-1.5707963267948966192313216916398_8*p))
        cTmp=cTmp*sqrt(0.63661977236758134307553505349006_8/x)
        AsyHank1=cTmp*cHio
    EndIf
End Function AsyHank1

!=====================================================
!   Hankel functions of second kind asymptotics
!=====================================================
Double Complex Function AsyHank2(p,x)
    Implicit None
    Integer:: k
    Double Precision:: p, p2, Tmp=0._8
    Double Complex:: cCk, cTmp, cHio, x
    If (Aimag(x).GE.0._8) Then
        AsyHank2=conjg(AsyHank1(p,x))
    Else
        cCk=(1._8,0._8)
        cHio=cCk
        p2=p*p
        Do k=1,60
            cTmp=Dcmplx(k-0.5_8,0._8)
            cCk=cCk*(0._8,-0.5_8)*(p2-cTmp*cTmp)/(k*x)
            cHio=cHio+cCk
        EndDo
        cTmp=Exp((0._8,-1._8)*(x-0.78539816339744830961566084581988_8-1.5707963267948966192313216916398_8*p))
        cTmp=cTmp*sqrt(0.63661977236758134307553505349006_8/x)
        AsyHank2=cTmp*cHio
    EndIf
End Function AsyHank2

End Module Specfun
Module Series

!=====================================================
!   Interface for subprograms which make series summation by
!   Euler transform
!   SerSumEuler(Summa, Cterm, No)
!   where Summa, Cterm - is the double precision or double complex
!   variables (summa and current term), No - integer variable (current term).
!=====================================================
Interface SerSumEuler
    Module procedure SerSumEulReal, SerSumEulComp
End Interface SerSumEuler


Contains
!====================================================================
!Subroutine of Eiler Series Transformation for quick Summ Computing (double precision terms)
!====================================================================
Subroutine SerSumEulReal(Sum,KTerm,kNum)
Implicit None
Double Precision, Save:: Coef
Double Precision, Dimension(0:150), Save:: SeriesTerm
Double Precision:: DS, Psum, Sum, KTerm, AJ
Integer:: kNum, J
    If (kNum.EQ.0) Then
        SeriesTerm=0._8
        SeriesTerm(0)=KTerm
        Coef=0.5_8
        DS=1._8
        PSum=SeriesTerm(0)
        Sum=0._8
    Else
        SeriesTerm(kNum)=KTerm
        Coef=Coef*0.5_8
        DS=1._8
        PSum=SeriesTerm(kNum)
        Do J=1, kNum
            DS=DS*Dble(kNum-J+1)/Dble(J)
            PSum=PSum+DS*SeriesTerm(kNum-J)
        EndDo
    EndIf
    Sum=Sum+Psum*Coef
    Return
End Subroutine SerSumEulReal

!====================================================================
!Subroutine of Eiler Series Transformation for quick Summ Computing (double complex terms)
!====================================================================
Subroutine SerSumEulComp(Sum,KTerm,kNum)
Implicit None
Double Complex, Save:: Coef
Double Complex, Dimension(0:150), Save:: SeriesTerm
Double Complex:: DS, Psum, Sum, KTerm, AJ
Integer:: kNum, J
    If (kNum.EQ.0) Then
        SeriesTerm=(0._8,0._8)
        SeriesTerm(0)=KTerm
        Coef=(0.5_8,0._8)
        DS=(1._8,0._8)
        PSum=SeriesTerm(0)
        Sum=(0._8,0._8)
    Else
        SeriesTerm(kNum)=KTerm
        Coef=Coef*(0.5_8,0._8)
        DS=(1._8,0._8)
        PSum=SeriesTerm(kNum)
        Do J=1, kNum
            DS=DS*Dcmplx(Dble(kNum-J+1)/Dble(J),0._8)
            PSum=PSum+DS*SeriesTerm(kNum-J)
        EndDo
    EndIf
    Sum=Sum+Psum*Coef
    Return
End Subroutine SerSumEulComp

End Module Series
Module Service

!=====================================================
!   Interface for subprograms which gives the number of digits in integer number
!   K=NumOfDig(I)
!   where I, K - is the integer variables.
!=====================================================
Interface NumOfDig
    Module procedure NumOfDig
End Interface NumOfDig
!=====================================================
!   Interface for subprograms which convert an integer variable to string
!   String=NumToString(I[, len])
!   where I - is the operang (integer variable), len - is the output string length 
!   (intrger, optional parameter). If len is specified, result string would be appended
!   by left zeros to len symbols. If len is not specified, String has len=10 and no zeros.
!=====================================================
Interface NumToString
    Module procedure NumToString, NumToStringZero
End Interface NumToString
!=====================================================
!   Interface for subprograms which convert a digits string to double precision
!   X=StringToReal(String)
!   where String - is the string consists of digits and (maybe) decimal dot. The 
!   length of String may be arbitrary. X - is the double precision result of conversion.
!=====================================================
Interface StringToReal
    Module procedure StringToReal
End Interface StringToReal
!=====================================================
!   Interface for subprograms which convert a string of time (from DATE_AND_TIME
!   subroutine, generally) to double precision number
!   X=TimeToReal(String)
!   where String - is the string of time such as '145515.785' (length = 10!). X -
!   is the double precision result of conversion.
!=====================================================
Interface TimeToReal
    Module procedure TimeToReal
End Interface TimeToReal
!=====================================================
!   Interface for subprograms which compute the time interval into the day for
!   two string times.
!   X=TimeInterval(Start, Finish)
!   where Start and Finish - is the strings of time such as '145515.785' (length = 10!)
!   for start and finish of any process. If start>finish, the result should be increased
!   to 24 hour. X - is the double precision result of conversion.
!=====================================================
Interface TimeInterval
    Module procedure TimeInterval
End Interface TimeInterval
!=====================================================
!   Interface for subprograms which compute the polar angle for point in Cartesian
!   coordinate system
!   Angle=PolarAngle(X, Y)
!   where X, Y - is the Cartesian coordinates (double precision variables). Angle is
!   polar angle in radians (double precision), counted out X axis. 0<=Angle<=2 pi
!=====================================================
Interface PolarAngle
    Module Procedure PolarAngle
End Interface PolarAngle
!=====================================================
!   Interface for subprograms which compute the polar angle for point in Cartesian
!   coordinate system
!   Angle=PolarAngleSym(X, Y)
!   where X, Y - is the Cartesian coordinates (double precision variables). Angle is
!   polar angle in radians (double precision), counted out X axis. -pi <=Angle<= pi
!=====================================================
Interface PolarAngleSym
    Module Procedure PolarAngleSym
End Interface PolarAngleSym


!=====================================================
!   Module procedures.
!=====================================================
Contains

!=====================================================
!   The number of digits
!=====================================================
Integer Function NumOfDig(Inum)
    Implicit None
    Integer:: Curr, Inum, J
    J=1
    Curr=Inum
    Do while (Curr-Mod(Curr, 10).NE.0)
        Curr=(Curr-Mod(Curr, 10))/10
        J=J+1
    EndDo
    NumOfDig=J
End Function NumOfDig

!=====================================================
!   The number to string conversion
!=====================================================
Function NumToString(Inum)
    Implicit None
    Integer:: Alen, Curr, Tmp, Next, Inum, J, K
    Character(Len=10):: NumToString
    NumToString=' '
    J=NumOfDig(Inum)
    Curr=Inum
    Do K=1, J
        Tmp=Curr
        Curr=Mod(Curr,10)
        Next=(Tmp-Curr)/10
        NumToString=Achar(Curr+48)//Trim(NumToString)
        Curr=Next
    EndDo
End Function NumToString

!=====================================================
!   The number to string conversion with zeros appending
!=====================================================
Function NumToStringZero(Inum, Digs)
    Implicit None
    Integer:: Alen, Curr, Tmp, Next, Inum, J, K, Digs
    Character(Digs):: NumToStringZero
    NumToStringZero=' '
    J=NumOfDig(Inum)
    Curr=Inum
    Do K=1, J
        Tmp=Curr
        Curr=Mod(Curr,10)
        Next=(Tmp-Curr)/10
        NumToStringZero=Achar(Curr+48)//Trim(NumToStringZero)
        Curr=Next
    EndDo
    Do K=1, Digs-J
        NumToStringZero='0'//Trim(NumToStringZero)
    EndDo
End Function NumToStringZero

!=====================================================
!   The string to real conversion
!=====================================================
Double Precision Function StringToReal(InString)
    Implicit none
    Character (len=*):: InString
    Integer:: Pdot, Slen, I, J, K
    Double Precision:: Tmp
    Pdot=Index(InString,'.')
    !Write(*,*) InString, Pdot
    Slen=Len_trim(InString)
    If (Pdot.NE.0) Then
        J=0
        Do K=Slen, Pdot+1, -1
            I=Iachar(InString(K:K))-48
            J=J+I*10**(Slen-K)
        EndDo
        Tmp=J*(0.1_8**(Slen-Pdot))
    Else
        Tmp=0._8
        Pdot=Slen+1
    EndIf
    J=0
    Do K=1,Pdot-1
        I=Iachar(InString(K:K))-48
        J=J+I*10**(Pdot-1-K)
    EndDo
    StringToReal=Dble(J)+Tmp
End function StringToReal

!=====================================================
!   Time string to seconds from midnigth (in double precision) conversion
!=====================================================
Double Precision Function TimeToReal(InDate)
    Implicit none
    Character (len=10) Indate
    TimeToReal=(StringToReal(InDate(1:2))*60._8+StringToReal(InDate(3:4)))*60._8+StringToReal(InDate(5:10))
End function TimeToReal

!=====================================================
!   Time interval computing
!=====================================================
Double Precision Function TimeInterval(Starttime, Finishtime)
    Implicit none
    Character (len=10) Finishtime, Starttime
    Double Precision:: Ftime, Stime, Tmp
    Ftime=TimeToReal(Finishtime)
    Stime=TimeToReal(Starttime)
    If (Ftime.LT.Stime) Then
        Tmp=TimeToReal('240000.000')
    Else
        Tmp=0._8
    EndIf
    TimeInterval=Tmp+Ftime-Stime
End function TimeInterval

!=====================================================
!   The polar angle (from zero to 2 pi) computing
!=====================================================
Double Precision Function PolarAngle(X, Y)
    Implicit None
    Double Precision:: X, Y, Phiw
        If (Y.GE.0._8) Then
            Phiw = PolarAngleSym (X,Y)
        Else
            Phiw = 6.2831853071795862319959269370883703232_8-Acos(X/Sqrt(X*X+Y*Y))
        EndIf
        PolarAngle=Phiw
End Function PolarAngle

!=====================================================
!   The polar angle (from -pi to pi) computing
!=====================================================
Double Precision Function PolarAngleSym(X, Y)
    Implicit None
    Double Precision:: X, Y, Phiw
    If ((Abs(X).LE.Epsilon(X)).AND.(Abs(Y).LE.Epsilon(Y))) Then
        Phiw=0._8
    Else
        Phiw= Sign(Acos(X/Sqrt(X*X+Y*Y)), Y)
    EndIf   
    PolarAngleSym = Phiw
End Function PolarAngleSym
End Module Service