subroutine print_debug()
    use variables 
    implicit none
!    integer(8) :: i,j,k,naux
!    
!      do i=1,np
!         do j=1,nv1(i)
!             print *, "i,nv1,vv1(i,j) = ",i,nv1(i), vv1(i,j)
!         end do            
!       end do

!       do i=1,np
!          do j=1,nv(i)
!              print *, "i,nv,vv(i,j) = ",i,nv(i), vv(i,j)
!          end do            
!        end do


!       do i=1,np
!         do j=1,nv_t(i)
!             print *, "i,nv_t,vv_t(i,j) = ",i,nv_t(i), vv_t(i,j)
!         end do            
!       end do
! 
!       do i=1,np
!         naux=0
!         do j=1,nv_t(i)
!             do k=j+1,nv_t(i)
!                 naux=naux+1
!                 print *, "i,j,k,nv_t,costhetaijk0(i,naux) = ",i,vv_t(i,j),vv_t(i,k),nv_t(i),costhetaijk0(i,naux)
!             end do            
!         end do
!       end do

      
!      do i=1,np
!       viz(i)%vv1(1)=i
!       print *, "viz(i)%vv1(1) = ", viz(i)%vv1(1)
!      end do
! !    i=1
! !         print *, "pos(i,1:3) = ", pos(i,1),pos(i,2),pos(i,3)
    
    return 
end

subroutine rk1()
    use variables 
    implicit none
    real(kind=8):: tempo0

    tempo0=tempo
    pos0=pos
    vel0=vel
    
    call define_ace()
    k1v=vel
    k1a=ace

    pos=pos0+k1v*dt
    vel=vel0+k1a*dt
    tempo=tempo0+dt    
    
    return 
end

subroutine rk2()
    use variables 
    implicit none
    real(kind=8):: tempo0

    tempo0=tempo
    pos0=pos
    vel0=vel
    
    call define_ace()
    k1v=vel
    k1a=ace

    pos=pos0+k1v*dt2
    vel=vel0+k1a*dt2
    tempo=tempo0+dt2    
    
    call define_ace()
    k2v=vel
    k2a=ace

    pos=pos0+(k1v+k2v)*dt2
    vel=vel0+(k1a+k2a)*dt2
    tempo=tempo0+dt    
    
    return 
end

subroutine rk4()
    use variables 
    implicit none
    real(kind=8):: tempo0

    tempo0=tempo
    pos0=pos
    vel0=vel
    
    call define_ace()
    k1v=vel
    k1a=ace

    pos=pos0+k1v*dt2
    vel=vel0+k1a*dt2
    tempo=tempo0+dt2    
    
    call define_ace()
    k2v=vel
    k2a=ace

    pos=pos0+k2v*dt2
    vel=vel0+k2a*dt2
    tempo=tempo0+dt2    
    
    call define_ace()
    k3v=vel
    k3a=ace

    pos=pos0+k3v*dt
    vel=vel0+k3a*dt
    tempo=tempo0+dt    
    
    call define_ace()
    k4v=vel
    k4a=ace
    
    pos=pos0+(k1v+2.d0*k2v+2.d0*k3v+k4v)*dt6
    vel=vel0+(k1a+2.d0*k2a+2.d0*k3a+k4a)*dt6
    tempo=tempo0+dt    
    

    call define_ace()
 
    pkn3=pkn2 ; pkn2=pkn1 ; pkn1=pkn ; pkn=vel
    vkn3=vkn2 ; vkn2=vkn1 ; vkn1=vkn ; vkn=ace    
    
    return 
end


subroutine PreditorCorretor()
    use variables 
    implicit none
    real(kind=8):: tempo0

    tempo0=tempo
    call define_ace()

    pos0=pos ; vel0=vel
    vkn=ace  ; pkn=vel 
!   
    pos=pos0+(55*pkn-59*pkn1+37*pkn2-9*pkn3)*dt24
    vel=vel0+(55*vkn-59*vkn1+37*vkn2-9*vkn3)*dt24
!
    tempo=tempo0+dt
    call define_ace()
!
    pos=pos0+(9*vel+19*pkn-5*pkn1+pkn2)*dt24
    vel=vel0+(9*ace+19*vkn-5*vkn1+vkn2)*dt24
