//
// This file is part of the Bones source-to-source compiler examples. The C-code
// is largely identical in terms of functionality and variable naming to the code
// found in PolyBench/C version 3.2. For more information on PolyBench/C or Bones
// please use the contact information below.
//
// == More information on PolyBench/C
// Contact............Louis-Noel Pouchet <pouchet@cse.ohio-state.edu>
// Web address........http://polybench.sourceforge.net/
// 
// == More information on Bones
// Contact............Cedric Nugteren <c.nugteren@tue.nl>
// Web address........http://parse.ele.tue.nl/bones/
//
// == File information
// Filename...........benchmark/fdtd-2d-apml.c
// Author.............Cedric Nugteren
// Last modified on...30-May-2012
//

#include "common.h"

// This is 'fdtd-2d-apml', a 2D finite different time domain kernel using an anisotropic perfectly matched layer
int main(void) {
	int i,j,k;
	int ix,iy,iz;
	float clf, tmp;
	float czm_iz, czp_iz;
	float cymh_iy, cyph_iy;
	
	// Declare arrays on the stack
	float czm[CZ+1];
	float czp[CZ+1];
	float cxmh[CXM+1];
	float cxph[CXM+1];
	float cymh[CYM+1];
	float cyph[CYM+1];
	float Ry[CZ+1][CYM+1];
	float Ax[CZ+1][CYM+1];
	float Ex[CZ+1][CYM+1][CXM+1];
	float Ey[CZ+1][CYM+1][CXM+1];
	float Hz[CZ+1][CYM+1][CXM+1];
	float Bza[CZ+1][CYM+1][CXM+1];
	
	// Set the constants
	float mui = 2341;
	float ch = 42;
	
	// Set the input data
	for (i=0; i<=CZ; i++) {
		czm[i] = ((float) i+1) / CXM;
		czp[i] = ((float) i+2) / CXM;
	}
	for (i=0; i<=CXM; i++) {
		cxmh[i] = ((float) i+3) / CXM;
		cxph[i] = ((float) i+4) / CXM;
	}
	for (i=0; i<=CYM; i++) {
		cymh[i] = ((float) i+5) / CXM;
		cyph[i] = ((float) i+6) / CXM;
	}
	for (i=0; i<=CZ; i++) {
		for (j=0; j<=CYM; j++) {
			Ry[i][j] = ((float) i*(j+1) + 10) / CYM;
			Ax[i][j] = ((float) i*(j+2) + 11) / CYM;
			for (k=0; k<=CXM; k++) {
				Ex[i][j][k] = ((float) i*(j+3) + k+1) / CXM;
				Ey[i][j][k] = ((float) i*(j+4) + k+2) / CYM;
				Hz[i][j][k] = ((float) i*(j+5) + k+3) / CZ;
			}
		}
	}
	
	// Perform the computation
	#pragma scop
	for (iz=0; iz<CZ; iz++) {
		for (iy=0; iy<CYM; iy++) {
			czm_iz = czm[iz];
			czp_iz = czp[iz];
			cymh_iy = cymh[iy];
			cyph_iy = cyph[iy];
			
			#pragma species kernel iz:iz,iy:iy+1,0:CXM-1|chunk(0:0,0:1,0:0) ^ iz:iz,iy:iy,0:CXM-1|neighbourhood(0:0,0:0,0:1) ^ iz:iz,iy:iy,0:CXM-1|element  ^ iz:iz,iy:iy,0:CXM-1|element ^ 0:CXM-1|element ^ 0:CXM-1|element -> iz:iz,iy:iy,0:CXM-1|element ^ iz:iz,iy:iy,0:CXM-1|element
			for (ix=0; ix<CXM; ix++) {
				clf = Ex[iz][iy][ix] - Ex[iz][iy+1][ix] + Ey[iz][iy][ix+1] - Ey[iz][iy][ix];
				tmp = (cymh_iy/cyph_iy) * Bza[iz][iy][ix] - (ch/cyph_iy) * clf;
				Hz[iz][iy][ix] = Hz[iz][iy][ix] * (cxmh[ix]/cxph[ix]) + (mui * czp_iz/cxph[ix]) * tmp - (mui * czm_iz/cxph[ix]) * Bza[iz][iy][ix];
				Bza[iz][iy][ix] = tmp;
			}
			#pragma species endkernel fdtd-2d-apml
			
			clf = Ex[iz][iy][CXM] - Ex[iz][iy+1][CXM] + Ry[iz][iy] - Ey[iz][iy][CXM];
			tmp = (cymh_iy/cyph_iy) * Bza[iz][iy][CXM] - (ch/cyph_iy) * clf;
			Hz[iz][iy][CXM] = (cxmh[CXM]/cxph[CXM]) * Hz[iz][iy][CXM] + (mui * czp_iz/cxph[CXM]) * tmp - (mui * czm_iz/cxph[CXM]) * Bza[iz][iy][CXM];
			Bza[iz][iy][CXM] = tmp;
			
			for (ix=0; ix<CXM; ix++) {
				clf = Ex[iz][CYM][ix] - Ax[iz][ix] + Ey[iz][CYM][ix+1] - Ey[iz][CYM][ix];
				tmp = (cymh[CYM]/cyph_iy) * Bza[iz][iy][ix] - (ch/cyph_iy) * clf;
				Hz[iz][CYM][ix] = (cxmh[ix]/cxph[ix]) * Hz[iz][CYM][ix] + (mui * czp_iz/cxph[ix]) * tmp - (mui * czm_iz/cxph[ix]) * Bza[iz][CYM][ix];
				Bza[iz][CYM][ix] = tmp;
			}
			
			clf = Ex[iz][CYM][CXM] - Ax[iz][CXM] + Ry[iz][CYM] - Ey[iz][CYM][CXM];
			tmp = (cymh[CYM]/cyph[CYM]) * Bza[iz][CYM][CXM] - (ch/cyph[CYM]) * clf;
			Hz[iz][CYM][CXM] = (cxmh[CXM]/cxph[CXM]) * Hz[iz][CYM][CXM] + (mui * czp_iz/cxph[CXM]) * tmp - (mui * czm_iz/cxph[CXM]) * Bza[iz][CYM][CXM];
			Bza[iz][CYM][CXM] = tmp;
		}
	}
	#pragma endscop
	
	// Clean-up and exit the function
	fflush(stdout);
	Hz[8][9][3] = Hz[8][9][3];
	return 0;
}