% an example of contouring
clear all
loadjg('//usjgofs.whoi.edu/jgofs/eqpac/tt007/ammonium')
whos;

% remap event to yearday
day=eventtoday(event)';
depth=depth';nh4=nh4';

% set up grid
d=[0:10:160]';e=[38:2:70];[eg,dg]=meshgrid(e,d);

% first estimate
average=sum(nh4)/length(nh4);zg=average*ones(size(eg));wg=zeros(size(zg));

% deviation
nh4r=nh4-average;rmserror=sqrt(sum(nh4r.^2)/(length(nh4)-1))

% coarse mapping
radius=10;ehalf=0.0;
[zgn,wg,rmserror,nh4rn]=barnes(e,d,zg,day,depth,nh4,nh4r,radius,ehalf,2.0,1.0);
rmserror

% finer mapping
zg=zgn;nh4r=nh4rn;radius=5;
[zgn,wg,rmserror,nh4rn]=barnes(e,d,zg,day,depth,nh4,nh4r,radius,ehalf,2.0,1.0);
rmserror

% display
figure(1)
cl=contour(e,d,zgn);clabel(cl);
xlabel('Year day');ylabel('Depth');axis('ij');
figure(2)
surf(e,d,zgn);view(20,50);axis('ij');axis([38,70,0,160,0,1.2]);
xlabel('Year day');ylabel('Depth');zlabel('NH4');