!
    pkn3=pkn2 ; pkn2=pkn1 ; pkn1=pkn
    vkn3=vkn2 ; vkn2=vkn1 ; vkn1=vkn
!
    tempo=tempo0+dt

    return 
end



subroutine define_ace()
 use variables 
 implicit none
 integer(kind=8) :: i,j,k,jv,kv,naux,SomaTipo!,i1,indexi,indexj
 real(kind=8) :: xi,yi,zi,V11,V12,F11,F12,V0,Gx
 real(kind=8) :: xij,yij,zij,rij2,rij,rij1,forcamola,forcalj,sigrij2,sigrij6,sigrij8,sig,epi
 real(kind=8) :: xik,yik,zik,rik2,rik,rik1,costhetaijk,rijdotrij,costhetaijkMthetaijk0
 real(kind=8) :: rijr02,forcaFENE,auxKfene,auxRfene,auxDRfene,auxmas,yiw1,yiw2,riw12,riw22
! !
    V0=1.0
    Epot=0.0D0
    Gx=d2pi*0.5D0

    ace(:,1)=0.D0
    ace(:,2)=0.D0
    ace(:,3)=0.D0

    forca(:,1)=0.0D0
    forca(:,2)=0.0D0
    forca(:,3)=0.0D0

    do i=1,np

        xi=pos(i,1) ; yi=pos(i,2) ; zi=pos(i,3)
        yiw1=yi-wall1
        yiw2=yi-wall2
        
        riw12=1.D0/yiw1**2
        riw22=1.D0/yiw2**2
        
        V11=V0*riw12**3
        V12=V0*riw22**3
        
        F11=V0*(6*yiw1*riw12**4)
        F12=V0*(6*yiw2*riw22**4)
        
        forca(i,1)=forca(i,1)        +V12*dsin(Gx*xi)*Gx
        forca(i,2)=forca(i,2)+F11+F12+F12*dcos(Gx*xi)          
        ! 

!  loop LENNARD JONES 
        do jv=1,nv(i)
!              print *, "i,nv,vv(i,j) = ",i,nv(i), vv(i,j)
            j=vv(i,jv)
            xij=pos(j,1)-xi ; yij=pos(j,2)-yi ; zij=pos(j,3)-zi
            rij2=xij**2 + yij**2 + zij**2

!           rij=DSQRT(rij2)
            SomaTipo=tipo(i)+tipo(j)
            sig=sigij(SomaTipo)
            epi=epiij(SomaTipo)
!     
            if ( tipo(i).ne.tipo(j) ) then
                sigrij2=sig*sig/rij2
                sigrij6=sigrij2*sigrij2*sigrij2
                sigrij8=sigrij6*sigrij2
                Epot=Epot+epi*sigrij6*(sigrij6-2)
                forcalj=12*epi*sigrij8*(1-sigrij6)/(sig*sig)
            else
                if (rij2 <= sig) then
                    sigrij2=sig*sig/rij2
                    sigrij6=sigrij2*sigrij2*sigrij2
                    sigrij8=sigrij6*sigrij2
                    Epot=Epot+epi*sigrij6*(sigrij6-2)+epi
                    forcalj=12*epi*sigrij8*(1-sigrij6)/(sig*sig)
                else
                    Epot=0.D0
                    forcalj=0.D0
                end if
            end if

            forca(i,1)=forca(i,1)+forcalj*xij
            forca(i,2)=forca(i,2)+forcalj*yij
            forca(i,3)=forca(i,3)+forcalj*zij

            forca(j,1)=forca(j,1)-forcalj*xij
            forca(j,2)=forca(j,2)-forcalj*yij
            forca(j,3)=forca(j,3)-forcalj*zij

        end do 
!  
 !  loop FENE
        do jv=1,nv1(i)
            j=vv1(i,jv)
            xij=pos(j,1)-xi ; yij=pos(j,2)-yi ; zij=pos(j,3)-zi
            rij2=xij**2 + yij**2 + zij**2
            rij=DSQRT(rij2)
!           
            SomaTipo=tipo(i)+tipo(j)
