/*!\file:  DataInterpx.cpp
 * \brief  "c" core code for interpolating values onto a mesh.
 */ 

#include "./DataInterpx.h"
#include "../shared/shared.h"

#undef __FUNCT__ 
#define __FUNCT__ "DataInterpx"

int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);

int DataInterpx( Vec* pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods) {

	int i;
	int m,n;
	int data_mesh_m; //local number of rows for vector data_mesh.

	/*output: */
	Vec data_mesh=NULL;
	
	double x_grid,y_grid;
	double xi,eta;
	double G1,G2,G3,G4,data_value;


	/*Some checks on arguments: */
	if ((M<=0)|| (N<=0) || (nods<=0)){
		throw ErrorException(__FUNCT__," error message: nothing to be done according to the dimensions of input matrices and vectors.");
	}
	if(M!=(y_rows-1)){
		throw ErrorException(__FUNCT__," error message: y vector length should be 1 more than data number of rows.");
	}
	if(N!=(x_rows-1)){
		throw ErrorException(__FUNCT__," error message: x vector length should be 1 more than data number of columns.");
	}

	/*Allocate output vector: */
	data_mesh=NewVec(nods);

	/*Nearest interpolation method for now: */
	for ( i=MPI_Lowerrow(nods); i<MPI_Upperrow(nods); i++) {

		x_grid=*(x_mesh+i);
		y_grid=*(y_mesh+i);

		/*Find indices m and n into y and x, for which  y(m)<=y_grids<=y(m+1) and x(n)<=x_grid<=x(n+1)*/
		if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
			
			/*Do a bilinear interpolation. First find parameter coordinates: */
			xi=-1.0+2*(x_grid-x[n])/(x[n+1]-x[n]);
			eta=-1.0+2*(y_grid-y[m])/(y[m+1]-y[m]);
			/*Then find values of data at each grid: (xi,eta)=(0,0),(1,0),(1,1),(0,1) */
			G1=*(data+m*N+n);
			G2=*(data+m*N+n+1);
			G3=*(data+(m+1)*N+n+1);
			G4=*(data+(m+1)*N+n);
			/*Find data_value, by bilinear interpolation: */
			data_value=G1*(1+xi)*(1+eta)/4.0+ G2*(1-xi)*(1+eta)/4.0+G3*(1-xi)*(1-eta)/4.0+G4*(1+xi)*(1-eta)/4.0;

		}
		else{
			data_value=-9999;
		}
		VecSetValues(data_mesh,1,&i,&data_value,INSERT_VALUES);
	}

	/*Assign output pointers:*/
	*pdata_mesh=data_mesh;
}

int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid){

	int foundx=0;
	int foundy=0;
	int i;
	int m=-1;
	int n=-1;

	for (i=0;i<x_rows-1;i++){
		if ( (*(x+i)<=xgrid) && (xgrid<*(x+i+1)) ){
			m=i;
			foundx= 1;
			break;
		}
	}
	if(*(x+x_rows-1)==xgrid){
		m=x_rows-2;
		foundx=1;
	}
	
	for (i=0;i<y_rows-1;i++){
		if ( (*(y+i)<=ygrid) && (ygrid<*(y+i+1)) ){
			n=i;
			foundy= 1;
			break;
		}
	}
	if(*(y+y_rows-1)==ygrid){
		m=y_rows-2;
		foundy=1;
	}

	/*Assign output pointers:*/
	*pm=m;
	*pn=n;
	return foundx*foundy;
}
