#include "matrixMul.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

//works................................

////////////////////////////////////////////////////////////////////////////////
// export C interface
extern "C"

void cpuComputation(double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double);

////////////////////////////////////////////////////////////////////////////////////////////////////////
//
////////////////////////////////////////////////////////////////////////////////////////////////////////
void
cpuComputation(double* Vx, double* Vz, double* c1, double* c2, double* c3, double* c4, double* Txx, double* Tzz, double* Txz, double* source, double dd)
{
    for(int n=0;n<MAX;n++)
	{
		
		for (int j = 2; j<=JE; j++)
		{
			for (int i = 2; i<=IE; i++)
			{	
				Vx[INDEX(i,j,IE)] = Vx[INDEX(i,j,IE)] + c1[INDEX(i,j,IB)] * ( (Txx[INDEX(i,j,IB)]-Txx[INDEX(i-1,j,IB)]) + (Txz[INDEX(i,j,IB)]-Txz[INDEX(i,j-1,IB)]) );
				Vz[INDEX(i,j,IE)] = Vz[INDEX(i,j,IE)] + c1[INDEX(i,j,IB)] * ( (Txz[INDEX(i,j,IB)]-Txz[INDEX(i-1,j,IB)]) + (Tzz[INDEX(i,j,IB)]-Tzz[INDEX(i,j-1,IB)]) );
			}
		}
		for (int j = 2;j<=JE-1;j++)
		{
			for (int i = 2; i<=IE-1; i++)
			{
				Txx[INDEX(i,j,IB)] = Txx[INDEX(i,j,IB)] + c2[INDEX(i,j,IB)] *(Vx[INDEX(i+1,j,IE)]-Vx[INDEX(i,j,IE)]) + c3[INDEX(i,j,IB)] *(Vz[INDEX(i,j+1,IE)]-Vz[INDEX(i,j,IE)]);
				Tzz[INDEX(i,j,IB)] = Tzz[INDEX(i,j,IB)] + c2[INDEX(i,j,IB)] *(Vz[INDEX(i,j+1,IE)]-Vz[INDEX(i,j,IE)]) + c3[INDEX(i,j,IB)] *(Vx[INDEX(i+1,j,IE)]-Vx[INDEX(i,j,IE)]);
				Txz[INDEX(i,j,IB)] = Txz[INDEX(i,j,IB)] + c4[INDEX(i,j,IB)] * ( (Vx[INDEX(i,j+1,IE)]-Vx[INDEX(i,j,IE)]) + (Vz[INDEX(i+1,j,IE)]-Vz[INDEX(i,j,IE)]) );
			}
		}

		//Inject pulse (Both Txx and Tzz are situated in the same place)
		if (n < dd)
		{			
			Tzz[INDEX((IE/2),(JE/2),IB)] =  source[n];			
			Txx[INDEX((IE/2),(JE/2),IB)] =  source[n];
		}		
	}	
}

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