!            print*,i,j,tipo(i),tipo(j),SomaTipo
            auxKfene=Kfene(SomaTipo)
            auxRfene=rfene(SomaTipo)
            auxDRfene=drfene(SomaTipo)

            rijr02=1.D0-((rij-auxRfene)/auxDRfene)**2
   
            Epot=Epot-0.5D0*auxKfene*dlog(rijr02)*auxDRfene**2
            forcaFENE=auxKfene*(1.D0-auxRfene/rij)/rijr02
!            print*,i,j,rij2,tipo(i),tipo(j),auxKfene,auxRfene,auxDRfene
!            print*,rijr02,dlog(rijr02),forcaFENE
!   
            forca(i,1)=forca(i,1)+forcaFENE*xij
            forca(i,2)=forca(i,2)+forcaFENE*yij
            forca(i,3)=forca(i,3)+forcaFENE*zij

            forca(j,1)=forca(j,1)-forcaFENE*xij
            forca(j,2)=forca(j,2)-forcaFENE*yij
            forca(j,3)=forca(j,3)-forcaFENE*zij

    end do 
    
     !  loop bending
     naux=0
     do j=1,nv_t(i)
      jv=vv_t(i,j)
      xij=pos(jv,1)-xi ; yij=pos(jv,2)-yi ; zij=pos(jv,3)-zi
      rij2=xij**2 + yij**2 + zij**2
      rij=DSQRT(rij2)
      rij1=1.D0/rij
!      
      do k=j+1,nv_t(i)
       kv=vv_t(i,k)
       xik=pos(kv,1)-xi ; yik=pos(kv,2)-yi ; zik=pos(kv,3)-zi
       rik2=xik**2 + yik**2 + zik**2
       rik=DSQRT(rik2)
       rik1=1.D0/rik
       rijdotrij=xij*xik+yij*yik+zij*zik
       costhetaijk=rijdotrij*rij1*rik1
       naux=naux+1
       costhetaijkMthetaijk0=costhetaijk-costhetaijk0(i,naux)
!       
       SomaTipo=tipo(i)+tipo(jv)
!       print*,i,jv,kv,tipo(i),tipo(jv),tipo(kv),SomaTipo,0.5D0*kmolaAng(SomaTipo)*costhetaijkMthetaijk0**2
       Epot=Epot+0.5D0*kmolaAng(SomaTipo)*costhetaijkMthetaijk0**2
       forcamola=kmolaAng(SomaTipo)*costhetaijkMthetaijk0*rij1*rik1
! 
       forca(i,1)=forca(i,1)+forcamola*(xij+xik-rijdotrij*(rij1*rij1*xij+rik1*rik1*xik) )
       forca(i,2)=forca(i,2)+forcamola*(yij+yik-rijdotrij*(rij1*rij1*yij+rik1*rik1*yik) )
       forca(i,3)=forca(i,3)+forcamola*(zij+zik-rijdotrij*(rij1*rij1*zij+rik1*rik1*zik) )
! 
       forca(jv,1)=forca(jv,1)-forcamola*(xik-rijdotrij*rij1*rij1*xij)
       forca(jv,2)=forca(jv,2)-forcamola*(yik-rijdotrij*rij1*rij1*yij)
       forca(jv,3)=forca(jv,3)-forcamola*(zik-rijdotrij*rij1*rij1*zij)

       forca(kv,1)=forca(kv,1)-forcamola*(xij-rijdotrij*rik1*rik1*xik)
       forca(kv,2)=forca(kv,2)-forcamola*(yij-rijdotrij*rik1*rik1*yik)
       forca(kv,3)=forca(kv,3)-forcamola*(zij-rijdotrij*rik1*rik1*zik)
      end do
     end do 
!
    end do
