/*
 *  Update E-field, including any PMLs.
 *
 *  Template parameters:
 *   common_header: Rendered contents of common.cl
 *   pmls: [{'axis': 'x', 'polarity': 'n', 'thickness': 8}, ...] list of pml dicts containing
        axes, polarities, and thicknesses.
 *
 *  OpenCL args:
 *   E, H, dt, eps, [p{01}e{np}, Psi_{xyz}{np}_E]
 */



typedef float ftype;


/*
 * Field size info
 */
const size_t sx = 395;
const size_t sy = 345;
const size_t sz = 73;
const size_t field_size = sx * sy * sz;

//Since we use i to index into Ex[], Ey[], ... rather than E[], do nothing if
// i is outside the bounds of Ex[].
if (i >= field_size) {
    PYOPENCL_ELWISE_CONTINUE;
}


/*
 * Array indexing
 */
// Given a linear index i and shape (sx, sy, sz), defines x, y, and z
//  as the 3D indices of the current element (i).
// (ie, converts linear index [i] to field indices (x, y, z)
const size_t x = i / (sz * sy);
const size_t y = (i - x * sz * sy) / sz;
const size_t z = (i - y * sz - x * sz * sy);

// Calculate linear index offsets corresponding to offsets in 3D
// (ie, if E[i] <-> E(x, y, z), then E[i + diy] <-> E(x, y + 1, z)
const size_t dix = sz * sy;
const size_t diy = sz;
const size_t diz = 1;


/*
 * Pointer math
 */
//Pointer offsets into the components of a linearized vector-field
// (eg. Hx = H + XX, where H and Hx are pointers)
const size_t XX = 0;
const size_t YY = field_size;
const size_t ZZ = field_size * 2;

//Define pointers to vector components of each field (eg. Hx = H + XX)
__global ftype *Ex = E + XX;
__global ftype *Ey = E + YY;
__global ftype *Ez = E + ZZ;

__global ftype *Hx = H + XX;
__global ftype *Hy = H + YY;
__global ftype *Hz = H + ZZ;


/*
 * Implement periodic boundary conditions
 *
 * mx ([m]inus [x]) gives the index offset of the adjacent cell in the minus-x direction.
 * In the event that we start at x == 0, we actually want to wrap around and grab the cell
 * x_{-1} == (sx - 1) instead, ie. mx = (sx - 1) * dix .
 *
 * px ([p]lus [x]) gives the index offset of the adjacent cell in the plus-x direction.
 * In the event that we start at x == (sx - 1), we actually want to wrap around and grab
 * the cell x_{+1} == 0 instead, ie. px = -(sx - 1) * dix .
 */

int mx = -dix;
int px = +dix;
int wrap_x = (sx - 1) * dix;
if ( x == 0 ) {
  mx = wrap_x;
} else if ( x == sx - 1 ) {
  px = -wrap_x;
} 

int my = -diy;
int py = +diy;
int wrap_y = (sy - 1) * diy;
if ( y == 0 ) {
  my = wrap_y;
} else if ( y == sy - 1 ) {
  py = -wrap_y;
} 

int mz = -diz;
int pz = +diz;
int wrap_z = (sz - 1) * diz;
if ( z == 0 ) {
  mz = wrap_z;
} else if ( z == sz - 1 ) {
  pz = -wrap_z;
} 


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

__global ftype *epsx = eps + XX;
__global ftype *epsy = eps + YY;
__global ftype *epsz = eps + ZZ;


/*
 *   Precalclate derivatives
 */
ftype dHxy = Hx[i] - Hx[i + my];
ftype dHxz = Hx[i] - Hx[i + mz];

ftype dHyx = Hy[i] - Hy[i + mx];
ftype dHyz = Hy[i] - Hy[i + mz];

ftype dHzx = Hz[i] - Hz[i + mx];
ftype dHzy = Hz[i] - Hz[i + my];

/*
 *   PML Update
 */
// PML effects on E
ftype pExi = 0;
ftype pEyi = 0;
ftype pEzi = 0;

int pml_xn_thickness = 8;

