module util
 
 use model
 use lnb_param

 contains

!
!_______DERIVS_________________________________________________________________
!
 subroutine derivs(x,y,dydx)
 
 implicit none
 double complex  :: la, mu
 double precision :: x,ro,g
 double complex :: flm,rlm,rm0,rg0,frh,fgr,fn
 logical :: solid
 integer :: iy,id,n2,ny
 double complex, dimension(6)   :: y,dydx,f
 double complex, dimension(6,6) :: dy

 call earth_nlayers_viscoelas(x,ro,la,mu,g,solid)
 !print*,la,mu
!
!_______Expressions
!
 flm=(la+2.d0*mu)
 rlm=(3.d0*la+2.d0*mu)/(la+2.d0*mu)
 rm0=mu/mu0
 rg0=g/go0
 frh=ro*g*ra/mu0
 fgr=4.d0*Pi*GG*ro*ra/go0
 fn=dble(deg*(deg+1))
 
!
!_______Systeme differentiel yi dans un solide
!
 if(solid) then
 
 ny = 6

 dy(1,1)= (-2.d0*la/flm)/x
 dy(1,2)= mu0/flm
 dy(1,3)= (la*fn/flm)/x
 dy(1,4)= 0.d0
 dy(1,5)= 0.d0
 dy(1,6)= 0.d0

 dy(2,1)=  4.d0*(-frh+rm0*rlm/x)/x + ro*freq**2*ra/mu0
 dy(2,2)=(-4.d0*mu/flm)/x
 dy(2,3)= fn*(frh-2.d0*rm0*rlm/x)/x
 dy(2,4)= fn/x
 dy(2,5)= 0.d0
 dy(2,6)= -frh/rg0

 dy(3,1)= -1.d0/x
 dy(3,2)= 0.d0
 dy(3,3)= 1.d0/x
 dy(3,4)= 1.d0/rm0
 dy(3,5)= 0.d0
 dy(3,6)= 0.d0

 dy(4,1)= (frh-2.d0*rm0*rlm/x)/x
 dy(4,2)= ( -la/flm)/x
 dy(4,3)= (2.d0*rm0*(la*(2.d0*fn-1.d0)+2.d0*mu*(fn-1.d0))/flm)/(x*x) + ro*freq**2*ra/mu0
 dy(4,4)= -3.d0/x
 dy(4,5)= -(frh/rg0)/x
 dy(4,6)= 0.d0

 dy(5,1)= fgr
 dy(5,2)= 0.d0
 dy(5,3)= 0.d0
 dy(5,4)= 0.d0
 dy(5,5)= 0.d0
 dy(5,6)= 1.d0

 dy(6,1)= 0.d0
 dy(6,2)= 0.d0
 dy(6,3)= -(fgr*fn)/x
 dy(6,4)= 0.d0
 dy(6,5)= fn/(x*x)
 dy(6,6)= -2.d0/x
  
 else
 
 ny = 2

 dy(:,:)= 0.d0

 dy(1,1)= fgr/rg0
 dy(1,2)= 1.d0
 dy(2,1)= (-4.d0*(fgr/rg0)+fn/x)/x
 dy(2,2)= -2.d0/x-fgr/rg0
 
 end if
!
!_______Derivees des yi
! 
 dydx(:)=0.d0 
 do id=1,ny
  do iy=1,ny
     dydx(id)=dydx(id)+dy(id,iy)*y(iy)
  end do
 end do

    end subroutine derivs


    subroutine propagFDS(ystart,xmin,xmax,nstep)
     implicit none
     double precision :: xmin,xmax,dr,x
     integer :: nstep,i
     double complex, dimension(6)   :: ystart,dydx,y,y1

     dr = (xmax -xmin)/nstep
     y = ystart
     x = xmin

     do i = 1,nstep

  call derivs(x,y,dydx)
  y1 = y + dr*dydx
  x = x + dr
  y = y1

     end do

     ystart = y

 end subroutine propagFDS

!_____LUDCMP___________________________________________________________________
!
      SUBROUTINE ludcmp(a,n,np,indx,d)
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      INTEGER n,np,indx(n),NMAX
      REAL*8 d,a(np,np),TINY
      PARAMETER (NMAX=500,TINY=1.0D-30)
      INTEGER i,imax,j,k
      REAL*8 aamax,dum,sum,vv(NMAX)
      d=1d0
      do 12 i=1,n
        aamax=0d0
        do 11 j=1,n
          if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j))
11      continue
        if (aamax.eq.0d0) stop 'singular matrix in ludcmp'
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0d0
        do 16 i=j,n
          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*dabs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(a(j,j).eq.0d0)a(j,j)=TINY
        if(j.ne.n)then
          dum=1d0/a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      return
      END SUBROUTINE ludcmp
!
!_____LUBKSB___________________________________________________________________
!
      SUBROUTINE lubksb(a,n,np,indx,b)
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      INTEGER n,np,indx(n)
      REAL*8 a(np,np),b(np)
      INTEGER i,ii,j,ll
      REAL*8 sum
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0d0) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        do 13 j=i+1,n
          sum=sum-a(i,j)*b(j)
13      continue
        b(i)=sum/a(i,i)
14    continue
      return
      END SUBROUTINE lubksb



end module util
