"""
The potential fields of a triaxial ellipsoid.
"""

from __future__ import division
import numpy as np
from scipy import linalg

from ..constants import SI2MGAL, G, CM, T2NT, SI2EOTVOS, PERM_FREE_SPACE
from .. import utils

import scipy.special

def bx_c(xp,yp,zp,inten,inc,dec,ellipsoids):
    '''
    Calculates the X component of the magnetic field generated by n-ellipsoid.
    
    Parameters:
    
    * xp,yp,zp: arrays
        Grid of observation points x, y, and z.
    * inten,inc,dec: floats
        Intensity, inclination and declination of the Earth's magnetic field.
    * ellipsoids: list of :class:`fatiando.mesher.Sphere` 
        Ellipsoid model.
    
    Returns:
    
    * bx: array
        The X component of the magnetic field generated by n-ellipsoid.
    '''
    
    if xp.shape != yp.shape != zp.shape:
        raise ValueError("Input arrays xp, yp, and zp must have same shape!")
    size = len(xp)
    res = np.zeros(size, dtype=np.float)
    ctemag = 100.
    
    for ellipsoid in ellipsoids:
        b1,b2,b3,V,N1,N2,N3,JRD_ang = ellipsoid_def (xp,yp,zp,inten,inc,dec,ellipsoid)
        bx = b1*V[0,0]+b2*V[0,1]+b3*V[0,2]
        res += bx 
    res = res*ctemag
    return res
    
def by_c(xp,yp,zp,inten,inc,dec,ellipsoids):
    '''
    Calculates the Y component of the magnetic field generated by n-ellipsoid.
    
    Parameters:
    
    * xp,yp,zp: arrays
        Grid of observation points x, y, and z.
    * inten,inc,dec: floats
        Intensity, inclination and declination of the Earth's magnetic field.
    * ellipsoids: list of :class:`fatiando.mesher.Sphere` 
        Ellipsoid model.
    
    Returns:
    
    * by: array
        The Y component of the magnetic field generated by n-ellipsoid.
    '''
    
    if xp.shape != yp.shape != zp.shape:
        raise ValueError("Input arrays xp, yp, and zp must have same shape!")
    size = len(xp)
    res = np.zeros(size, dtype=np.float)
    ctemag = 100.
    
    for ellipsoid in ellipsoids:
        b1,b2,b3,V,N1,N2,N3,JRD_ang = ellipsoid_def (xp,yp,zp,inten,inc,dec,ellipsoid)
        by = b1*V[1,0]+b2*V[1,1]+b3*V[1,2]
        res += by
    res = res*ctemag
    return res

def bz_c(xp,yp,zp,inten,inc,dec,ellipsoids):
    '''
    Calculates the Z component of the magnetic field generated by n-ellipsoid.
    
    Parameters:
    
    * xp,yp,zp: arrays
        Grid of observation points x, y, and z.
    * inten,inc,dec: floats
        Intensity, inclination and declination of the Earth's magnetic field.
    * ellipsoids: list of :class:`fatiando.mesher.Sphere` 
        Ellipsoid model.
    
    Returns:
    
    * bz: array
        The Z component of the magnetic field generated by n-ellipsoid.
    '''
    
    if xp.shape != yp.shape != zp.shape:
        raise ValueError("Input arrays xp, yp, and zp must have same shape!")
    size = len(xp)
    res = np.zeros(size, dtype=np.float)
    ctemag = 100.
    
    for ellipsoid in ellipsoids:
        b1,b2,b3,V,N1,N2,N3,JRD_ang = ellipsoid_def (xp,yp,zp,inten,inc,dec,ellipsoid)
        bz = b1*V[2,0]+b2*V[2,1]+b3*V[2,2]
        res += bz
    res = res*ctemag
    return res
    
