module mod_potential
    use mod_parameters
    use mod_variables
    use omp_lib
    implicit none
    type :: other_moons
        real(8) :: mass, ang_vel, semimajor_a, period, rij, nij, C20, C22, S22, cospij, sinpij
        real(8) :: cos2pij, sin2pij
    end type other_moons
    type(other_moons) :: Europa, Io
contains

subroutine induced_potential
    implicit none
    integer    :: j, m, il, ir, jm
    real(8)    :: cc, jj, cj
    complex(8) :: ind_pot

    ind_pot = dzero
    do il = 1, numoflayers
        layer(il)%pot = dzero
    end do
    cc = 4.0d0 * Grav_Const * Pi

    !$omp parallel do private(j, m, ir, jm, jj, cj, ind_pot) shared(cc)
    do il = 1, numoflayers
        do j = 2, min(jfinal, 8), 2
            jj = dble(j)
            cj = 1.0d0/(2.0d0*jj+1.0d0)
            do m = 0, min(2,j), 2
                jm = jm2(j,m)
                do ir = 1, layer(il)%n
                    call integrator(layer(il)%r_i(ir), j, m, il, ir, ind_pot)
                    layer(il)%pot(ir, jm) = ind_pot * layer(il)%r_i(ir) * cc * cj
                end do
            end do
        end do
    end do
    !$omp end parallel do
end subroutine induced_potential

subroutine self_gravity
    implicit none
    integer :: j, m, il, ir, jm1, jp1, jm
    real(8) :: jj, dc1, dc2, c1, c2, hi, hup, hlo

    do il = 1, numoflayers
        do j = 2, min(jfinal, 8), 2
            jj  = dble(j)
            dc1 =  dsqrt( jj       /(2.0d0*jj+1.0d0))
            dc2 = -dsqrt((jj+1.0d0)/(2.0d0*jj+1.0d0))
            c1  =  (jj+1.0d0)*dc1
            c2  = -(jj      )*dc2
            do m = 0, min(2,j), 2
                jm1 = jml2(j,m,j-1)
                jp1 = jml2(j,m,j+1)
                jm  = jm2 (j,m)
                do ir = 1, layer(il)%n-1
                    hi  =   layer(il)%r_i(ir+1)  - layer(il)%r_i(ir)
                    hup =  (layer(il)%r_i(ir+1)  - layer(il)%r_l(ir)) / hi
                    hlo =  (layer(il)%r_l(ir  )  - layer(il)%r_i(ir)) / hi
                    layer(il)%grad_pot(ir, jm1) = dc1 * (layer(il)%pot(ir+1, jm)     - layer(il)%pot(ir, jm)    )/hi + &
                                                   c1 * (layer(il)%pot(ir+1, jm)*hlo + layer(il)%pot(ir, jm)*hup)/(layer(il)%r_l(ir+1))
                    layer(il)%grad_pot(ir, jp1) = dc2 * (layer(il)%pot(ir+1, jm)     - layer(il)%pot(ir, jm)    )/hi + &
                                                   c2 * (layer(il)%pot(ir+1, jm)*hlo + layer(il)%pot(ir, jm)*hup)/(layer(il)%r_l(ir+1))
                end do
            end do
        end do
    end do
end subroutine self_gravity


