import numpy as np
from constants import *
from get_matrixes import *

def calc_kn(n,r,rho):

   grav_acc = __get_g__(r,rho)
   lmax = len(r)
   
   # Create full propagation matrix (Full_Mat)
   # Layer (L) matrix and inverse
   L,Linv = fundamental_matrix(n)
   # Core (Lcore) matrix (inverse not needed)
   Lcore = core_propagator_matrix(n,r[0]/r[-1])
   Full_Mat = Lcore

   # Cycle over layers above the core
   for j in range(1,lmax):
      # Boundary matrix
      delrho   = rho[j-1] - rho[j]
      Blayer   = boundary_matrix(grav_acc[j-1],delrho,r[j-1])
      Full_Mat = np.matmul(Blayer,Full_Mat)
      # Layer matrix
      X        = radial_matrix(n, r[j]/r[j-1])
      Player   = np.matmul(L,np.matmul(X,Linv))
      Full_Mat = np.matmul(Player,Full_Mat)

   # Build final matrix equality
   lhs = np.zeros((2,2))
   lhs[0,0] = Full_Mat[0,0]
   lhs[1,0] = Full_Mat[1,0]
   lhs[0,1] = -1.
   lhs[1,1] = n + 1. - 4*np.pi*grav_const*rho[-1]*r[-1]/grav_acc[-1]
   rhs      = np.zeros((2,1))
   rhs[1,0] = (2*n + 1) 
   # Invert to find solution
   Solution = np.matmul(np.linalg.inv(lhs),rhs)
   kn = Solution[1] - 1.
   return kn


def __get_g__(r,rho):
   # Calculate gravity vector
   pref     = (4*np.pi/3)*grav_const 
   diff3    = np.diff(np.append(0,r)**3.)
   g = pref*np.divide(np.cumsum(np.multiply(rho,diff3)),r**2.)
   return g