def tf_c(xp,yp,zp,inten,inc,dec,ellipsoids):
    '''
    Calculates the approximated total-field anomaly generated by n-ellipsoid.
    
    Parameters:
    
    * xp,yp,zp: arrays
        Grid of observation points x, y, and z.
    * inten,inc,dec: floats
        Intensity, inclination and declination of the Earth's magnetic field.
    * ellipsoids: list of :class:`fatiando.mesher.Sphere` 
        Ellipsoid model.
    
    Returnss:

    * tf : array
        The total-field anomaly
    '''
    
    if xp.shape != yp.shape != zp.shape:
        raise ValueError("Input arrays xp, yp, and zp must have same shape!")
    size = len(xp)
    res = np.zeros(size, dtype=np.float)
    ctemag = 100.
    
    for ellipsoid in ellipsoids:
        b1,b2,b3,V,N1,N2,N3,JRD_ang = ellipsoid_def (xp,yp,zp,inten,inc,dec,ellipsoid)
        bx = (b1*V[0,0]+b2*V[0,1]+b3*V[0,2])*ctemag
        by = (b1*V[1,0]+b2*V[1,1]+b3*V[1,2])*ctemag
        bz = (b1*V[2,0]+b2*V[2,1]+b3*V[2,2])*ctemag
        tf = bx*np.cos(np.deg2rad(inc))*np.cos(np.deg2rad(dec)) + by*np.cos\
        (np.deg2rad(inc))*np.sin(np.deg2rad(dec)) + bz*np.sin(np.deg2rad(inc))
        res += tf
    res = res
    return res,N1,N2,N3,JRD_ang
    
def ellipsoid_def (xp,yp,zp,inten,inc,dec,ellipsoid):
    '''
    Calculate the potential fields of a homogeneous ellipsoid.

    **Magnetic**

    Calculates the magnetic effect produced by a triaxial, a prolate or/and 
    an oblate ellipsoid. The functions are
    based on Clark et al. (1986).
    '''

    axis = ellipsoid.axis
    center = ellipsoid.center
    alpha, delta, gamma = structural_angles(ellipsoid.strike, ellipsoid.dip, ellipsoid.rake)
    V = ellipsoid.V(angles = [alpha, delta, gamma])
    
    # Remanence values
    intensity_rem = ellipsoid.props['remanence'][0]
    incli_rem = ellipsoid.props['remanence'][1]
    decli_rem = ellipsoid.props['remanence'][2]
    ln, mn, nn = utils.dircos(incli_rem, decli_rem)
    
    k_int = np.array([ellipsoid.props['k'][0],ellipsoid.props['k'][1],\
    ellipsoid.props['k'][2]])
    #isotropic case
    if k_int[0] == (k_int[1] and k_int[2]):
        km = np.diag(k_int)
    #anisotropic case
    else:
        U = ellipsoid.V(angles = [ellipsoid.props['k'][3],\
        ellipsoid.props['k'][4],ellipsoid.props['k'][5]])
        km = k_matrix(U,V,np.diag(k_int))
    print km    
    # Ellipsoid cartesian body coordinates
    x1,x2,x3 = x_e(xp,yp,zp,center,V)
    
    # Largest real root of the cubic equation (Lambda)
    lamb,teta,q,p,p2,p1,p0 = lamb_T(axis,x1,x2,x3)
    
    # Derivatves of lambda
    dlambx1,dlambx2,dlambx3 = dlambx_T(axis,x1,x2,x3,lamb)
    
    # Calculate the eliptical integral parameters
    F,E,F2,E2,k,theta_l = legendre_integrals(axis,lamb)
    N1,N2,N3 = N_desmagT(axis,F2,E2)
    
    # Integrals calculations
    A, B, C = potential_integrals(axis,k,theta_l,F,E)
    
    # Geometry for the magnetic field
    m11,m12,m13,m21,m22,m23,m31,m32,m33, cte, V1, V2, V3 = mx(axis,x1,x2,x3,\
    dlambx1,dlambx2,dlambx3,A,B,C,lamb)
    
    # Earth's field and total body magnetization (including demagnetization) 
    #in the body's coordinate
    JN = JN_e (intensity_rem,ln,mn,nn,V)
    lt,mt,nt = utils.dircos(inc, dec)
    Ft = F_e (inten,lt,mt,nt,V)
    JR = JR_e (km,JN,Ft)
    JRD = JRD_e (km,N1,N2,N3,JR)
    JRD_carte = (V).dot(JRD)
    JRD_ang = utils.vec2ang(JRD_carte)
    
    # Components of the magnetic field in the body coordinates
    B1 = B1_e (m11,m12,m13,JRD,axis[0],axis[1],axis[2])
    B2 = B2_e (m21,m22,m23,JRD,axis[0],axis[1],axis[2])
    B3 = B3_e (m31,m32,m33,JRD,axis[0],axis[1],axis[2])
    return B1,B2,B3,V,N1,N2,N3,JRD_ang

