function data_prime=griddata_mesh_to_mesh_3d(index,x,y,z,data,x_prime,y_prime,z_prime,varargin);
%GRIDDATA_MESH_TO_MESH_3D - interpolates data from a 3d mesh to a 3d mesh
%
%   griddata_mesh_to_mesh: interpolate data (from mesh with triangulation index,x,y,z)
%   into mesh with triangulation (x_prime,y_prime,z_prime).
%
%   Usage:
%      data_prime=griddata_mesh_to_mesh_3d(index,x,y,z,data,x_prime,y_prime,z_prime,varargin)

	global client_server_mode IMDataCounter uset

	%restrict mesh_prime grids to the ones which are inside the mesh hull 
	in=inhull([x_prime y_prime z_prime],[x y z]);
	pos_hull=find(in);

	%keep only mesh_prime grids inside the hull
	x_prime_hull=x_prime(pos_hull);
	y_prime_hull=y_prime(pos_hull);
	z_prime_hull=z_prime(pos_hull);

	%call griddata_mesh_to_mesh_hulled on this new sub mesh from mesh 2
	if strcmpi(client_server_mode,'yes'),
		r_data_prime_hull=GriddataMeshToMeshHulled_client(index,x,y,data,x_prime_hull,y_prime_hull);
		data_prime_hull=IMdb('select matrix from r_data_prime_hull');
		IMdb('drop r_data_prime_hull');
	else
		data_prime_hull=griddata_mesh_to_mesh_hulled(index,x,y,z,data,x_prime_hull,y_prime_hull,z_prime_hull,varargin);
	end

	%return output
	data_prime=zeros(length(x_prime),1); data_prime(:)=NaN;
	data_prime(pos_hull)=data_prime_hull;
end

function data_prime=griddata_mesh_to_mesh_hulled(index,x,y,z,data,x_prime,y_prime,z_prime,varargin);

	%Analyse data: is it given per node or per element?
	if isempty(varargin{1})
		element_data=(length(data)==size(index,1));%data is given for a given element
		grid_data=(length(data)==length(x));       %data is given for a given grid
		if element_data==0 & grid_data==0, error('griddata_mesh_to_mesh_hulled error message: unknown size of data, specify element or node'), end
	else
		if strcmpi(varargin{1},'node')
			element_data=0;
			grid_data=1;
		elseif strcmpi(varargin{1},'element')
			element_data=1;
			grid_data=0;
		else
			error('griddata_mesh_to_mesh_hulled error message: wrong input argument. Must be ''node'' or ''element''')
		end
	end

	%Some parameters and initialization of data_prime
	nel=length(index);
	data_prime=zeros(length(x_prime),1); data_prime(:)=NaN;

	%compute some quantities that will speed up the process
	x3x1=x(index(:,1))-x(index(:,3));
	y3y1=y(index(:,1))-y(index(:,3));
	x3x2=x(index(:,2))-x(index(:,3));
	y3y2=y(index(:,2))-y(index(:,3));
	x3=x(index(:,3));
	y3=y(index(:,3));

	%triangles jacobian determinant (2*area)
	delta=x(index(:,2)).*y(index(:,3))-y(index(:,2)).*x(index(:,3))-x(index(:,1)).*y(index(:,3))+y(index(:,1)).*x(index(:,3))+x(index(:,1)).*y(index(:,2))-y(index(:,1)).*x(index(:,2));

	%Get area coordinates:
	for n=1:nel,
		if mod(n,10000)==0,
			disp(sprintf('\r%s%.2f%s','      interpolation progress:  ',n/nel*100,' %'));
		end
		%first area coordinate
		area_1=(y3y2(n)*(x_prime-x3(n))-x3x2(n)*(y_prime-y3(n)))/delta(n);
		%second area coordinate
		area_2=(x3x1(n)*(y_prime-y3(n))-y3y1(n)*(x_prime-x3(n)))/delta(n);
		%third area coordinate
		area_3=1-area_1-area_2;

		%get grids in mesh prime that have all area coordinates positive (meaning they are on the same extruded column):
		pos_horiz=find((area_1>=0) & (area_2>=0) & (area_3>=0));
		if ~isempty(pos_horiz),

			%Figure out in which extruded penta it is located
			lower=area_1(pos_horiz)*z(index(n,1))+area_2(pos_horiz)*z(index(n,2))+area_3(pos_horiz)*z(index(n,3));
			upper=area_1(pos_horiz)*z(index(n,4))+area_2(pos_horiz)*z(index(n,5))+area_3(pos_horiz)*z(index(n,6));
			xi=2*(z_prime(pos_horiz)-lower)./(upper-lower)-1;
			pos=find(z_prime(pos_horiz)>=lower & z_prime(pos_horiz)<=upper);
			pos_vert=pos_horiz(pos);
			xi=xi(pos);

			if ~isempty(pos_vert)
				if grid_data
					%mesh 2 grids in pos belong to mesh 1 element n. interpolate value linearly.
					data_prime(pos_vert)=area_1(pos_vert).*(1-xi)/2*data(index(n,1))+area_2(pos_vert).*(1-xi)/2*data(index(n,2))+area_3(pos_vert).*(1-xi)/2*data(index(n,3))...
							+area_1(pos_vert).*(1+xi)/2*data(index(n,4))+area_2(pos_vert).*(1+xi)/2*data(index(n,5))+area_3(pos_vert).*(1+xi)/2*data(index(n,6));
				elseif element_data
					data_prime(pos_vert)=data(n);
				end
			end
		end
	end
	if nel>10000,
		disp(sprintf('\r%s%.2f%s','      interpolation progress:  ',100,' %'));
	end
end