!  
!  
! 
!  
! !    if (tempo.lt.tpuxar) then
! !     x0mola=atomo(iext)%pos%x
! !     forcaMolaexterna_x=0.0D0
! !    end if
! !    if (tempo.gt.tpuxar) then
! !     gamm=0.0D0
! !     forcaMolaexterna_x=-kmolaext*(atomo(iext)%pos%x-x0mola)
! !     atomo(iext)%forca%x=atomo(iext)%forca%x+forcaMolaexterna_x
! !     x0mola=x0mola+0.5*dt
! !    end if
!  
 Fxcm=0.D0 ; Fycm=0.D0 ; Fzcm=0.D0
 do i=1,Np
  auxmas=1.D0/massa(i)
  ace(i,1)=forca(i,1)*auxmas!-gamm*vel(i,1)
  ace(i,2)=forca(i,2)*auxmas!-gamm*vel(i,2)
  ace(i,3)=forca(i,3)*auxmas!-gamm*vel(i,3)
 
  Fxcm=Fxcm+forca(i,1)*massa(i)
  Fycm=Fycm+forca(i,2)*massa(i)
  Fzcm=Fzcm+forca(i,3)*massa(i)
 end do

  Fxcm=Fxcm/MassT
  Fycm=Fycm/MassT
  Fzcm=Fzcm/MassT
 
! ! 
!  do i=1,Np
!   ace(i,1)=ace(i,1)-gamm*vel(i,1)
!   ace(i,2)=ace(i,2)-gamm*vel(i,2)
!   ace(i,3)=ace(i,3)-gamm*vel(i,3)
!  end do
! !
    
    return 
end


subroutine enercin()
    use variables 
    implicit none
    integer(8) :: i
    real(kind=8) :: vxi,vyi,vzi
    
    Vxcm=0.D0 ; Vycm=0.D0 ; Vzcm=0.D0
    Ecin=0.0D0
    do i=1,np
        vxi = vel(i,1)
        vyi = vel(i,2)
        vzi = vel(i,3)
        Ecin = Ecin+massa(i)*(vxi*vxi+vyi*vyi+vzi*vzi)*0.5D0
        Vxcm = Vxcm + vxi*massa(i)
        Vycm = Vycm + vyi*massa(i)
        Vzcm = Vzcm + vzi*massa(i)
    end do
!    print*,Ecin
     Vxcm=Vxcm/MassT
     Vycm=Vycm/MassT
     Vzcm=Vzcm/MassT 
    
    return 
end


subroutine renormvel()
    use variables 
    implicit none
    integer(8) :: i
    real(8) :: alpha,Tt
    
   
    Ecin=0.0D0
    do i=1,np
        Ecin=Ecin+massa(i)*(vel(i,1)**2+vel(i,2)**2+vel(i,3)**2)*0.5D0
    end do
    Tt=(2.D0*Ecin/(3.D0*np))
    
    alpha=sqrt(Temp/Tt)
    
    vel=vel*alpha
    
!     do i=1,np
!         vel(i,1)=vel(i,1)*alpha
!         vel(i,2)=vel(i,2)*alpha
!         vel(i,3)=vel(i,3)*alpha
!     end do
!    print*,Ecin,alpha
    
    return 
end

! 
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C SUB-ROTINA DEFINE TODAS AS VARIAVEIS DE ENTRADA
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine PotStelle(X_p,Y_p,Z_p,i1,j1)
    use variables
    implicit none
    real(kind=8), intent(in) :: X_p,Y_p,Z_p
    integer(kind=8), intent(in) :: i1,j1
    real(kind=8) :: sigma1,epsilon1

    integer(kind=8) :: ii,m1 

    
    real(kind=8) :: x_l,y_l,z_l,a1,r0,d76,dz_l
    real(kind=8) :: as,q1,as_l,A3,A6,A630,A11,A12,gn_l
    real(kind=8) :: dcox,dcoy,dsix,dsiy,dcoxy
    real(kind=8) :: t20,t30,t40,t50,t11,t13,t14,t41,t31,t42,v11,v12,v13
    real(kind=8) :: V0z,V0z_l,v14,v21,v22,v23,Vz,Vz_l
    real(kind=8) :: bes0,bes1,bes2,bes3,bes4,bes5
    real(kind=8) :: BESSK,BESSK1,BESSK0
    
    ! 
    d76=2.0D0**(7.0D0/6.0D0)
    r0=d76/DSQRT(3.0D0)

    sigma1=sigij(i1+j1)
    epsilon1=epiij(i1+j1)
    a1=r0
    dz_l=1.0D0/(2.0D0*a1)
    as=a1*a1
    q1=1.0D0
    as_l=as/(a1*a1)
    A3=sigma1/a1
