function v=CARPETinterp3(d,x,y,z)
%
% v=CARPETinterp3(d,x,y,z)
%
% Gives the interpolated values of d in the 3D parallelepiped 
% defined by the vectors x,y,z
% using the 'best' ref. level which contains each point.
%
% REMEBER: For points between the border of a refinement level and a point 
% of the previous (coarser) ref.level the coarser ref.level is used
%
% EXAMPLE:
%  str = CARPETListDatasetHDF5('data/hdf5','phi',1);
%  d   = CARPETLoadDatasetHDF5(str,2);
%  matrixofvalues=CARPETinterp3(d,[1:0.3:5],[-7:0.2:4],[2 2.3 6])
%
% REMEMBER:
% See in the help of matlab the difference between ndgrid end meshgrid.
% Meshgrid switch the first 2 indices of the matrix of coordinates
% and is used for visualization, pay attention making comparisons.
%
% EXAMPLE:
% [xg,yg,zg]=ndgrid([1:0.3:5],[-7:0.2:4],[2 2.3 6]);
% test=xg.*matrixofvalue; 
% It's OK
%
% [xg,yg,zg]=meshgrid([1:0.3:5],[-7:0.2:4],[2 2.3 6]);
% test=xg.*matrixofvalue; 
% It produces an ERROR
%
% Author: Gian Mario Manca, manca@fis.unipr.it
%         http://www.fis.unipr.it/~manca
%
% These scripts are distributed under the GNU General Public License (GPL)
% 11 October 2006

 v=zeros(length(x),length(y),length(z));
 
 for a=[d.rlmin:d.rl]
 %for a=[1:1]
      for b=1:d.c(a)
        
        n=d.num2name(a,b);
        
        xi=find(x>=d.dt(n).Xr(1)&x<=d.dt(n).Xr(end));
        yi=find(y>=d.dt(n).Yr(1)&y<=d.dt(n).Yr(end));
        zi=find(z>=d.dt(n).Zr(1)&z<=d.dt(n).Zr(end));
                
        [xgi,ygi,zgi]=ndgrid(x(xi),y(yi),z(zi));
                
         if length(xi)*length(yi)*length(zi)>=1
           v(xi,yi,zi)=interp3(d.dt(n).Xr,d.dt(n).Yr,d.dt(n).Zr,d.dt(n).val,xgi,ygi,zgi,'linear',0);
           %v(xi,yi,zi)=a; %TEST to see which Ref Level is used
           %v(xi,yi,zi)=b; %TEST to see which Component is used
           %disp(d.names(n));
         end
                   
      end
 end
  
         
end