//CL//
/*
 *  Update H-field, including any PMLs.
 *   Also precalculate values for poynting vector if necessary.
 *
 *  Template parameters:
 *   common_header: Rendered contents of common.cl
 *   pmls: [('x', 'n'), ('z', 'p'),...] list of pml axes and polarities
 *   pml_thickness: Number of cells (integer)
 *   do_poynting: Whether to precalculate poynting vector components (boolean)
 *
 *  OpenCL args:
 *   E, H, dt, [p{01}h{np}, Psi_{xyz}{np}_H], [oS]
 */

#pragma OPENCL EXTENSION cl_khr_3d_image_writes : enable
#define read_field(E, x, y, z) read_imagef(E, sampler, (int4)(z, y, x, 0))
#define write_field(E, x, y, z, E__) write_imagef(E, (int4)(z, y, x, 0), E__)


const sampler_t sampler = CLK_FILTER_NEAREST | \
                          CLK_NORMALIZED_COORDS_FALSE | \
                          CLK_ADDRESS_NONE;


__kernel void update_h(
        __read_only image3d_t E,
        __global float* H,
        const float dt,
        const size_t n
        {% if pmls -%}
        __global float* S,
        {%- endif %}
        ) {

    size_t lid = get_local_id(0);
    size_t gsize = get_global_size(0);
    size_t work_group_start = get_local_size(0) * get_group_id(0);
    size_t i;

    for (i = work_group_start + lid; i < n; i += gsize) {

        {{common_header | indent(8, False)}}

        ////////////////////////////////////////////////////////////////////////////

        float4 E__ = read_field(E, x, y, z);
        float4 Epx = read_field(E, x + px, y, z);
        float4 Epy = read_field(E, x, y + py, z);
        float4 Epz = read_field(E, x, y, z + pz);

        {% if pmls -%}
        const int pml_thickness = {{pml_thickness}};
        {%- endif %}

        /*
         *  Precalculate derivatives
         */
        float dExy = Epy.x - E__.x;
        float dExz = Epz.x - E__.x;

        float dEyx = Epx.y - E__.y;
        float dEyz = Epz.y - E__.y;

        float dEzx = Epx.z - E__.z;
        float dEzy = Epy.z - E__.z;


        /*
         *  PML Update
         */
        // PML contributions to H
        float pHxi = 0;
        float pHyi = 0;
        float pHzi = 0;

        {%- for r, p in pmls -%}
            {%- set u, v = ['x', 'y', 'z'] | reject('equalto', r) -%}
            {%- set psi = 'Psi_' ~ r ~ p ~ '_H' -%}
            {%- if r != 'y' -%}
                {%- set se, sh = '-', '+' -%}
            {%- else -%}
                {%- set se, sh = '+', '-' -%}
            {%- endif -%}

            {%- if p == 'n' %}

        if ( {{r}} < pml_thickness ) {
            const size_t ir = {{r}};                          // index into pml parameters

            {%- elif p == 'p' %}

        if ( s{{r}} > {{r}} && {{r}} >= s{{r}} - pml_thickness ) {
            const size_t ir = (s{{r}} - 1) - {{r}};                     // index into pml parameters

            {%- endif %}
            const size_t ip = {{v}} + {{u}} * s{{v}} + ir * s{{v}} * s{{u}};  // linear index into Psi
            {{psi ~ u}}[ip] = p0h{{p}}[ir] * {{psi ~ u}}[ip] + p1h{{p}}[ir] * dE{{v ~ r}};
            {{psi ~ v}}[ip] = p0h{{p}}[ir] * {{psi ~ v}}[ip] + p1h{{p}}[ir] * dE{{u ~ r}};
            pH{{u}}i {{sh}}= {{psi ~ u}}[ip];
            pH{{v}}i {{se}}= {{psi ~ v}}[ip];
        }
        {%- endfor %}

        /*
         *  Update H
         */
        {% if do_poynting -%}
        // Save old H for averaging
        float Hx_old = Hx[i];
        float Hy_old = Hy[i];
        float Hz_old = Hz[i];
        {%- endif %}

        // H update equations
        float4 H__ = vload4(i, H);
        H__.x -= dt * (dEzy - dEyz + pHxi);
        H__.y -= dt * (dExz - dEzx + pHyi);
        H__.z -= dt * (dEyx - dExy + pHzi);
        vstore4(H__, i, H);

        {% if do_poynting -%}
        // Average H across timesteps
        float aHxt = Hx[i] + Hx_old;
        float aHyt = Hy[i] + Hy_old;
        float aHzt = Hz[i] + Hz_old;

        /*
         *  Calculate unscaled S components at H locations
         */
        float4 S__;
        S__.x = Epx.y * aHzt - Epx.z * aHyt;
        S__.y = Epy.z * aHxt - Epy.x * aHzt;
        S__.z = Epz.x * aHyt - Epz.y * aHxt;

        write_field(S, x, y, z, S__);
        {%- endif -%}
    }
}