subroutine integrator(r, j, m, layer_number, radial_number, ind_pot)
    implicit none
    integer    :: layer_number, radial_number, il, ir, j, m, jm, nr, nn
    real(8)    :: cint0, cint1, cext0, cext1, r, drho, dr
    complex(8) :: ind_pot


    !! INTEGRATION OF THE PARTS THAT ARE BELOW THE GIVEN POINT
    jm = jm2(j, m)
    ind_pot = 0.0d0
    do il = 1, layer_number
        if(layer_number > 1) then
            if(il < layer_number) then
                nr = layer(il)%n
            else 
                nr = radial_number
            end if
        else
            nr = radial_number
        end if
        do ir = 1, nr-1
            dr = layer(il)%r_i(ir+1) - layer(il)%r_i(ir)
            cint0 = (layer(il)%r_i(ir)  /r)  ** (j+2)
            cint1 = (layer(il)%r_i(ir+1)/r)  ** (j+2)
            ind_pot = ind_pot + 0.5d0 * dr*  (cint0 * layer(il)%rhoe(ir, jm) + cint1 * layer(il)%rhoe(ir+1, jm))
        end do
    end do
    if(layer_number > 1) then
        do il = 1, layer_number-1
            nn      = layer(il)%n
            drho    = layer(il)%rho_i(nn) - layer(il+1)%rho_i(1)
            ind_pot = ind_pot + drho * (layer(il)%ur_i(nn, jm)) * (layer(il)%r_i(nn)  /r)  ** (j+2)
        end do
    end if

    !! INTEGRATION OF THE PARTS THAT ARE ABOVE THE GIVEN POINT
    do il = layer_number, numoflayers
        if(layer_number < numoflayers) then
            if(il > layer_number) then
                nr = 1
            else
                nr = radial_number
            end if
        else
            nr = radial_number
        end if
        do ir = nr, layer(il)%n-1
            dr = layer(il)%r_i(ir+1) - layer(il)%r_i(ir)
            cext0 = (r/layer(il)%r_i(ir)  )** (j-1)
            cext1 = (r/layer(il)%r_i(ir+1))** (j-1)
            ind_pot = ind_pot + 0.5d0 * dr * (cext0 * layer(il)%rhoe(ir, jm) + cext1 * layer(il)%rhoe(ir+1, jm))
        end do
    end do
    
    do il = layer_number, numoflayers
        nn = layer(il)%n
        if((il<numoflayers)) then
            drho    = layer(il)%rho_i(nn) - layer(il+1)%rho_i(1)
        elseif(il==numoflayers)then
            drho    = layer(il)%rho_i(nn)
        end if
        ind_pot = ind_pot + drho * (layer(il)%ur_i(nn, jm)) * (r/layer(il)%r_i(nn)  )** (j-1)
    end do
end subroutine integrator


subroutine love_numbers
    implicit none
    integer    :: nn, j
    complex(8) :: surface_pot

    nn = layer(numoflayers)%n

    surface_pot = layer(numoflayers)%pot(nn, jm2(2,0))
    k20 = real(surface_pot)

    surface_pot = layer(numoflayers)%pot(nn, jm2(2,2))
    real_k22    = real(surface_pot)
    imag_k22    = imag(surface_pot)

    h20         = real(layer(numoflayers)%ur_i(nn, jm2(2,0)))!/real(forcing_pot/layer(numoflayers)%g_i(nn))
    real_h22    = real(layer(numoflayers)%ur_i(nn, jm2(2,2)))
    imag_h22    = imag(layer(numoflayers)%ur_i(nn, jm2(2,2)))

    if(SWITCH_velocity==1)then
    l20         = real(layer(numoflayers)%dis(nn, jml2(2,0,1))*dsqrt(3.0d0/5.0d0) + layer(numoflayers)%dis(nn, jml2(2,0,3))*dsqrt(2.0d0/5.0d0))
    real_l22    = real(layer(numoflayers)%dis(nn, jml2(2,2,1))*dsqrt(3.0d0/5.0d0) + layer(numoflayers)%dis(nn, jml2(2,2,3))*dsqrt(2.0d0/5.0d0))
    imag_l22    = imag(layer(numoflayers)%dis(nn, jml2(2,2,1))*dsqrt(3.0d0/5.0d0) + layer(numoflayers)%dis(nn, jml2(2,2,3))*dsqrt(2.0d0/5.0d0))
    else
    l20         = real(layer(numoflayers)%vel(nn, jml2(2,0,1))*dsqrt(3.0d0/5.0d0) + layer(numoflayers)%vel(nn, jml2(2,0,3))*dsqrt(2.0d0/5.0d0))
    real_l22    = real(layer(numoflayers)%vel(nn, jml2(2,2,1))*dsqrt(3.0d0/5.0d0) + layer(numoflayers)%vel(nn, jml2(2,2,3))*dsqrt(2.0d0/5.0d0))
    imag_l22    = imag(layer(numoflayers)%vel(nn, jml2(2,2,1))*dsqrt(3.0d0/5.0d0) + layer(numoflayers)%vel(nn, jml2(2,2,3))*dsqrt(2.0d0/5.0d0))
    end if


    do j = 1, 3
        kh(3*(j-1)+1) = real(layer(numoflayers)%pot(nn, jm2(2*(j+1), 0))) !!4,0 
        kh(3*(j-1)+2) = real(layer(numoflayers)%pot(nn, jm2(2*(j+1), 2))) !!4,2, real 
        kh(3*(j-1)+3) = imag(layer(numoflayers)%pot(nn, jm2(2*(j+1), 2))) !!4,2, imag 
    end do
