Corrections needed for compute_power_density function

[dreme]

The compute_power_density function needs some corrections:

Firstly, the field level input should be specified as S0 rather than E0 and S0 does not need to be included in the function outputs.

Secondly, the function should accept a numpy array input for z, rather than assuming z = np.linspace(-2, 0, 2001)

Thirdly, the SE and SE outputs are not aligning with the FEKO results - e.g. see below:

For another comparator, the planeEMRyong function below can be used as it has been validated against the FEKO results:

def planeEMRyong(E0,z,zi,εr,μr,σ,f,θ,pol):
    '''Calculates power density (SE = ½|E|²/Z0, SH = ½|H|²Z0) for a plane wave
       travelling through a multilayered infinite planar medium.
       The first and last layers have infinite thickness.
       Theory on which the algorithm is based can be found e.g. in Weng
       Cho Chew, "Waves and Fields in Inhomogeneous Media", IEEE PRESS
       Series on Electromagnetic Waves. This code was adapted from
       a Matlab script developed by Kimmo Karkkainen.

       INPUTS:
        E0 = E-field level of incident plane wave (V/m peak)
         z = list of z coords where |E|, SAR field levels are evaluated (m)
        zi = list of z coords for interfaces between planar layers (m)
        εr = list of relative permittivity for each layer
        μr = list of relative permeability for each layer
         σ = list of conductivity (S/m) for each layer
         f = frequency (Hz) of the planewave
         θ = incoming angle of the propagation vector (k)
             relative to the normal of the first interface. 
       pol = polarization of the plane wave:
               'TM' = TM polarized wave (H parallel to interface)
               'TE' = TE polarized wave (E parallel to interface)

       OUTPUTS:       
          SE = Equivalent plane wave power density (SE = ½|E|²/Z0) at specified points (z)
          SH = Equivalent plane wave power density (SE = ½Z0|H|²) at specified points (z)
'''

    # Set constants
    Z0 = 376.730313668           # free space impedance
    ε0 = 8.854187817620389e-12   # permittivity of free space
    µ0 = 1.2566370614359173e-06  # permeability of free space
    π = np.pi

    # Initialise variables    
    z = np.array(z).round(8)
    θ = θ * π / 180.   # convert θ from degrees to radians
    zi.append(1e9)     # add a very large z coord for 'infinite' last layer
    N = len(εr)        # N is the no. of layers
    w = 2. * π * f     # angular frequency 
    i = 1j             # i is imaginary unit

    # Wavenumber and its z-component in layers
    ε = [er*ε0 + s/(i*w) for er, s in zip(εr,σ)]
    μ = [ur*µ0 for ur in μr]
    K = [w*np.sqrt(e*µ0) for e in ε]
    Kz = [0]*N
    Kz[0] = K[0] * np.cos(θ)
    Kt    = K[0] * np.sin(θ)

    for L in range(1, N):
        SQ = K[L]**2 - Kt**2
        SQr = np.real(SQ)
        SQi = np.imag(SQ)
        if (SQi == 0) & (SQr < 0):
            Kz[L] = -np.sqrt(SQ)
        else:
            Kz[L] = np.sqrt(SQ)

    # Calculate reflection and transmission coefficients for interfaces
    R = np.zeros((N,N), dtype=complex)
    T = R.copy()
    for k in range(N-1):
        if pol == 'TM':
            # Reflection coefficient for magnetic field
            e1, e2 = ε[k+1] * Kz[k], ε[k] * Kz[k+1]
            R[k,k+1] = (e1 - e2) / (e1 + e2)
        elif pol == 'TE':
            # Reflection coefficient for electric field
            m1, m2 = μ[k+1] * Kz[k],  μ[k] * Kz[k+1]
            R[k,k+1] = (m1 - m2) / (m1 + m2)
        else:
            raise Exception(f'pol ({pol}) must be either TE or TM')

        R[k+1,k] = -R[k,k+1] 
        T[k,k+1] = 1 + R[k,k+1]
        T[k+1,k] = 1 + R[k+1,k]

    # Calculate generalized reflection coefficients for interfaces:
    gR = np.zeros(N, dtype=complex)
    for k in range(N-1)[::-1]:
        thickness = zi[k+1] - zi[k]  # layer thickness
        g = gR[k+1] * np.exp(-2*i*Kz[k+1] * thickness)
        gR[k] = (R[k,k+1] + g) / (1 + R[k,k+1] * g)

    A = np.zeros(N, dtype=complex)
    if pol == 'TM':
        A[0] = np.abs(K[0] / (w*μ[0]))  # Amplitude of magnetic field in layer 1
    else:
        A[0] = 1                        # Amplitude of electric field in layer 1

    # Calculate amplitudes in other layers:
    for k in range(1, N):
        A[k] = A[k-1] * np.exp(-i*Kz[k-1]*zi[k-1])*T[k-1,k] \
            /(1-R[k,k-1] * gR[k] * np.exp(-2*i*Kz[k]*(zi[k]-zi[k-1]))) \
            /np.exp(-i*Kz[k]*zi[k-1])

    # Form a vector that tells in which layer the calculation points are located:
    zl = []
    layer = 0
    for zp in z:
        while zp >= zi[layer]:
            layer += 1
        zl.append(layer)

    # Calculate E-field:
    Azl = np.array([A[z] for z in zl])
    Kzzl = np.array([Kz[z] for z in zl])
    Kzl = np.array([K[z] for z in zl])
    gRzl = np.array([gR[z] for z in zl])
    εzl = np.array([ε[z] for z in zl])
    zizl = np.array([zi[z] for z in zl])
    σzl = np.array([σ[z] for z in zl])

    if pol == 'TM':
        # The forward and backward Hf and Hb have only y component

        Hf = Azl * np.exp(-i*Kzzl*z)
        Hb = Azl * gRzl * np.exp(-2*i*Kzzl*zizl + i*Kzzl*z)
        Ex = Z0 * np.cos(θ) * (Hf-Hb)
        Ez = (-Z0) * np.sin(θ) * (Hf+Hb)
        Sh = 0.5 * E0 ** 2 * Z0 * (np.abs(Hf+Hb))**2
        Se = 0.5 * E0 ** 2 * (1/Z0) * (np.abs(Ex)**2 + np.abs(Ez)**2)

    else:
        # The forward and backward Ef and Eb have only y component
        Ef = Azl * np.exp(-i*Kzzl*z)
        Eb = Azl * gRzl * np.exp(-2*i*Kzzl*zizl + i*Kzzl*z)
        Hx = -(1/Z0) * np.cos(θ) * (Ef-Eb)
        Hz = (1/Z0) * np.sin(θ) * (Ef+Eb) 
        Sh = 0.5 * E0 ** 2 * Z0 * (np.abs(Hx)**2 + np.abs(Hz)**2)
        Se = 0.5 * E0 ** 2 * (1/Z0) * (np.abs(Ef+Eb))**2

    return Se, Sh