if ( x < pml_xn_thickness ) {
    const size_t ir = x;                          // index into pml parameters
    const size_t ip = z + y * sz + ir * sz * sy;  // linear index into Psi
    Psi_xn_Ey[ip] = px0en[ir] * Psi_xn_Ey[ip] + px1en[ir] * dHzx;
    Psi_xn_Ez[ip] = px0en[ir] * Psi_xn_Ez[ip] + px1en[ir] * dHyx;
    pEyi -= Psi_xn_Ey[ip];
    pEzi += Psi_xn_Ez[ip];
}int pml_xp_thickness = 8;

if ( sx > x && x >= sx - pml_xp_thickness ) {
    const size_t ir = (sx - 1) - x;                     // index into pml parameters
    const size_t ip = z + y * sz + ir * sz * sy;  // linear index into Psi
    Psi_xp_Ey[ip] = px0ep[ir] * Psi_xp_Ey[ip] + px1ep[ir] * dHzx;
    Psi_xp_Ez[ip] = px0ep[ir] * Psi_xp_Ez[ip] + px1ep[ir] * dHyx;
    pEyi -= Psi_xp_Ey[ip];
    pEzi += Psi_xp_Ez[ip];
}int pml_yn_thickness = 8;

if ( y < pml_yn_thickness ) {
    const size_t ir = y;                          // index into pml parameters
    const size_t ip = z + x * sz + ir * sz * sx;  // linear index into Psi
    Psi_yn_Ex[ip] = py0en[ir] * Psi_yn_Ex[ip] + py1en[ir] * dHzy;
    Psi_yn_Ez[ip] = py0en[ir] * Psi_yn_Ez[ip] + py1en[ir] * dHxy;
    pExi += Psi_yn_Ex[ip];
    pEzi -= Psi_yn_Ez[ip];
}int pml_yp_thickness = 8;

if ( sy > y && y >= sy - pml_yp_thickness ) {
    const size_t ir = (sy - 1) - y;                     // index into pml parameters
    const size_t ip = z + x * sz + ir * sz * sx;  // linear index into Psi
    Psi_yp_Ex[ip] = py0ep[ir] * Psi_yp_Ex[ip] + py1ep[ir] * dHzy;
    Psi_yp_Ez[ip] = py0ep[ir] * Psi_yp_Ez[ip] + py1ep[ir] * dHxy;
    pExi += Psi_yp_Ex[ip];
    pEzi -= Psi_yp_Ez[ip];
}int pml_zn_thickness = 8;

if ( z < pml_zn_thickness ) {
    const size_t ir = z;                          // index into pml parameters
    const size_t ip = y + x * sy + ir * sy * sx;  // linear index into Psi
    Psi_zn_Ex[ip] = pz0en[ir] * Psi_zn_Ex[ip] + pz1en[ir] * dHyz;
    Psi_zn_Ey[ip] = pz0en[ir] * Psi_zn_Ey[ip] + pz1en[ir] * dHxz;
    pExi -= Psi_zn_Ex[ip];
    pEyi += Psi_zn_Ey[ip];
}int pml_zp_thickness = 8;

if ( sz > z && z >= sz - pml_zp_thickness ) {
    const size_t ir = (sz - 1) - z;                     // index into pml parameters
    const size_t ip = y + x * sy + ir * sy * sx;  // linear index into Psi
    Psi_zp_Ex[ip] = pz0ep[ir] * Psi_zp_Ex[ip] + pz1ep[ir] * dHyz;
    Psi_zp_Ey[ip] = pz0ep[ir] * Psi_zp_Ey[ip] + pz1ep[ir] * dHxz;
    pExi -= Psi_zp_Ex[ip];
    pEyi += Psi_zp_Ey[ip];
}

/*
 *  Update E
 */
Ex[i] += dt / epsx[i] * (dHzy - dHyz + pExi);
Ey[i] += dt / epsy[i] * (dHxz - dHzx + pEyi);
Ez[i] += dt / epsz[i] * (dHyx - dHxy + pEzi);
==========================================
/*
 *  Update H-field, including any PMLs.
 *   Also precalculate values for poynting vector if necessary.
 *
 *  Template parameters:
 *   common_header: Rendered contents of common.cl
 *   pmls: [{'axis': 'x', 'polarity': 'n', 'thickness': 8}, ...] list of pml dicts containing
        axes, polarities, and thicknesses.
 *   do_poynting: Whether to precalculate poynting vector components (boolean)
 *
 *  OpenCL args:
 *   E, H, dt, [p{xyz}{01}h{np}, Psi_{xyz}{np}_H], [oS]
 */