! 
    A6=A3*A3*A3*A3*A3*A3
    A630=A6/30.0D0
    A11=2.0D0*A6/5.0D0
    A12=d2pi*A6/as_l
    m1=10
    gn_l=d2pi
! 
    x_l=X_p/a1
    y_l=Y_p/a1
    z_l=Z_p/a1
! 
    dcox=dcos(d2pi*x_l)
    dcoy=dcos(d2pi*y_l)
    dsix=dsin(d2pi*x_l)
    dsiy=dsin(d2pi*y_l)
    dcoxy=dcox+dcoy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 


    t30=0.0D0
    t40=0.0D0
    t50=0.0D0
    do ii=0,m1
        t11=z_l+ii*dz_l
        t13=1.0D0/(t11*t11)
        t14=t13*t13
        t20=(A11*t14*t13-1.0D0)*t14
        t30=t30+t20
 
        t41=(A6*t14*t13-1.0D0)*t14/t11
        t40=t40+t41
    end do
! 
! 
! 
! cccccccccc E_0 CCCCCCCCCC
    t31=t30*A12*q1
! cccccccccc E'_0 CCCCCCCCCC
    t42=-4.0D0*t40*A12*q1
! 
    bes0=BESSK0(gn_l*z_l)
    bes1=BESSK1(gn_l*z_l)
    bes2=BESSK(2,gn_l*z_l)
    bes3=BESSK(3,gn_l*z_l)
    bes4=BESSK(4,gn_l*z_l)
    bes5=BESSK(5,gn_l*z_l)
! 
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C
! CCCCCCCCCCCCCCCCC E_n CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    v11=0.5D0*(gn_l/z_l)
    v12=v11*V11
    v13=v12*V12*V11
! 
    v14=A12*(A630*v13*bes5-2.0D0*bes2*v12)
! CCCCCCCCCCCCCCCCC E'_n CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    v21=A630*v13*(10.0D0*bes5/z_l+gn_l*bes4)
    v22=2.0D0*v12*(4.0D0*bes2/z_l+gn_l*bes1)
    v23=-A12*(v21-v22)
! 
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C	Termos que dependem de q_z
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    V0z=epsilon1*t31/a1
    V0z_l=epsilon1*t42/as
! 
    Vz=epsilon1*2.0D0*v14/a1
    Vz_l=epsilon1*2.0D0*v23/as
! 
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C	1_a derivadas
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    Vint=V0z+Vz*dcoxy
! 
    Vint_x=-Vz*dsix*d2pi/a1
    Vint_y=-Vz*dsiy*d2pi/a1
    Vint_z=V0z_l+Vz_l*dcoxy
    
!    write(*,*) Vint

! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
    return
end




      FUNCTION BESSK(N,X)
      IMPLICIT NONE
      INTEGER N,J
      REAL *8 X,BESSK,BESSK0,BESSK1,TOX,BK,BKM,BKP
! ------------------------------------------------------------------------
!     CE SOUS-PROGRAMME CALCULE LA FONCTION BESSEL MODIFIFIEE 3E ESPECE
!     D'ORDRE N ENTIER POUR TOUT X REEL POSITIF > 0.  ON UTILISE ICI LA
!     FORMULE DE RECURRENCE CLASSIQUE EN PARTANT DE BESSK0 ET BESSK1.
!
!     THIS ROUTINE CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD
!     KIND OF INTEGER ORDER, N FOR ANY POSITIVE REAL ARGUMENT, X. THE
!     CLASSICAL RECURSION FORMULA IS USED, STARTING FROM BESSK0 AND BESSK1.
! ------------------------------------------------------------------------ 
!     REFERENCE:
!     C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
!     MATHEMATICAL TABLES, VOL.5, 1962.
! ------------------------------------------------------------------------
      IF (N.EQ.0) THEN
      BESSK = BESSK0(X)
      RETURN
      ENDIF
      IF (N.EQ.1) THEN
      BESSK = BESSK1(X)
      RETURN
      ENDIF
      IF (X.EQ.0.D0) THEN
      BESSK = 1.D30
      RETURN
      ENDIF
      TOX = 2.D0/X
      BK  = BESSK1(X)
      BKM = BESSK0(X)
      DO 11 J=1,N-1
      BKP = BKM+DFLOAT(J)*TOX*BK
      BKM = BK
      BK  = BKP
   11 CONTINUE
      BESSK = BK
      RETURN
      END
