import Mesh.MeshGen as MeshGen
import Mesh.View as View
import numpy as np
from Greens import EMConst
from pymiecoated import Mie
from Greens.Greens import *

# some useful constants
eps0 = EMConst.eps0
mu0 = EMConst.mu0
c0 = EMConst.c0

def Ein_x(omega, xyz_sphere):
    k0 = omega*np.sqrt(mu0 * eps0)
    N_p = xyz_sphere.shape[0]
    Ein = np.zeros((N_p*3,1), complex)  # the incident field
    for n in range(N_p):
        Ein[n] =  np.exp(-1j*k0*xyz_sphere[n,0])
    return Ein


def GMatGen(omega, dl_new, xyz_sphere, order):
    N_p = xyz_sphere.shape[0] #the number of points
    sample, weight = np.polynomial.legendre.leggauss(order)
    Gxx = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gyx = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gzx = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gxy = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gyy = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gzy = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gxz = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gyz = np.asmatrix(np.zeros((N_p, N_p), complex))
    Gzz = np.asmatrix(np.zeros((N_p, N_p), complex))


    # this is the matrix shows the field generated by x direction of current
    direction = np.array([1,0,0])
    for m in range(N_p):
        for n in range(N_p):
            if m == n:
                Gxx[m,n] = LocalE(omega, dl_new, sample, weight)
                Gyx[m,n] = 0.0 + 0.0j
                Gzx[m,n] = 0.0 + 0.0j
            else:
                Gxx[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[0]
                Gyx[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[1]
                Gzx[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[2]
    Gx = np.concatenate((Gxx, Gyx, Gzx), axis=1)
    # this is the matrix shows the field generated by x direction of current
    direction = np.array([0,1,0])
    for m in range(N_p):
        for n in range(N_p):
            if m == n:
                Gyy[m,n] = LocalE(omega, dl_new, sample, weight)
                Gxy[m,n] = 0.0 + 0.0j
                Gzy[m,n] = 0.0 + 0.0j
            else:
                Gxy[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[0]
                Gyy[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[1]
                Gzy[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[2]
    Gy = np.concatenate((Gxy, Gyy, Gzy), axis=1)


    # this is the matrix shows the field generated by x direction of current
    direction = np.array([0,0,1])
    for m in range(N_p):
        for n in range(N_p):
            if m == n:
                Gzz[m,n] = LocalE(omega, dl_new, sample, weight)
                Gyz[m,n] = 0.0 + 0.0j
                Gxz[m,n] = 0.0 + 0.0j
            else:
                Gxz[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[0]
                Gyz[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[1]
                Gzz[m,n] = RemoteE(omega, xyz_sphere[m, :], xyz_sphere[n, :],
                                   dl_new, sample, weight, direction)[2]
    Gx = np.concatenate((Gxx,Gxy,Gxz), axis = 1)
    Gy = np.concatenate((Gyx,Gyy,Gyz), axis = 1)
    Gz = np.concatenate((Gzx,Gzy,Gzz), axis = 1)
    G = np.vstack((Gx, Gy, Gz))
    return G

if __name__ == '__main__':

    order =5
    freq = 100e6
    radii = 0.10
    eps_r = 2 - 10j
    delta_eps = (eps_r - 1)*eps0
    omega = freq * 2 * np.pi
    sigma = -np.imag(eps_r*eps0) * omega


    c = c0 / np.sqrt(abs(eps_r))
    l_lam = c / freq   # the wave length in object
    dl = l_lam / 45

    def Power_abs_mie(radii, freq, eps_r):
        l_lam_0 = c0 / freq
        size = 2 * PI * radii / l_lam_0
        mie = Mie(x = size, eps = np.conjugate(eps_r))
        Qabs = mie.qabs()
        Power = 0.5 / np.sqrt(mu0 / eps0)
        Power_abs = Power * Qabs * PI * radii**2
        return Power_abs



    #mesh generation
    xyz_sphere, dl_new = MeshGen.SphereMeshGen(radii, dl)
    N_p = xyz_sphere.shape[0]
    dV = dl_new ** 3
    size = 10  # the dot size
    View.ViewPoints(xyz_sphere, size)

    #Generate the G matrix
    GJ = GMatGen(omega, dl_new, xyz_sphere, order)
    G_E = 1j*omega*delta_eps*GJ - np.identity(3*N_p, complex)
    Ein = -Ein_x(omega, xyz_sphere)
    E = np.linalg.solve(G_E, Ein)
    Power_abs = 0.5 * sigma * np.sum(np.abs(E)**2)*dV



    Volume = 4.0 / 3.0 * np.pi * radii**3
    Volume_num = dV * N_p
    Power_abs_correct = Power_abs * ((Volume / N_p) / dV)
    #print Volume
    #print Volume_num
    print Power_abs
    print Power_abs_correct
    print Power_abs_mie(radii, freq, eps_r)
    #print radii, freq, eps_r