def structural_angles(strike, dip, rake):
    '''
    Calculates the orientation angles alpha, gamma
    and delta (Clark et al., 1986)
    as functions of the geological angles strike, dip and
    rake (Clark et al., 1986; Allmendinger et al., 2012).
    The function implements the formulas presented by
    Clark et al. (1986).

    Parameters:

    *strike: float
             strike direction (in degrees).
    *dip: float
          true dip (in degrees).
    *rake: float
           angle between the strike and the semi-axis a
           of the body (in degrees).

    Returns:

    *alpha, gamma, delta: float, float, float
            orientation angles (in radians) defined according
            to Clark et al. (1986).

    References:

    Clark, D., Saul, S., and Emerson, D.: Magnetic and gravity anomalies
    of a triaxial ellipsoid, Exploration Geophysics, 17, 189-200, 1986.

    Allmendinger, R., Cardozo, N., and Fisher, D. M.:
    Structural geology algorithms : vectors and tensors,
    Cambridge University Press, 2012.
    '''

    strike_r = np.deg2rad(strike)
    cos_dip = np.cos(np.deg2rad(dip))
    sin_dip = np.sin(np.deg2rad(dip))
    cos_rake = np.cos(np.deg2rad(rake))
    sin_rake = np.sin(np.deg2rad(rake))

    aux = sin_dip*sin_rake
    aux1 = cos_rake/np.sqrt(1 - aux*aux)
    aux2 = sin_dip*cos_rake

    if aux1 > 1.:
        aux1 = 1.
    if aux1 < -1.:
        aux1 = -1.

    alpha = strike_r - np.arccos(aux1)
    if aux2 != 0:
        gamma = -np.arctan(cos_dip/aux2)
    else:
        if cos_dip >= 0:
            gamma = np.pi/2
        if cos_dip <= 0:
            gamma = -np.pi/2
    delta = np.arcsin(aux)

    assert delta <= np.pi/2, 'delta must be lower than \
or equalt to 90 degrees'

    assert (gamma >= -np.pi/2) and (gamma <= np.pi/2), 'gamma must lie between \
-90 and 90 degrees.'

    return alpha, gamma, delta
	
def x_e (xp,yp,zp,center,V):
        '''
        Calculates the new coordinates with origin at the center 
        of the ellipsoid.

        Parameters:
        
        * xp,yp,zp: arrays
            Grid of observation points x, y, and z.
        * center: float
            Origin of the center of the ellipsoid.
        * V: array
            Matrix of conversion.
        
        Returns:
        
        * x1, x2, x3: arrays
            The three grid points of the body's coordinates.
        '''
        x1 = (xp-center[0])*V[0,0]+(yp-center[1])*V[1,0]-(zp+center[2])*V[2,0]
        x2 = (xp-center[0])*V[0,1]+(yp-center[1])*V[1,1]-(zp+center[2])*V[2,1]
        x3 = (xp-center[0])*V[0,2]+(yp-center[1])*V[1,2]-(zp+center[2])*V[2,2]
        return x1, x2, x3
        
def JN_e (intensity_rem,ln,mn,nn,V):
        '''
        Changes the remanent magnetization vector to the body coordinate.
        
        Parameters:
        
        * ln,nn,mn:
            Direction cosines of the remanent magnetization vector.
        * V:
            Matrix of conversion.
        
        Returns:
        
        * JN:
            Remanent magnetization vector in the body coordinate.         
        '''
        
        JN = intensity_rem*np.array([[(ln*V[0,0]+mn*V[1,0]+nn*V[2,0])], \
        [(ln*V[0,1]+mn*V[1,1]+nn*V[2,1])], [(ln*V[0,2]+mn*V[1,2]+nn*V[2,2])]])
        return JN

