% Copyright (c) 2010 Kristopher L. Kuhlman (klkuhlm at sandia dot gov)
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%
% this function is part of the set of MATLAB scripts 
% comprising the implementation of a universal kriging
% program that computes estimates of potential via
% cokriging, taking into account boundary conditions.
%
% Kuhlman, K.L., and E. Pardo-Igúzquiza, 2010. Universal
% cokriging of hydraulic heads accounting for boundary 
% conditions, Journal of Hydrology, 384(1–2), 14–25.
% http://dx.doi.org/10.1016/j.jhydrol.2010.01.002

%% generaete data for synthetic example
clear;

plot_figs = 1;

% a 10x10 domain 
xmin = 0.0;         xmax = 10.0;
ymin = 0.0;         ymax = 10.0;

% head data inside the box
xdatamin = 1.0;     xdatamax = 9.0;
ydatamin = 1.0;     ydatamax = 9.0;
npts = 20;

header1 = 'test01';
header2 = 'test02';

% latin hypercube sampling - taken from SCEM-UA
headdata = lhsu([xdatamin,ydatamin],[xdatamax,ydatamax],npts);  

%% heads
headdata(:,3) = sum(headdata(:,1:2),2) + randn(size(headdata,1),1);  % head = x+y

if plot_figs == 1
    figure
    subplot(1,3,1)
    scatter3(headdata(:,1),headdata(:,2),headdata(:,3),'k*');
    daspect([1,1,2]);
    axis([xmin,xmax,ymin,ymax,xmin+ymin,xmax+ymax]);
    title('random sample locations')
    xlabel('X')
    ylabel('Y')
    zlabel('head')
end

save('-ASCII', [header1,'_head.dat'],  'headdata')
save('-ASCII', [header2,'_head.dat'],  'headdata')  % same head for both

%% boundary vectors 
nvec = 5; % this many on each side
r = 0.5;  % move vectors in from edge

% locations of noflow boundary vectors (bottom & top)
nsveclocs = [linspace(xmin+r,xmax-r,nvec)', zeros(nvec,1)+ymin;
           linspace(xmin+r,xmax-r,nvec)', zeros(nvec,1)+ymax];

% vectors along top & bottom
% --------------------------
noflow = [zeros(2*nvec,1), [ones(nvec,1);-ones(nvec,1)]];      

rot = [0,-1;1,0];  % rotation 90 CCW (Strang Sect. 2.6, p116)
chead = noflow*rot;

if plot_figs == 1
    subplot(1,3,2)
    plot(headdata(:,1),headdata(:,2),'k*')
    axis([xmin,xmax,ymin,ymax]);
    hold on;
    quiver(nsveclocs(:,1),nsveclocs(:,2),noflow(:,1),noflow(:,2),0.25,'g');
    quiver(nsveclocs(:,1),nsveclocs(:,2),chead(:,1),chead(:,2),0.25,'r');
    daspect([1,1,1]);
    title('top and bottom no-flow (green) and constant head (red) data')
end

save('-ASCII',[header1,'_nslocs.dat'],  'nsveclocs')
save('-ASCII',[header1,'_nfvecs.dat'],  'noflow')
save('-ASCII',[header1,'_chvecs.dat'],  'chead')

% vectors along left & right sides
% ---------------------------------
ewveclocs = [zeros(nvec,1)+xmin, linspace(ymin+r,ymax-r,nvec)';
             zeros(nvec,1)+ymax, linspace(xmin+r,xmax-r,nvec)'];

noflow = [[ones(nvec,1);-ones(nvec,1)], zeros(2*nvec,1)];      
chead = noflow*rot;

if plot_figs == 1
    subplot(1,3,3)
    plot(headdata(:,1),headdata(:,2),'k*')
    axis([xmin,xmax,ymin,ymax]);
    hold on;
    quiver(ewveclocs(:,1),ewveclocs(:,2),noflow(:,1),noflow(:,2),0.25,'g');
    quiver(ewveclocs(:,1),ewveclocs(:,2),chead(:,1),chead(:,2),0.25,'r');
    daspect([1,1,1]);
    title('left and right no-flow (green) and constant head (red) data')
end

save('-ASCII',[header2,'_ewlocs.dat'],  'ewveclocs')
save('-ASCII',[header2,'_nfvecs.dat'],  'noflow')
save('-ASCII',[header2,'_chvecs.dat'],  'chead')