typedef float ftype;


/*
 * Field size info
 */
const size_t sx = 395;
const size_t sy = 345;
const size_t sz = 73;
const size_t field_size = sx * sy * sz;

//Since we use i to index into Ex[], Ey[], ... rather than E[], do nothing if
// i is outside the bounds of Ex[].
if (i >= field_size) {
    PYOPENCL_ELWISE_CONTINUE;
}


/*
 * Array indexing
 */
// Given a linear index i and shape (sx, sy, sz), defines x, y, and z
//  as the 3D indices of the current element (i).
// (ie, converts linear index [i] to field indices (x, y, z)
const size_t x = i / (sz * sy);
const size_t y = (i - x * sz * sy) / sz;
const size_t z = (i - y * sz - x * sz * sy);

// Calculate linear index offsets corresponding to offsets in 3D
// (ie, if E[i] <-> E(x, y, z), then E[i + diy] <-> E(x, y + 1, z)
const size_t dix = sz * sy;
const size_t diy = sz;
const size_t diz = 1;


/*
 * Pointer math
 */
//Pointer offsets into the components of a linearized vector-field
// (eg. Hx = H + XX, where H and Hx are pointers)
const size_t XX = 0;
const size_t YY = field_size;
const size_t ZZ = field_size * 2;

//Define pointers to vector components of each field (eg. Hx = H + XX)
__global ftype *Ex = E + XX;
__global ftype *Ey = E + YY;
__global ftype *Ez = E + ZZ;

__global ftype *Hx = H + XX;
__global ftype *Hy = H + YY;
__global ftype *Hz = H + ZZ;


/*
 * Implement periodic boundary conditions
 *
 * mx ([m]inus [x]) gives the index offset of the adjacent cell in the minus-x direction.
 * In the event that we start at x == 0, we actually want to wrap around and grab the cell
 * x_{-1} == (sx - 1) instead, ie. mx = (sx - 1) * dix .
 *
 * px ([p]lus [x]) gives the index offset of the adjacent cell in the plus-x direction.
 * In the event that we start at x == (sx - 1), we actually want to wrap around and grab
 * the cell x_{+1} == 0 instead, ie. px = -(sx - 1) * dix .
 */

int mx = -dix;
int px = +dix;
int wrap_x = (sx - 1) * dix;
if ( x == 0 ) {
  mx = wrap_x;
} else if ( x == sx - 1 ) {
  px = -wrap_x;
} 

int my = -diy;
int py = +diy;
int wrap_y = (sy - 1) * diy;
if ( y == 0 ) {
  my = wrap_y;
} else if ( y == sy - 1 ) {
  py = -wrap_y;
} 

int mz = -diz;
int pz = +diz;
int wrap_z = (sz - 1) * diz;
if ( z == 0 ) {
  mz = wrap_z;
} else if ( z == sz - 1 ) {
  pz = -wrap_z;
} 


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

/*
 *  Precalculate derivatives
 */
ftype dExy = Ex[i + py] - Ex[i];
ftype dExz = Ex[i + pz] - Ex[i];

ftype dEyx = Ey[i + px] - Ey[i];
ftype dEyz = Ey[i + pz] - Ey[i];

ftype dEzx = Ez[i + px] - Ez[i];
ftype dEzy = Ez[i + py] - Ez[i];


/*
 *  Precalculate averaged E
 */
ftype aExy = Ex[i + py] + Ex[i];
ftype aExz = Ex[i + pz] + Ex[i];

ftype aEyx = Ey[i + px] + Ey[i];
ftype aEyz = Ey[i + pz] + Ey[i];

ftype aEzx = Ez[i + px] + Ez[i];
ftype aEzy = Ez[i + py] + Ez[i];


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

int pml_xn_thickness = 8;

if ( x < pml_xn_thickness ) {
    const size_t ir = x;                          // index into pml parameters
    const size_t ip = z + y * sz + ir * sz * sy;  // linear index into Psi
    Psi_xn_Hy[ip] = px0hn[ir] * Psi_xn_Hy[ip] + px1hn[ir] * dEzx;
    Psi_xn_Hz[ip] = px0hn[ir] * Psi_xn_Hz[ip] + px1hn[ir] * dEyx;
    pHyi += Psi_xn_Hy[ip];
    pHzi -= Psi_xn_Hz[ip];
}int pml_xp_thickness = 8;