[yc1976Python]

@dreme, can you please check and specify the incident angle values for the plot for SE?

[dreme]

@yc1976Python, see updated plots with the correct θ labels. The line plots remain the same though.

[yc1976Python]

Firstly, the field level input should be specified as S0 rather than E0 and S0 does not need to be included in the function outputs. Yong: OK. In the next code update, we can use S0 as the direct input and agree that no need to plot this in the outputs.

Secondly, the function should accept a numpy array input for z, rather than assuming z = np.linspace(-2, 0, 2001) Yong: Yes, this can be done with no problem. We can specify the number of points as per the formula discussed. We can either pass Nt as an input or pass the whole z array z = np.linspace(-2, 0, Nt) as an input for the function.

Thirdly, the SE and SE outputs are not aligning with the FEKO results - e.g. see below: Yong: Hi Vitas, the code in the main branch is the latest one and gives the correct result. The result can be plot from 0 to +2m. The code you attached is before the refactoring. One important point is that the array needs to be constructed as z = np.linspace(-2, 0, Nt) to be used for the function. Otherwise, it will give us wrong result

Below is the new with z arranged from 0 to 2m (f=1000MHz, at 0 and 89 degrees of incidence angles respectively).

[dreme]

@yc1976Python I would still like to see S0 in the plots but there is no need to include S0 as an output of the compute_power_density function since it is already known as an input to the function.