def N_desmagT (axis,F2,E2):
        '''
        Calculates the three demagnetization factor along major, 
        intermediate and minor axis.
        
        Parameters:
        
        * a,b,c: float
            Major, intermediate and minor axis, respectively.
        * F2, E2: float
            Lagrange's normal eliptic integrals of first and second order.
        
        Returns:
        
        * N1, N2, N3: floats
            Major, intermediate and minor demagnetization factors, respectively.        
        '''
        
        N1 = ((axis[0]*axis[1]*axis[2])/((axis[0]**2-axis[1]**2)*\
        (axis[0]**2-axis[2]**2)**0.5)) * (F2-E2)
        
        N2 = (((axis[0]*axis[1]*axis[2])*\
        (axis[0]**2-axis[2]**2)**0.5)/((axis[0]**2-axis[1]**2)*\
        (axis[1]**2-axis[2]**2))) * (E2-((axis[1]**2-axis[2]**2)\
        /(axis[0]**2-axis[2]**2)) * F2-((axis[2]*(axis[0]**2-axis[1]**2))\
        /(axis[0]*axis[1]*(axis[0]**2-axis[2]**2)**0.5)))
        
        N3 = ((axis[0]*axis[1]*axis[2])/((axis[1]**2-axis[2]**2)*\
        (axis[0]**2-axis[2]**2)**0.5)) * (((axis[1]*(axis[0]**2-axis[2]**2)**0.5)\
        /(axis[0]*axis[2]))-E2)

        return N1, N2, N3

def k_matrix (U,V,K):
        '''
        Build the susceptibility tensor for the anisotropic case.
        
        Parameters:
        
        * U: array
            Direction cosines of the susceptibilities.
        * V: array
            Matrix of coordinates conversion.
        * K: array
            Diagonal matrix with k1,k2,k3 (intensity of the susceptibilities).
        
        Returns:
        
        * km: array
            Susceptibility tensors matrix.
        '''
        
        km = np.dot(np.dot(np.dot(V.T,U), K), np.dot(U.T,V))
        return km

def lamb_T (axis,x1,x2,x3):
        '''
        Calculates the larger root of the cubic equation: 
        s^3 + p2*s^2 + p1*s + p0 = 0.
        
        Parameters:
        
        * a, b, c: floats
            Major, intermediate and minor axis, respectively.
        * x1, x2, x3: arrays
            Axis of the body coordinate system.
        
        Returns:
        
        * lamb: array
            Larger root of the cubic equation: s^3 + p2*s^2 + p1*s + p0 = 0.
        * teta, q, p, p2, p1, p0: arrays
            Constants of the cubic equation.        
        '''
        
        p0 = (axis[0]*axis[1]*axis[2])**2-(axis[1]*axis[2]*x1)**2-\
        (axis[2]*axis[0]*x2)**2-(axis[0]*axis[1]*x3)**2
        
        p1 = (axis[0]*axis[1])**2+(axis[1]*axis[2])**2+(axis[2]*axis[0])**2-\
        (axis[1]**2+axis[2]**2)*x1**2-(axis[2]**2+axis[0]**2)*x2**2-\
        (axis[0]**2+axis[1]**2)*x3**2
        
        p2 = axis[0]**2+axis[1]**2+axis[2]**2-x1**2-x2**2-x3**2
        p = p1-(p2**2)/3.
        q = p0-((p1*p2)/3.)+2*(p2/3.)**3
        p3 = (-q/(2*np.sqrt((-p/3.)**3)))
        for i in range (len(p3)):
            if p3[i] > 1.:
                p3[i] = 1.
        teta = np.arccos(p3)
        lamb = 2.*((-p/3.)**0.5)*np.cos(teta/3.)-(p2/3.)
        return lamb, teta, q, p, p2, p1, p0