if ( sx > x && x >= sx - pml_xp_thickness ) {
    const size_t ir = (sx - 1) - x;                     // index into pml parameters
    const size_t ip = z + y * sz + ir * sz * sy;  // linear index into Psi
    Psi_xp_Hy[ip] = px0hp[ir] * Psi_xp_Hy[ip] + px1hp[ir] * dEzx;
    Psi_xp_Hz[ip] = px0hp[ir] * Psi_xp_Hz[ip] + px1hp[ir] * dEyx;
    pHyi += Psi_xp_Hy[ip];
    pHzi -= Psi_xp_Hz[ip];
}int pml_yn_thickness = 8;

if ( y < pml_yn_thickness ) {
    const size_t ir = y;                          // index into pml parameters
    const size_t ip = z + x * sz + ir * sz * sx;  // linear index into Psi
    Psi_yn_Hx[ip] = py0hn[ir] * Psi_yn_Hx[ip] + py1hn[ir] * dEzy;
    Psi_yn_Hz[ip] = py0hn[ir] * Psi_yn_Hz[ip] + py1hn[ir] * dExy;
    pHxi -= Psi_yn_Hx[ip];
    pHzi += Psi_yn_Hz[ip];
}int pml_yp_thickness = 8;

if ( sy > y && y >= sy - pml_yp_thickness ) {
    const size_t ir = (sy - 1) - y;                     // index into pml parameters
    const size_t ip = z + x * sz + ir * sz * sx;  // linear index into Psi
    Psi_yp_Hx[ip] = py0hp[ir] * Psi_yp_Hx[ip] + py1hp[ir] * dEzy;
    Psi_yp_Hz[ip] = py0hp[ir] * Psi_yp_Hz[ip] + py1hp[ir] * dExy;
    pHxi -= Psi_yp_Hx[ip];
    pHzi += Psi_yp_Hz[ip];
}int pml_zn_thickness = 8;

if ( z < pml_zn_thickness ) {
    const size_t ir = z;                          // index into pml parameters
    const size_t ip = y + x * sy + ir * sy * sx;  // linear index into Psi
    Psi_zn_Hx[ip] = pz0hn[ir] * Psi_zn_Hx[ip] + pz1hn[ir] * dEyz;
    Psi_zn_Hy[ip] = pz0hn[ir] * Psi_zn_Hy[ip] + pz1hn[ir] * dExz;
    pHxi += Psi_zn_Hx[ip];
    pHyi -= Psi_zn_Hy[ip];
}int pml_zp_thickness = 8;

if ( sz > z && z >= sz - pml_zp_thickness ) {
    const size_t ir = (sz - 1) - z;                     // index into pml parameters
    const size_t ip = y + x * sy + ir * sy * sx;  // linear index into Psi
    Psi_zp_Hx[ip] = pz0hp[ir] * Psi_zp_Hx[ip] + pz1hp[ir] * dEyz;
    Psi_zp_Hy[ip] = pz0hp[ir] * Psi_zp_Hy[ip] + pz1hp[ir] * dExz;
    pHxi += Psi_zp_Hx[ip];
    pHyi -= Psi_zp_Hy[ip];
}

/*
 *  Update H
 */
// Save old H for averaging
ftype Hx_old = Hx[i];
ftype Hy_old = Hy[i];
ftype Hz_old = Hz[i];

// H update equations
Hx[i] -= dt * (dEzy - dEyz + pHxi);
Hy[i] -= dt * (dExz - dEzx + pHyi);
Hz[i] -= dt * (dEyx - dExy + pHzi);

// Average H across timesteps
ftype aHxt = Hx[i] + Hx_old;
ftype aHyt = Hy[i] + Hy_old;
ftype aHzt = Hz[i] + Hz_old;

/*
 *  Calculate unscaled S components at H locations
 */
__global ftype *oSxy = oS + 0 * field_size;
__global ftype *oSyz = oS + 1 * field_size;
__global ftype *oSzx = oS + 2 * field_size;
__global ftype *oSxz = oS + 3 * field_size;
__global ftype *oSyx = oS + 4 * field_size;
__global ftype *oSzy = oS + 5 * field_size;

