# -*- coding: utf-8 -*-
# planewave.py
# plane wave incidence

import numpy as np

def planewave(pos, theta, phi, pol, a, r, k0, gnd):
    # ei = E [V/m], hi = Z0 * H [V/m]
    # ei, hi, pos : 3 vector
    ei = np.zeros(3, complex)
    hi = np.zeros(3, complex)

    if (gnd < 0) or (gnd > 1):
        return ei, hi
    if (pol < 1) or (pol > 5):
        return ei, hi
    if (gnd == 1) and (k0 * pos[2] < -1e-8):
        return ei, hi

    # unit vector in r, theta and phi
    sint = np.sin(np.deg2rad(theta))
    cost = np.cos(np.deg2rad(theta))
    sinp = np.sin(np.deg2rad(phi))
    cosp = np.cos(np.deg2rad(phi))
    r1 = [+ sint * cosp, + sint * sinp, + cost]
    t1 = [+ cost * cosp, + cost * sinp, - sint]
    p1 = [       - sinp,        + cosp,      0]

    # E/H field
    e = _einc(t1, p1, pol, a, r) * np.exp(1j * k0 * np.dot(pos, r1))
    h = -np.cross(r1, e)
    ei = e
    hi = h

    # add ground image
    if gnd == 1:
        r1[2] *= -1
        t1[2] *= -1
        e2 = -_einc(t1, p1, pol, a, r) * np.exp(1j * k0 * np.dot(pos, r1))
        h2 = -np.cross(r1, e2)
        ei = ei + e2
        hi = hi + h2

    return ei, hi

# E-incidence (private)
def _einc(t1, p1, pol, a, r):
    # ei, t1, p1 : 3 vector
    # pol, a, r : scalar
    ei = np.zeros(3, complex)

    if   pol == 1:
        # V
        #ei = complex(-t1, 0)
        for i in range(3):
            ei[i] = complex(-t1[i], 0)
    elif pol == 2:
        # H
        #ei = complex(+p1, 0)
        for i in range(3):
            ei[i] = complex(+p1[i], 0)
    elif pol == 3:
        # RHCP
        #ei = np.sqrt(0.5) * complex(+t1, +p1)
        for i in range(3):
            ei[i] = np.sqrt(0.5) * complex(+t1[i], +p1[i])
    elif pol == 4:
        # LHCP
        #ei = np.sqrt(0.5) * complex(+t1, -p1)
        for i in range(3):
            ei[i] = np.sqrt(0.5) * complex(+t1[i], -p1[i])
    elif pol == 5:
        # elliptical
        major = + t1 * np.cos(np.deg2rad(a)) + p1 * np.sin(np.deg2rad(a))
        minor = - t1 * np.sin(np.deg2rad(a)) + p1 * np.cos(np.deg2rad(a))
        #ei = (1 / np.sqrt(1 + r**2)) * complex(major, r * minor)
        for i in range(3):
            ei[i] = (1 / np.sqrt(1 + r**2)) * complex(major[i], r * minor[i])
    #print(pol, ei)

    return ei