! ----------------------------------------------------------------------
      FUNCTION BESSK0(X)
!     CALCUL DE LA FONCTION BESSEL MODIFIEE DU 3EME ESPECE D'ORDRE 0
!     POUR TOUT X REEL NON NUL.
!
!     CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
!     ORDER ZERO FOR ANY POSITIVE REAL ARGUMENT, X.
! ----------------------------------------------------------------------
      IMPLICIT NONE
      REAL*8 X,BESSK0,Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,    &
      BESSI0
      DATA P1,P2,P3,P4,P5,P6,P7/-0.57721566D0,0.42278420D0,0.23069756D0, &
      0.3488590D-1,0.262698D-2,0.10750D-3,0.74D-5/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,-0.7832358D-1,0.2189568D-1, & 
      -0.1062446D-1,0.587872D-2,-0.251540D-2,0.53208D-3/
      IF(X.EQ.0.D0) THEN
      BESSK0=1.D30
      RETURN
      ENDIF
      IF(X.LE.2.D0) THEN
      Y=X*X/4.D0
      AX=-LOG(X/2.D0)*BESSI0(X)
      BESSK0=AX+(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
      Y=(2.D0/X)
      AX=EXP(-X)/DSQRT(X)
      BESSK0=AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      ENDIF
      RETURN
      END
! ----------------------------------------------------------------------
      FUNCTION BESSK1(X)
!     CALCUL DE LA FONCTION BESSEL MODIFIEE DE 3EME ESPECE D'ORDRE 1
!     POUR TOUT X REEL POSITF NON NUL.
!
!     CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
!     ORDER ONE FOR ANY POSITIVE REAL ARGUMENT, X.
! ----------------------------------------------------------------------
      IMPLICIT NONE
      REAL*8 X,BESSK1,Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,BESSI1
      DATA P1,P2,P3,P4,P5,P6,P7/1.D0,0.15443144D0,-0.67278579D0,  &
      -0.18156897D0,-0.1919402D-1,-0.110404D-2,-0.4686D-4/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,0.23498619D0,-0.3655620D-1, &
      0.1504268D-1,-0.780353D-2,0.325614D-2,-0.68245D-3/
      IF(X.EQ.0.D0) THEN
      BESSK1=1.D32
      RETURN
      ENDIF
      IF(X.LE.2.D0) THEN
      Y=X*X/4.D0
      AX=LOG(X/2.D0)*BESSI1(X)
      BESSK1=AX+(1.D0/X)*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
      Y=(2.D0/X)
      AX=EXP(-X)/DSQRT(X)
      BESSK1=AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      ENDIF
      RETURN
      END
!
!     Bessel Function of the 1st kind of order zero.
!
      FUNCTION BESSI0(X)
      IMPLICIT NONE
      REAL *8 X,BESSI0,Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
      DATA P1,P2,P3,P4,P5,P6,P7/1.D0,3.5156229D0,3.0899424D0,1.2067429D0,  &
      0.2659732D0,0.360768D-1,0.45813D-2/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1,  &
      0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,  &
      0.2635537D-1,-0.1647633D-1,0.392377D-2/
      IF(ABS(X).LT.3.75D0) THEN
      Y=(X/3.75D0)**2
      BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
      ELSE
      AX=ABS(X)
      Y=3.75D0/AX
      BX=EXP(AX)/DSQRT(AX)
      AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
      BESSI0=AX*BX
      ENDIF
      RETURN
      END
!
!     Bessel Function of the 1st kind of order one.
!
      FUNCTION BESSI1(X)
      IMPLICIT NONE
      REAL *8 X,BESSI1,Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
      DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,  &
      0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, &
      -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,        &
      -0.2895312D-1,0.1787654D-1,-0.420059D-2/
      IF(ABS(X).LT.3.75D0) THEN
      Y=(X/3.75D0)**2
      BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
      AX=ABS(X)
      Y=3.75D0/AX
      BX=EXP(AX)/DSQRT(AX)
      AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
      BESSI1=AX*BX
      ENDIF
      RETURN
      END