oSxy[i] = aEyx * aHzt;
oSxz[i] = -aEzx * aHyt;
oSyz[i] = aEzy * aHxt;
oSyx[i] = -aExy * aHzt;
oSzx[i] = aExz * aHyt;
oSzy[i] = -aEyz * aHxt;
==========================================
/*
 *  Update E-field, including any PMLs.
 *
 *  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)
 *
 *  OpenCL args:
 *   E, H, dt, S, oS
 */



typedef float ftype;


/*
 * Field size info
 */
const size_t sx = 395;
const size_t sy = 345;
const size_t sz = 73;
const size_t field_size = sx * sy * sz;

//Since we use i to index into Ex[], Ey[], ... rather than E[], do nothing if
// i is outside the bounds of Ex[].
if (i >= field_size) {
    PYOPENCL_ELWISE_CONTINUE;
}


/*
 * Array indexing
 */
// Given a linear index i and shape (sx, sy, sz), defines x, y, and z
//  as the 3D indices of the current element (i).
// (ie, converts linear index [i] to field indices (x, y, z)
const size_t x = i / (sz * sy);
const size_t y = (i - x * sz * sy) / sz;
const size_t z = (i - y * sz - x * sz * sy);

// Calculate linear index offsets corresponding to offsets in 3D
// (ie, if E[i] <-> E(x, y, z), then E[i + diy] <-> E(x, y + 1, z)
const size_t dix = sz * sy;
const size_t diy = sz;
const size_t diz = 1;


/*
 * Pointer math
 */
//Pointer offsets into the components of a linearized vector-field
// (eg. Hx = H + XX, where H and Hx are pointers)
const size_t XX = 0;
const size_t YY = field_size;
const size_t ZZ = field_size * 2;

//Define pointers to vector components of each field (eg. Hx = H + XX)
__global ftype *Ex = E + XX;
__global ftype *Ey = E + YY;
__global ftype *Ez = E + ZZ;

__global ftype *Hx = H + XX;
__global ftype *Hy = H + YY;
__global ftype *Hz = H + ZZ;


/*
 * Implement periodic boundary conditions
 *
 * mx ([m]inus [x]) gives the index offset of the adjacent cell in the minus-x direction.
 * In the event that we start at x == 0, we actually want to wrap around and grab the cell
 * x_{-1} == (sx - 1) instead, ie. mx = (sx - 1) * dix .
 *
 * px ([p]lus [x]) gives the index offset of the adjacent cell in the plus-x direction.
 * In the event that we start at x == (sx - 1), we actually want to wrap around and grab
 * the cell x_{+1} == 0 instead, ie. px = -(sx - 1) * dix .
 */

int mx = -dix;
int px = +dix;
int wrap_x = (sx - 1) * dix;
if ( x == 0 ) {
  mx = wrap_x;
} else if ( x == sx - 1 ) {
  px = -wrap_x;
} 

int my = -diy;
int py = +diy;
int wrap_y = (sy - 1) * diy;
if ( y == 0 ) {
  my = wrap_y;
} else if ( y == sy - 1 ) {
  py = -wrap_y;
} 

int mz = -diz;
int pz = +diz;
int wrap_z = (sz - 1) * diz;
if ( z == 0 ) {
  mz = wrap_z;
} else if ( z == sz - 1 ) {
  pz = -wrap_z;
} 


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


/*
 * Calculate S from oS (pre-calculated components)
 */
__global ftype *Sx = S + XX;
__global ftype *Sy = S + YY;
__global ftype *Sz = S + ZZ;

// Use unscaled S components from H locations
__global ftype *oSxy = oS + 0 * field_size;
__global ftype *oSyz = oS + 1 * field_size;
__global ftype *oSzx = oS + 2 * field_size;
__global ftype *oSxz = oS + 3 * field_size;
__global ftype *oSyx = oS + 4 * field_size;
__global ftype *oSzy = oS + 5 * field_size;

ftype s_factor = dt * 0.125;
Sx[i] = (oSxy[i] + oSxz[i] + oSxy[i + my] + oSxz[i + mz]) * s_factor;
Sy[i] = (oSyz[i] + oSyx[i] + oSyz[i + mz] + oSyx[i + mx]) * s_factor;
Sz[i] = (oSzx[i] + oSzy[i] + oSzx[i + mx] + oSzy[i + my]) * s_factor;