Probably better to pass the whole z array rather than Nt since this gives more flexibility (e.g. specifying different height range.

I'll download your latest code for the compute_power_density function from the main branch and check again.

[dreme]

@yc1976Python I just tried the compute_power_density code from the main branch and it is still giving the same results as before. However, the plots that you displayed above look fine. Hmmm 😕

Here is the code that I'm running for the comparison of the compute_power_density function to the FEKO results:

def compareScompute(thetas,gnd,er,σ,fMHz,pol,Stype):
    '''Compare FEKO, planeEMR & planeEMRyong results for two layer case'''

    E0 = 100
    z = np.linspace(-2,0,2001)
    fig, axes = plt.subplots(1,5, figsize=(10,6))

    for θ, ax in zip(thetas,axes):

        # Calculate compute_power_density SE & SH levels
        Se, Sh, S0 = compute_power_density(gnd, E0, fMHz, θ, pol)
        Scalc = Se if Stype=='SE' else Sh

       # Get FEKO S levels from Sm dataframe
        mask = (Sm.gnd==gnd) & (Sm.pol==pol) & (Sm.theta==θ) & \
               (Sm.fMHz==fMHz) & (Sm.z <= 2) & (Sm['S type']==Stype)
        df = Sm[mask]

        # Create plot
        ax.plot(df.S,df.z,'-',lw=4, alpha=0.7, label='FEKO',color='pink')
        ax.plot(Scalc,-z,'--',lw=1,label='calculated',color='darkblue')

        if θ==0: ax.set_ylabel('z (m)')
        ax.set_xlabel('S (W/m²)')
        ax.set_xlim(0, 60)
        ax.set_ylim(0, 2.4)
        ax.legend(fontsize=8)
        ax.grid(ls='--')
        ax.set_title(r'$\theta=$'+f'{θ}°')

    plt.suptitle(f'Compare FEKO to compute_power_density\n{fMHz:g} MHz, {gnd} ground, {Stype}, {pol}')    
    plt.tight_layout()
    plt.savefig(r'C:\Users\emeau\OneDrive\Desktop\compareS.svg')

Here also is the compute_power_density function that I copied from main

import numpy as np
from scipy.constants import epsilon_0 as eps0, mu_0 as mu0

__all__ = ["compute_power_density"]

def compute_power_density(ground_type, E0, f, theta, pol):
    """Returns power density for a plane wave traveling through a
    multilayered infinite planar medium where the first and last layers
    have infinite thickness.

    Parameters
    ----------
    ground_type : string
        Type of ground, either 'PEC Ground' or any other string for
        real ground.
    E0 : float
        Electric field magnitude of the incident plane wave (V/m).
    f : float
        Frequency (Hz) of the plane wave.
    theta : float
        Incoming angle of the propagation vector (k) relative to the
        normal of the first interface.
    pol : {'TM', 'TE'}
        Polarization of the plane wave.
        'TM' = TM polarized wave (H parallel to interface)
        'TE' = TE polarized wave (E parallel to interface)

    Returns
    -------
    tuple
        Containing the equivalent plane wave power density (SH, SE, S0)

    Notes
    -----
    The algorithm is based on the theory found in Weng Cho Chew,
    'Waves and Fields in Inhomogeneous Media,' IEEE PRESS Series on
    Electromagnetic Waves. The code was adapted from a Matlab script
    developed by Kimmo Karkkainen.
    """
    # initialize settings
    Z0 = np.sqrt(mu0 / eps0)  # free-space impedance
    z = np.linspace(-2, 0, 2001)  # z-direction coordinates
    zi = [0]  # interface level between layer 1 and 2
    epsr = [1, 1]  # relative permittivities of layers 1 and 2
    mur = [1, 1]  # relative permeability of layers 1 and 2
    # sigma = [0, 1e6]  # lossless conditions
    zi.append(1e9)  # add a very large z value to act as infinity
    N = len(epsr)  # number of layers
    w = 2.0 * np.pi * f * 1e6  # angular frequency
    theta = np.deg2rad(theta)

    S0 = (0.5*E0**2/Z0) * np.ones(len(z))

    if ground_type == 'PEC Ground':
        epsr = [1, 10]
        sigma = [0, 1e6]  # PEC Ground, lossless
    else:
        epsr = [1, 10]
        sigma = [0, 0.1]  # Real Ground

    # wavenumber
    eps = [er * eps0 + s / (1j * w) for er, s in zip(epsr, sigma)]
    mu = [ur * mu0 for ur in mur]
    K = [w * np.sqrt(e * mu0) for e in eps]
    Kz = [0] * N
    Kz[0] = K[0] * np.cos(theta)
    Kt = K[0] * np.sin(theta)
    for L in range(1, N):
        SQ = K[L] ** 2 - Kt**2
        SQr = np.real(SQ)
        SQi = np.imag(SQ)
        if (SQi == 0) & (SQr < 0):
            Kz[L] = -np.sqrt(SQ)
        else:
            Kz[L] = np.sqrt(SQ)

    # reflection and transmission coefficients for interfaces
    R = np.empty((N, N), dtype=complex)
    T = R.copy()
    for k in range(N - 1):
        if pol == "TM":  # reflection coefficient for the H field
            e1, e2 = eps[k + 1] * Kz[k], eps[k] * Kz[k + 1]
            R[k, k + 1] = (e1 - e2) / (e1 + e2)
        elif pol == "TE":  # reflection coefficient for the E field
            m1, m2 = mu[k + 1] * Kz[k], mu[k] * Kz[k + 1]
            R[k, k + 1] = (m1 - m2) / (m1 + m2)
        else:
            raise Exception(f"pol ({pol}) must be either TE or TM")
        R[k + 1, k] = -R[k, k + 1]
        T[k, k + 1] = 1 + R[k, k + 1]
        T[k + 1, k] = 1 + R[k + 1, k]

    # generalized reflection coefficients for interfaces
    gR = np.zeros(N, dtype=complex)
    for k in range(N - 1)[::-1]:
        thickness = zi[k + 1] - zi[k]  # layer thickness
        g = gR[k + 1] * np.exp(-2 * 1j * Kz[k + 1] * thickness)
        gR[k] = (R[k, k + 1] + g) / (1 + R[k, k + 1] * g)

    # amplitudes
    A = np.zeros(N, dtype=complex)
    if pol == "TM":  # amplitude of the H field in layer 1
        A[0] = np.abs(K[0] / (w * mu[0]))
    else:  # amplitude of the E field in layer 1
        A[0] = 1
    for k in range(1, N):  # amplitudes in other layers
        A[k] = (
            A[k - 1]
            * np.exp(-1j * Kz[k - 1] * zi[k - 1])
            * T[k - 1, k]
            / (1 - R[k, k - 1] * gR[k]
               * np.exp(-2 * 1j * Kz[k] * (zi[k] - zi[k - 1])))
            / np.exp(-1j * Kz[k] * zi[k - 1])
        )

    # form a vector that tells in which layer the calculation points are located
    zl = []
    layer = 0
    for zp in z:
        while zp >= zi[layer]:
            layer += 1
        zl.append(layer)

    # E-field and H-field components
    Azl = np.array([A[z] for z in zl])
    Kzzl = np.array([Kz[z] for z in zl])
    Kzl = np.array([K[z] for z in zl])
    gRzl = np.array([gR[z] for z in zl])
    epszl = np.array([eps[z] for z in zl])
    zizl = np.array([zi[z] for z in zl])
    sigmazl = np.array([sigma[z] for z in zl])

    if pol == "TM":  # the forward and backward Hf and Hb have only y component
        Hf = Azl * np.exp(-1j * Kzzl * z)
        Hb = Azl * gRzl * np.exp(-2 * 1j * Kzzl * zizl + 1j * Kzzl * z)
        Ex = Z0 * np.cos(theta) * (Hf - Hb)
        Ez = (-Z0) * np.sin(theta) * (Hf + Hb)
        SH = 0.5 * E0**2 * Z0 * (np.abs(Hf + Hb)) ** 2
        SE = 0.5 * E0**2 * (1 / Z0) * (np.abs(Ex) ** 2 + np.abs(Ez) ** 2)
    else:  # the forward and backward Ef and Eb have only y component
        Ef = Azl * np.exp(-1j * Kzzl * z)
        Eb = Azl * gRzl * np.exp(-2 * 1j * Kzzl * zizl + 1j * Kzzl * z)
        Hx = -(1 / Z0) * np.cos(theta) * (Ef - Eb)
        Hz = (1 / Z0) * np.sin(theta) * (Ef + Eb)
        SH = 0.5 * E0**2 * Z0 * (np.abs(Hx) ** 2 + np.abs(Hz) ** 2)
        SE = 0.5 * E0**2 * (1 / Z0) * (np.abs(Ef + Eb)) ** 2
    return SH, SE, S0

[yc1976Python]

@dreme Hi Vitas, your three comments are addressed and please have a look and let me know if you want me to pull a request.