C     =================================================================
C     File: libmath.f
C     =================================================================

C     =================================================================
C     Module: Auxiliary math subroutines
C     =================================================================

C     Last update of any of the component of this module: 
C
C     Mar 18, 2007.

C     ******************************************************************
C     ******************************************************************

C Calculo de lnGama(x)
C fonte: Numerical Recipes in Fortran
C http://www.phy.bg.ac.yu/ipcf/docs/Numerical_Recipes/bookfpdf.html

      double precision function gammaln(xx)

C     ARGUMENTS SCALARS

      double precision xx

c     LOCAL SCALARS

      integer j
      double precision ser,stp,tmp,x,y

C     LOCAL ARRAYS

      double precision cof(6)

      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     +24.01409824083091d0,-1.231739572450155d0,0.1208650973866179d-2,
     +-0.5395239384953d-5,2.5066282746310005d0/

      x = xx
      y = x
      tmp = x+5.5d0
      tmp = (x + 0.5d0) * log(tmp) - tmp
      ser = 1.000000000190015d0

      do j = 1,6

        y = y + 1.0d0
        ser = ser + cof(j) / y

      end do

      gammaln = tmp + log(stp * ser / x)

      return

      end

C     ******************************************************************
C     ******************************************************************

       double precision function fNormaVetor(N,pVetor)
C
C      Funcao para calcular a norma quadratica de um vetor
C
C      Parametros de Entrada:
C
C      N integer
C        numero de variaveis
C      pVetor double precision pVetor(N)
C             vetor que sera calculada a norma
C

C      SCALAR ARGUMENTS
       integer N

C      ARRAY ARGUMENTS
       double precision pVetor(N)

C      LOCAL SCALARS
       integer i
       double precision vNorma

       vNorma = 0.0d0

       do i=1,N
           vNorma = vNorma + pVetor(i)*pVetor(i)
       end do

       fNormaVetor = sqrt(vNorma)

       return

       end

C     ******************************************************************
C     ******************************************************************

       double precision function fNormaVetorTrunc(N,pVetor,pTrunc)
C
C      Funcao para calcular a norma quadratica de um vetor
C
C      Parametros de Entrada:
C
C      N integer
C        numero de variaveis
C      pVetor double precision pVetor(N)
C             vetor que sera calculada a norma
C
C      pTrunc integer
C            define o numero de pontos flutuantes relevantes

C      SCALAR ARGUMENTS
       integer N,pTrunc

C      ARRAY ARGUMENTS
       double precision pVetor(N)

C      LOCAL SCALARS
       integer i
       double precision vNorma,sNum

       vNorma = 0.0d0

       if( pTrunc .gt. 0 ) then

           do i=1,N
               sNum = aint( pVetor(i) * ( 10.0d0 ** pTrunc ) )
     +                                / ( 10.0d0 ** pTrunc )
           end do

       else

           do i=1,N
               vNorma = vNorma + pVetor(i)*pVetor(i)
           end do

       end if

       fNormaVetorTrunc = sqrt(vNorma)

       return

       end

C     ******************************************************************
C     ******************************************************************

       double precision function fTrunc(pNum,pTrunc)
C
C      Funcao para truncar um numero
C
C      Parametros de Entrada:
C
C      pNum double precision
C             numero que sera truncado
C
C      pTrunc integer
C            define o numero de pontos flutuantes relevantes

C      SCALAR ARGUMENTS
       integer pTrunc
       double precision pNum

       if( pTrunc .gt. 0 ) then

           pNum = aint( pNum * ( 10.0d0 ** pTrunc ) )
     +                         / ( 10.0d0 ** pTrunc )

       end if

       fTrunc = pNum

       return

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para  subtracao de vetores
       subroutine sDiferencaVetor(N,pMatriz1,pMatriz2,pRes)
       integer N,i
       double precision pMatriz1(N),pMatriz2(N),pRes(N)

       do i=1,N
           pRes(i) = pMatriz1(i) - pMatriz2(i)
       end do

       end

C     ******************************************************************
C     ******************************************************************

       subroutine sMultTranspVetor(N,pV1,pV2,pRes)

C      rotina para multiplicacao vetor transposto x vetor
C
C      Parametros de Entrada:
C
C      N integer
C        numero de veriaveis
C      pV1 double precision  pV1(N)
C          vetor transposto
C      pV2 double precision pv2(N)
C          segundo vetor
C
C      Parametros de Saida:
C
C      pRes double precision
C           resultado da multiplicacao
C

C      SCALARS ARGUMENTS
       integer N
       double precision pRes

C      ARRAYS ARGUMENTS
       double precision pV1(N),pV2(N)

C      LOCAL SCALARS
       integer i

       pRes = 0.0d0

       do i=1,N
         pRes=pRes+(pV1(i)*pV2(i))
       end do

       end