def legendre_integrals(axis,lamb):
        '''
        Calculates parameters and the Legendre's normal elliptic integrals 
        of first and second order.
        
        Parameters:
        
        * a, b, c: floats
            Major, intermediate and minor axis, respectively.
        * lamb: array
            Larger root of the cubic equation: s^3 + p2*s^2 + p1*s + p0 = 0.
        
        Returns:
        
        * F: array
            Legendre's normal elliptic integrals of first order.
        * E: array
            Legendre's normal elliptic integrals of second order.
        * F2: array
            Legendre's normal elliptic integrals of first order.
        * E2: array
            Legendre's normal elliptic integrals of second order.
        * k: array
            Legendre's normal elliptic integrals parameter.
        * theta_l: array
            Legendre's normal elliptic integrals parameter.        
        '''
        
        k = np.zeros_like(lamb)
        k1 = ((axis[0]**2-axis[1]**2)/(axis[0]**2-axis[2]**2))**0.5
        k.fill(k1)
        theta_l = np.arcsin(((axis[0]**2-axis[2]**2)/(axis[0]**2+lamb))**0.5)
        theta_l2 = np.arccos(axis[2]/axis[0])
        F = scipy.special.ellipkinc(theta_l, k**2)
        E = scipy.special.ellipeinc(theta_l, k**2)
        F2 = scipy.special.ellipkinc(theta_l2, k1**2)
        E2 = scipy.special.ellipeinc(theta_l2, k1**2)
        return F,E,F2,E2,k,theta_l

def dlambx_T (axis,x1,x2,x3,lamb):
        '''
        Calculates the derivatives of the ellipsoid equation for 
        each body coordinates in realation to lambda.
        
        Parameters:
        
        * a, b, c: floats
            Major, intermediate and minor axis, respectively.
        * x1, x2, x3: array
            Axis of the body coordinate system.
        * lamb: array
            Larger root of the cubic equation: s^3 + p2*s^2 + p1*s + p0 = 0.
        
        Returns:
        
        * dlambx1,dlambx2,dlambx3: array
            Derivatives of the ellipsoid equation for each body coordinates 
            in realation to x1,x2 and x3.        
        '''
        
        dlambx1 = (2*x1/(axis[0]**2+lamb))/((x1/(axis[0]**2+lamb))**2+\
        (x2/(axis[1]**2+lamb))**2+((x3/(axis[2]**2+lamb))**2))        
        
        dlambx2 = (2*x2/(axis[1]**2+lamb))/((x1/(axis[0]**2+lamb))**2+\
        (x2/(axis[1]**2+lamb))**2+((x3/(axis[2]**2+lamb))**2))
        
        dlambx3 = (2*x3/(axis[2]**2+lamb))/((x1/(axis[0]**2+lamb))**2+\
        (x2/(axis[1]**2+lamb))**2+((x3/(axis[2]**2+lamb))**2))
        return dlambx1, dlambx2, dlambx3

def potential_integrals(axis,k,theta_l,F,E):
        '''
        Calculates the integrals which is part of the solution of 
        the potential field of an homogeneous ellipsoid (Dirichlet,1839).
        
        Parameters:
        
        * a, b, c: floats
            Major, intermediate and minor axis, respectively.
        * k: array
            Legendre's normal elliptic integrals parameter.
        * theta_l: array
            Legendre's normal elliptic integrals parameter.
        * F: array
            Legendre's normal elliptic integrals of first order.
        * E: array
            Legendre's normal elliptic integrals of second order.
        
        Returns:
        
        * A2,B2,C2: array
            Integrals of the potential field of an homogeneous ellipsoid.        
        '''
        
        A2 = (2/((axis[0]**2-axis[1]**2)*(axis[0]**2-axis[2]**2)**0.5))*(F-E)
        
        B2 = ((2*(axis[0]**2-axis[2]**2)**0.5)/((axis[0]**2-axis[1]**2)*\
        (axis[1]**2-axis[2]**2)))*(E-((axis[1]**2-axis[2]**2)\
        /(axis[0]**2-axis[2]**2))*F-((k**2*np.sin(theta_l)*np.cos(theta_l))\
        /(1-k**2*np.sin(theta_l)*np.sin(theta_l))**0.5))
        
        C2 = (2/((axis[1]**2-axis[2]**2)*(axis[0]**2-axis[2]**2)**0.5))\
        *(((np.sin(theta_l)*((1-k**2*np.sin(theta_l)*np.sin(theta_l))**0.5))\
        /np.cos(theta_l))-E)
        return A2,B2,C2

