# -*- coding: utf-8 -*-
"""
nearfield.py
near field
"""

import numpy as np
import sol.planewave

C = 2.99792458e8
ETA0 = C * 4 * np.pi * 1e-7
EPS = 1e-10

def field(Nfeed, Iplanewave, Planewave, E_posc, E_posm, E_posp, E_tan, E_lng, Iground, Xc, pos, ifreq, k0, noinc):
    # E [V/m], H [A/m], angle [degree]
    eh = np.zeros(14, float)

    if (Iground > 0) and ((k0 * pos[2]) < -EPS):
        return

    eposc = E_posc.copy()
    eposm = E_posm.copy()
    eposp = E_posp.copy()
    etan  = E_tan.copy()
    elng  = E_lng[:, 0].copy()

    e, h = _fcurrent(pos, Xc[ifreq], eposc, eposm, eposp, etan, elng, k0)

    # add ground image
    if Iground > 0:
        eposc[:, 2] *= -1
        eposm[:, 2] *= -1
        eposp[:, 2] *= -1
        etan[:, 2]  *= -1
        e1, h1 = _fcurrent(pos, Xc[ifreq], eposc, eposm, eposp, etan, elng, k0)
        e -= e1
        h -= h1

    # add plane wave incidence
    if (Nfeed == 0) and (noinc == 0):
        pol   = Iplanewave
        theta = Planewave[0]
        phi   = Planewave[1]
        a     = Planewave[2]
        r     = Planewave[3]
        e2, h2 = sol.planewave.planewave(pos, theta, phi, pol, a, r, k0, Iground)
        e += e2
        h += h2 / ETA0

    # component
    eh[ 0] = np.sqrt(abs(e[0])**2 + abs(e[1])**2 + abs(e[2])**2)
    eh[ 1] = abs(e[0])
    eh[ 2] = abs(e[1])
    eh[ 3] = abs(e[2])
    eh[ 4] = np.angle(e[0], True)
    eh[ 5] = np.angle(e[1], True)
    eh[ 6] = np.angle(e[2], True)
    eh[ 7] = np.sqrt(abs(h[0])**2 + abs(h[1])**2 + abs(h[2])**2)
    eh[ 8] = abs(h[0])
    eh[ 9] = abs(h[1])
    eh[10] = abs(h[2])
    eh[11] = np.angle(h[0], True)
    eh[12] = np.angle(h[1], True)
    eh[13] = np.angle(h[2], True)
    
    return eh

# array version (private)
# reshape : for broadcast !!
def _fcurrent(pos, xc, eposc, eposm, eposp, etan, elng, k0):
    xc = xc.reshape((-1, 1)) # for broadcast !!

    r0 = 0.25 * elng  # very near
    kl = (k0 * elng).reshape((-1, 1))

    dc = pos - eposc
    dm = pos - eposm
    dp = pos - eposp

    ldc = np.sqrt(dc[:, 0]**2 + dc[:, 1]**2 + dc[:, 2]**2)
    ldm = np.sqrt(dm[:, 0]**2 + dm[:, 1]**2 + dm[:, 2]**2)
    ldp = np.sqrt(dp[:, 0]**2 + dp[:, 1]**2 + dp[:, 2]**2)

    rc = np.maximum(ldc, r0).reshape((-1, 1))
    rm = np.maximum(ldm, r0).reshape((-1, 1))
    rp = np.maximum(ldp, r0).reshape((-1, 1))

    r1c = dc / rc
    r1m = dm / rm
    r1p = dp / rp

    krc = k0 * rc
    krm = k0 * rm
    krp = k0 * rp

    gc = (np.exp(-1j * krc) / krc).reshape((-1, 1))
    gm = (np.exp(-1j * krm) / krm).reshape((-1, 1))
    gp = (np.exp(-1j * krp) / krp).reshape((-1, 1))

    fc = ((1 / krc) + 1j).reshape((-1, 1))
    fm = ((1 / krm) + 1j).reshape((-1, 1))
    fp = ((1 / krp) + 1j).reshape((-1, 1))

    efctr = (1j * k0 * ETA0) / (4 * np.pi)
    hfctr = - k0 / (4 * np.pi)

    e1 = efctr * xc * (- (r1p * fp * gp) + (r1m * fm * gm) - (etan * kl * gc))
    h1 = hfctr * xc * np.cross(r1c, etan) * fc * kl * gc

    # e = E [V/m], h = H [A/m] (3 vector)
    e = np.sum(e1, axis=0)
    h = np.sum(h1, axis=0)
    
    return e, h
