function run_varkrig(x,y,val,predict,coord_new,predict_new)
%function run_varkrig(x,y,val,predict,coord_new,predict_new)

% Variogram
subplot(2,2,1)
v=variogram([x,y],val,'plotit',true);
% Variogram with fitting of theoretical model
subplot(2,2,2)
[dum,dum,dum,vstruct] = variogramfit(v.distance,v.val,[],[],[],'model','stable','solver','fminsearch');

[xg,yg]=meshgrid(linspace(min(x),max(x),300),linspace(min(y),max(y),300));
%Kriging on regular grid
subplot(2,2,3)
[Zhat,Zvar] = kriging(vstruct,x,y,val,xg,yg);
imagesc(xg(1,:),yg(:,1),Zhat); axis image; axis xy
title('interpolated values')
subplot(2,2,4)
contour(xg(1,:),yg(:,1),Zvar); axis image; axis xy
title('interpolated variance')

if nargin>3
   figure(2)
   [beta,r,s2]=regre2new(val,predict,0.05,[0, 0]');
   Ynew=[ones(size(predict_new,1),1),predict_new]*beta;
   figure(3)
   subplot(1,2,1)
   xfull=[x;coord_new(:,1)]; yfull=[y;coord_new(:,2)]; zfull=[val;Ynew];
   scatter(xfull,yfull,[],zfull,'filled');axis image
   title('scattered values, with linear regression on missing data')
   subplot(1,2,2)
   zinterp=griddata(xfull,yfull,zfull,xg,flipud(yg));
   imagesc(xg(1,:),yg(:,1),zinterp);axis image
   title('scattered values interpolated to grid')
end