end subroutine love_numbers


complex(8) function tides(den, rr, m, tt) result (tid)
    implicit none
    integer    :: m
    real(8)    :: c20, rec22, imc22, den, rr, tt

    c20    =  ang_vel * ang_vel * ecc * dsqrt(18.0d0 * Pi) * cos(ang_vel * tt)
    rec22  =  ang_vel * ang_vel * ecc * dsqrt(27.0d0 * Pi) * cos(ang_vel * tt)
    imc22  =  ang_vel * ang_vel * ecc * dsqrt(48.0d0 * Pi) * sin(ang_vel * tt)
    
    tid = dzero
    if    (m.eq.0)then
        tid =  den*rr*dcmplx(c20,0.0d0)
    elseif(m.eq.1)then
    else
        tid = -den*rr*dcmplx(rec22,-imc22)
    end if
end function tides

complex(8) function mm_tides(den, rr, m, tt) result (tid)
    implicit none
    integer    :: m
    real(8)    :: den, rr, tt
    
    Europa%mass    = 4.79984d22
    Europa%period  = rotation_period*0.5d0
    Europa%ang_vel = 2.0d0*pi/Europa%period
    Europa%semimajor_a = 671100.d3
    Io%mass        = 8.93d22
    Io%period      = rotation_period*0.25d0
    Io%ang_vel     = 2.0d0*pi/Io%period
    Io%semimajor_a = 421700.d3
    Europa%nij = -(Europa%ang_vel-ang_vel)
    Io%nij     = -(Io%ang_vel-ang_vel)
    Europa%rij = dsqrt(Europa%semimajor_a**2 + semimajor_a**2 - 2.0d0*Europa%semimajor_a*semimajor_a*cos(Europa%nij*tt))
    Io%rij     = dsqrt(Io%semimajor_a**2     + semimajor_a**2 - 2.0d0*Io%semimajor_a    *semimajor_a*cos(Io%nij*tt - Pi))!!with initial delay of Pi

    Europa%cospij     = (Europa%semimajor_a - semimajor_a * cos(Europa%nij * tt)) / Europa%rij
    Europa%sinpij     =                       semimajor_a * sin(Europa%nij * tt)  / Europa%rij

    Io%cospij         = (Io%semimajor_a - semimajor_a * cos(Io%nij * tt)) / Io%rij
    Io%sinpij     =                       semimajor_a * sin(Io%nij * tt - Pi)  / Io%rij !!with initial delay of Pi

    Europa%cos2pij = Europa%cospij**2 - Europa%sinpij**2
    Europa%sin2pij = Europa%sinpij*Europa%cospij*2.0d0

    Io%cos2pij = Io%cospij**2 - Io%sinpij**2
    Io%sin2pij = Io%sinpij*Io%cospij*2.0d0

    !!Compute the gradient...
    Europa%C20 = dsqrt(2.0d0 * Pi) * Grav_Const * Europa%mass / (Europa%rij**3) !dsqrt(pi/5.0d0)*Grav_Const*Europa%mass/(Europa%rij**3)
    Europa%C22 = dsqrt(3.0d0 * Pi) * Grav_Const * Europa%mass * Europa%cos2pij / (Europa%rij**3)
    Europa%S22 = dsqrt(3.0d0 * Pi) * Grav_Const * Europa%mass * Europa%sin2pij / (Europa%rij**3)

    Io%C20     = dsqrt(2.0d0 * Pi) * Grav_Const * Io%mass     / (Io%rij**3) 
    Io%C22     = dsqrt(3.0d0 * Pi) * Grav_Const * Io%mass * Io%cos2pij / (Io%rij**3)
    Io%S22     = dsqrt(3.0d0 * Pi) * Grav_Const * Io%mass * Io%sin2pij / (Io%rij**3)

    tid = dzero
    if    (m.eq.0)then
        tid =  ( den*rr*dcmplx(Europa%C20+Io%C20, 0.0d0))
    elseif(m.eq.1)then
    else
        tid =  (-den*rr*dcmplx(Europa%C22+Io%C22, -Europa%S22-Io%S22))
    end if
end function mm_tides

end module mod_potential