# -*- coding: utf-8 -*-
# farfield.py

#import math
import numpy as np
#from numba import jit

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

"""
# far field (Numba version)
@jit(cache=True, nopython=True)
def field(ifreq, freq, theta, phi, ffctr, Ne, E_posc, E_lng, E_tan, Iground, Xc):

    cfar = np.zeros(2, 'c16')  # complex e-theta/e-phi
    pfar = np.zeros(7, 'f8')    # power e-abs/e-theta/e-phi/e-major/e-minor/e-RHCP/e-LHCP

    pos = np.zeros(3, 'f8')
    tan = np.zeros(3, 'f8')

    kwave = (2 * math.pi * freq) / C

    sint = math.sin(np.deg2rad(theta))
    cost = math.cos(np.deg2rad(theta))
    sinp = math.sin(np.deg2rad(phi))
    cosp = math.cos(np.deg2rad(phi))

    # グラウンド板があるときは、下半分の遠方界は0
    if Iground and (cost < -EPS):
        return cfar, pfar

    # r,θ,φ方向の単位ベクトル
    r1 = np.array([+sint * cosp, +sint * sinp, +cost])
    t1 = np.array([+cost * cosp, +cost * sinp, -sint])
    p1 = np.array([       -sinp,        +cosp,     0])

    # 線分要素の和
    sumt = 0
    sump = 0
    for ie in range(Ne):
        c = Xc[ifreq, ie]
        lng = E_lng[ie, 0]

        for k in range(3):
            pos[k] = E_posc[ie, k]
            tan[k] = E_tan[ie, k]

        ctheta, cphi = _farcomp(c, pos, tan, lng, kwave, r1, t1, p1)
        sumt += ctheta
        sump += cphi

        # グラウンド板イメージ
        if Iground:
            pos[2] *= -1
            tan[2] *= -1

            ctheta, cphi = _farcomp(c, pos, tan, lng, kwave, r1, t1, p1)
            sumt -= ctheta
            sump -= cphi

    # 因子
    cfctr = -1j * kwave * math.sqrt(ETA0 / (8 * math.pi * ffctr))
    cfar[0] = cfctr * sumt
    cfar[1] = cfctr * sump

    # 遠方界電力
    pfar = _component(cfar[0], cfar[1])

    return cfar, pfar
"""

# far field (array version)
def field(ifreq, freq, theta, phi, ffctr, Ne, E_posc, E_lng, E_tan, Iground, Xc):

    Ne = Ne  # dummy
    cfar = np.zeros(2, complex)  # complex e-theta/e-phi
    pfar = np.zeros(7, float)    # power e-abs/e-theta/e-phi/e-major/e-minor/e-RHCP/e-LHCP

    k0 = (2 * np.pi * freq) / C

    sint = np.sin(np.deg2rad(theta))
    cost = np.cos(np.deg2rad(theta))
    sinp = np.sin(np.deg2rad(phi))
    cosp = np.cos(np.deg2rad(phi))

    # グラウンド板があるときは、下半分の遠方界は0
    if Iground and (cost < -EPS):
        return cfar, pfar

    # r,θ,φ方向の単位ベクトル
    r1 = [+ sint * cosp, + sint * sinp, + cost]
    t1 = [+ cost * cosp, + cost * sinp, - sint]
    p1 = [       - sinp,        + cosp,      0]

    # 線分要素の和
    cri  = Xc[ifreq, :]
    lng  = E_lng[:, 0]
    posx = E_posc[:, 0]
    posy = E_posc[:, 1]
    posz = E_posc[:, 2]
    tanx = E_tan[:, 0]
    tany = E_tan[:, 1]
    tanz = E_tan[:, 2]
        
    pr =  (posx * r1[0]) + (posy * r1[1]) + (posz * r1[2])
    ct = ((tanx * t1[0]) + (tany * t1[1]) + (tanz * t1[2])) * cri * lng
    cp = ((tanx * p1[0]) + (tany * p1[1]) + (tanz * p1[2])) * cri * lng
    cexp = np.exp(1j * k0 * pr)
    sumt = np.sum(ct * cexp)
    sump = np.sum(cp * cexp)

    # グラウンド板イメージ
    if Iground:
        posz = -posz
        tanz = -tanz

        pr =  (posx * r1[0]) + (posy * r1[1]) + (posz * r1[2])
        ct = ((tanx * t1[0]) + (tany * t1[1]) + (tanz * t1[2])) * cri * lng
        cp = ((tanx * p1[0]) + (tany * p1[1]) + (tanz * p1[2])) * cri * lng
        cexp = np.exp(1j * k0 * pr)
        sumt -= np.sum(ct * cexp)
        sump -= np.sum(cp * cexp)

    # far field factor (complex)
    cfctr = -1j * k0 * np.sqrt(ETA0 / (8 * np.pi * ffctr))

    # far field e-theta/e-phi (complex)
    cfar[0] = cfctr * sumt
    cfar[1] = cfctr * sump

    # power of far field
    pfar = _component(cfar[0], cfar[1])

    return cfar, pfar


# (private) far field components (power)
#@jit(cache=True, nopython=True)
def _component(eth, eph):
    p = np.zeros(7, float)
    
    # abs
    p[0] = abs(eth)**2 + abs(eph)**2

    # theta/phi
    p[1] = abs(eth)**2
    p[2] = abs(eph)**2

    # major/minor
    tmp = abs(eth**2 + eph**2)
    p[3] = (abs(eth)**2 + abs(eph)**2 + tmp) / 2
    p[4] = (abs(eth)**2 + abs(eph)**2 - tmp) / 2

    # RHCP/LHCP
    p[5] = abs(eth + (1j * eph))**2 / 2
    p[6] = abs(eth - (1j * eph))**2 / 2

    return p

# far field factor
def factor(Ne, Nfeed, E_ifeed, E_feed, Zin, Z0, ifreq, mloss):
    f = 0

    if Nfeed > 0:
        # feed
        ifeed = 0
        for ie in range(Ne):
            if E_ifeed[ie] > 0:
                zin = Zin[ifeed, ifreq]
                v = E_feed[ie, 0]
                power = 0.5 * (1 / zin).real * v**2
                # matching loss
                if mloss > 0:
                    gamma = abs((zin - Z0) / (zin + Z0))
                    power = power / (1 - gamma**2)
                f += power
                ifeed += 1
    else:
        # plane wave
        f = 1 / (2 * ETA0)

    return f
