#****************************************************************
#NAME:  ReflectedField
#TYPE:	Subroutine
#DESCRIPTION: The subroutine implements the formula (13-20) of the 
#book C. Balanis "Advanced Engineering Electromagnetics" for the reflected
#field in the GO approach. Here the assumptions are: flat and PEC
#surface, spherical spreading factor. 
#INPUT: EiS= complex vector cointaing the field from the source at the reflection point
#	flagS= integer to account for 1=the source is the actual source 0=the source is one of the reflections
#	flagP= integer to account for 1=the observation is the final point 0=the observation is one of the further reflections
#	mu= real vector of the dipole moment
#	norm= real vector of the normal to the surface at the reflection point
#	S= real vector of coordinates of the source
#	QR= real vector of coordinates of the reflection point
#	f= frequency (Hz)
#
#OUTPUT: Eref= complex vector containing the reflected field at the 
#				observation point P
#
#CREATED BY: Giacomo Donzelli
#DATE:			13/04/2010
#MODIFICATIONS:
#COMMENTS: - for now, we assume that the first iteration comes out from a DEC with moment mu, thus flagS=1, then
#			flagS=0 and the field from the source is the field on QR from the previous itaration
#		       - if flagP=1 the observation point is the final point, otherwise it will be iterated
#****************************************************************
import numpy as np
from math import *
from cmath import *

        
class ReflectedField(object):
    """Classe che valuta il campo riflesso"""
    
    #Input pubblici
    g_freq  = None #freq.
    g_flags = None
    g_flagp = None 
    g_norm  = np.zeros(3, dtype=float) #normale (x,y,z)
    g_S     = np.zeros(3, dtype=float) #sorgente(x,y,z)
    g_QR    = np.zeros(3, dtype=float) #punto di riflessione(x,y,z)
    g_P     = np.zeros(3, dtype=float) #punto osservazione(x,y,z)
    g_mu    = np.zeros(3, dtype=float) #verso orientato come il dipolo.
    g_Eis   = np.zeros(3, dtype=complex) 
    g_w     = None
    g_k     = None
    g_lambda = None
    #PEC
    g_ReflMatrix = np.zeros([2,2], float)
    g_ReflMatrix[0][0]= 1
    g_ReflMatrix[1][1]= -1
    eps_0      = 8.85 * 10**(-12)
    mu_0       = 4* pi + 10**(-7)
        
    #-Init degli Input-
    def __init__(self, r_flags, r_flagp, r_freq, r_norm, r_s, r_qr, r_p, r_mu, r_eis):
        self.g_flags	= r_flags   #0 always. 1 se sulla sorgente
        self.g_flagp	= r_flagp   # 0 always. 1 se sul punto di osservazione
        self.g_freq 	= r_freq    #(Hz)
        self.g_norm 	= r_norm
        self.g_S 	= r_s
        self.g_QR 	= r_qr
        self.g_P 	= r_p
        self.g_mu 	= r_mu
        self.g_Eis	= r_eis
        self.g_w 	= 2*pi*r_freq
        self.g_k 	= self.g_w/(3.0*10**8)
        self.g_lambda 	= 2.0*pi/self.g_k
        self.eps_0      = 8.85 * 10**(-12)
        self.mu_0       = 4* pi + 10**(-7)
        
  
    def EvalEis(self):
        """Incident Field: Metodo che restituisce Eis"""		    

        ri = sqrt((self.g_QR[0]-self.g_S[0])**2+(self.g_QR[1]-self.g_S[1])**2+(self.g_QR[2]-self.g_S[2])**2)
        uri = np.zeros(3, dtype=float)
        uri[0] = (self.g_QR[0]-self.g_S[0])/ri
        uri[1] = (self.g_QR[1]-self.g_S[1])/ri
        uri[2] = (self.g_QR[2]-self.g_S[2])/ri
        
        uprpi  = np.cross(self.g_norm, uri)
        uprpi  = uprpi/abs(uprpi)
        upari  = np.cross(uprpi, uri)
        
        if self.g_flags == 1: #Eis ritorna modificato solo nel caso in cuiflags = 1
            temp  = np.cross(self.g_mu, uri)
            temp  = temp//abs(temp)
            pol  = np.cross(temp, uri)
            Eis = exp(-1j*k*ri)/(ri*pol)				
            return Eis
        else:
            return g_Eis #return Eis not modified
    
    def EvalEref(self):
        """reflected Field: Metodo che restituisce Eref"""
    
        EiQR = np.zeros(2, dtype=complex) 
        Eiloc= np.zeros(2, dtype=complex)			
        uri[0] = (self.g_QR[0]-self.g_P[0])/ri
        uri[1] = (self.g_QR[1]-self.g_P[1])/ri
        uri[2] = (self.g_QR[2]-self.g_P[2])/ri
        
        uprpi  = np.cross(self.g_norm, uri)
        uprpi  = uprpi/abs(uprpi)
        upari  = np.cross(uprpi, uri)	
        
        r  = sqrt((self.g_QR[0]-self.g_P[0])**2+(self.g_QR[1]-self.g_P[1])**2+(self.g_QR[2]-self.g_P[2])**2)
        ri = sqrt((self.g_QR[0]-self.g_S[0])**2+(self.g_QR[1]-self.g_S[1])**2+(self.g_QR[2]-self.g_S[2])**2)
        
        if (abs(r/self.glambda)<1.0*10**(-9)):#da definire questa eccezione
            ur = self.g_norm
            uprp = uprpi
            upar = upari
        else:
            ur[0] = (P[0]-QR[0])/r
            ur[1] = (P[1]-QR[1])/r
            ur[2] = (P[2]-QR[2])/r
            uprp  = np.cross(self.g_norm, ur)
            uprp  = uprp//abs(uprp)
            upar  = np.cross(uprp, ur)
        
        EiQR[0]=self.g_Eis*upari
        EiQR[1]=self.g_Eis*uprpi

        if self.g_flagp == 1:
            Eloc = np.matrixmultiply(self.g_ReflMatrix,EiQR)*(ri/(r+ri))*exp(-1j*k*r)
        else:			
            Eloc = matrixmultiply(self.g_ReflMatrix,EiQR)
        
        Eref = Eloc[0]*upar+Eloc[1]*uprp
    
        return Eref 