def mx(axis,x1,x2,x3,dlambx1,dlambx2,dlambx3,A,B,C,lamb):
        '''
        Additional calculations for the ellipsoid magnetic field.
        
        Parameters:
        
        * a, b, c: floats
            Major, intermediate and minor axis, respectively.
        * x1, x2, x3: array
            Axis of the body coordinate system.
        * A,B,C: array
            Integrals of the potential field of an homogeneous ellipsoid
        * lamb: array
            Larger root of the cubic equation: s^3 + p2*s^2 + p1*s + p0 = 0.
        
        Returns: 
        
        * m11, m12, m13, m21, m22, m23, m31, m32, m33, cte, V1, V2, V3: array
            Calculus for the ellipsoid magnetic field.        
        '''
        
        cte = 1/np.sqrt((axis[0]**2+lamb)*(axis[1]**2+lamb)*(axis[2]**2+lamb))
        V1 = x1/(axis[0]**2+lamb)
        V2 = x2/(axis[1]**2+lamb)
        V3 = x3/(axis[2]**2+lamb)
        m11 = (cte*dlambx1*V1)-A
        m12 = cte*dlambx1*V2
        m13 = cte*dlambx1*V3
        m21 = cte*dlambx2*V1
        m22 = (cte*dlambx2*V2)-B
        m23 = cte*dlambx2*V3
        m31 = cte*dlambx3*V1
        m32 = cte*dlambx3*V2
        m33 = (cte*dlambx3*V3)-C
        return m11, m12, m13, m21, m22, m23, m31, m32, m33, cte, V1, V2, V3
        
def jrd_cartesiano (inten,inc,dec,ellipsoids):
    '''
    Calculates the intensity and direction of the resultant 
    vector of magnetization.
    
    Parameters:
    
    * inten:
        Intensity of the Earth's magnetic field.
    * inc:
        Inclination of the Earth's magnetic field.
    * dec:
        Declination of the Earth's magnetic field.
    * ellipsoid:
        magnetic ellipsoid model.
    
    Returns:
    
    * JRD_ang:
        Vector with intensity and direction of the resultant 
        vector of magnetization 
        in the cartesian coordinates(degrees).    
    '''
    
    inc = np.deg2rad(inc)
    dec = np.deg2rad(dec)
    lt,mt,nt = utils.dircos (dec, inc)
    Ft = []
    JR = []
    JRD = []
    JRD_carte = []
    JRD_ang = []
    for ellipsoid in ellipsoids:
    
        Ft.append(F_e (inten,lt,mt,nt,ellipsoids[i].mcon[0,0],\
        ellipsoids[i].mcon[1,0],ellipsoids[i].mcon[2,0],\
        ellipsoids[i].mcon[0,1],ellipsoids[i].mcon[1,1],\
        ellipsoids[i].mcon[2,1],ellipsoids[i].mcon[0,2],\
        ellipsoids[i].mcon[1,2],ellipsoids[i].mcon[2,2]))
        
        JR.append(JR_e (ellipsoids[i].km,ellipsoids[i].JN,Ft[i]))
        
        JRD.append(JRD_e (ellipsoids[i].km,ellipsoids[i].N1,\
        ellipsoids[i].N2,ellipsoids[i].N3,JR[i]))
        
        JRD_carte.append((ellipsoids[i].mconT).dot(JRD[i]))
        
        JRD_ang.append(utils.vec2ang(JRD_carte[i]))
    return JRD_ang
    
def F_e (inten,lt,mt,nt,V):
    '''
    Change the magnetization vetor of the Earth's field 
    to the body coordinates.
    
    Parameters:
    
    * inten: float
        Intensity of the Earth's magnetic field.
    * lt,mt,nt: floats
        Direction cosines of the Earth's magnetic field.
    * V: array
        Matrix of body coordinates change.
    
    Returns:
    
    * Ft: array
    The magnetization vetor of the Earth's field to the body coordinates.    
    '''
	
    intT= inten/(4*np.pi*100)
    Ft = intT*np.array([[(lt*V[0,0]+mt*V[1,0]+nt*V[2,0])], \
    [(lt*V[0,1]+mt*V[1,1]+nt*V[2,1])], [(lt*V[0,2]+mt*V[1,2]+nt*V[2,2])]])
    return Ft
    
def JR_e (km,JN,Ft):
    '''
    Calculates the resultant magnetization vector without 
    self-demagnetization correction.
    
    Parameters:
    
    * km: array
        matrix of susceptibilities tensor.
    * JN: array
        Remanent magnetization
    * Ft: array
        Magnetization vetor of the Earth's field in the body coordinates.
    
    Returns:
    
    * JR: array
        Resultant magnetization vector without 
        self-demagnetization correction.  
    '''
    
    JR = km.dot(Ft) + JN
    return JR
    
def JRD_e (km,N1,N2,N3,JR):
    '''
    Calculates resultant magnetization vector with 
    self-demagnetization correction.
    
    Parameters:
    
    * km: array
        matrix of susceptibilities tensor.
    * N1,N2,N3: floats
        Demagnetization factors in relation to a, b and c, respectively.
    * JR: array
        resultant magnetization vector without self-demagnetization correction.
    
    Returns:
    
    * JRD: array
        Resultant magnetization vector without self-demagnetization correction.   
    '''
    
    I = np.identity(3)
    kn0 = km[:,0]*N1
    kn1 = km[:,1]*N2
    kn2 = km[:,2]*N3
    kn = (np.vstack((kn0,kn1,kn2))).T
    A = I + kn
    JRD = (linalg.inv(A)).dot(JR)
    return JRD

def B1_e (m11,m12,m13,J,a,b,c):
    '''
    Calculates the B1 component of the magnetic field generated 
    by n-ellipsoids in the body coordinates.
    
    Parameters:
    
    * m21,m22,m23: array
        Calculus for the ellipsoid magnetic field.
    * J: array
        Resultant magnetization vector without self-demagnetization correction.
    * a,b,c: floats
        Major, intermediate and minor axis, respectively.
    
    Returns:
    
    * B1: array
        The B1 component of the magnetic field generated 
        by n-ellipsoids in the body coordinates.    
    '''
    
    B1 = 2*np.pi*a*b*c*(m11*J[0]+m12*J[1]+m13*J[2])
    return B1

def B2_e (m21,m22,m23,J,a,b,c):
    '''
    Calculates the B2 component of the magnetic field generated by 
    n-ellipsoids in the body coordinates.
    
    Parameters:
    
    * m21,m22,m23: array
        Calculus for the ellipsoid magnetic field.
    * J: array
        Resultant magnetization vector without self-demagnetization correction.
    * a,b,c: floats
        Major, intermediate and minor axis, respectively.
    
    Returns:
    
    * B2: array
        The B2 component of the magnetic field generated 
        by n-ellipsoids in the body coordinates.    
    '''
    
    B2 = 2*np.pi*a*b*c*(m21*J[0]+m22*J[1]+m23*J[2])
    return B2
    
def B3_e (m31,m32,m33,J,a,b,c):
    '''
    Calculates the B3 component of the magnetic field generated by 
    n-ellipsoids in the body coordinates.
    
    Parameters:
    
    * m31,m32,m33: array
        Calculus for the ellipsoid magnetic field.
    * J: array
        Resultant magnetization vector with self-demagnetization correction.
    * a,b,c: floats
        Major, intermediate and minor axis, respectively.
    
    Returns:
    
    * B3:  array
        The B3 component of the magnetic field generated 
        by n-ellipsoids in the body coordinates.    
    '''
    
    B3 = 2*np.pi*a*b*c*(m31*J[0]+m32*J[1]+m33*J[2])
    return B3