function varargout = getBCparams(varargin)
%GETBCPARAMS  water level, wave height and wave period relations
%
%   Routine to derive, for various locations along the Dutch coast, the
%   wave height based on the water level or the wave period based on the
%   wave height. The relations are obtained from HKV (2005).
%
%   Syntax:
%   varargout = getBCparams(varargin)
%
%   Input:
%   varargin =
%
%   Output:
%   varargout =
%
%   Example
%   getBCparams('ELD', 'Hsig', 'Tp')
%
%   See also

%   --------------------------------------------------------------------
%   Copyright (C) 2009 Delft University of Technology
%       Kees den Heijer
%
%       C.denHeijer@TUDelft.nl	
%
%       Faculty of Civil Engineering and Geosciences
%       P.O. Box 5048
%       2600 GA Delft
%       The Netherlands
%
%   This library is free software; you can redistribute it and/or
%   modify it under the terms of the GNU Lesser General Public
%   License as published by the Free Software Foundation; either
%   version 2.1 of the License, or (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library; if not, write to the Free Software
%   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
%   USA
%   or http://www.gnu.org/licenses/licenses.html, http://www.gnu.org/, http://www.fsf.org/
%   --------------------------------------------------------------------

% Created: 11 Feb 2009
% Created with Matlab version: 7.4.0.287 (R2007a)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords:

%%
OPT = struct(...
    'ncfile', 'HKV2005.nc');

if any(strcmp(fieldnames(OPT), varargin))
    id = find(strcmp(fieldnames(OPT), varargin));
    id = sort([id id+1]);
    OPT = setproperty(OPT, varargin{id});
    varargin(id) = [];
end

%% data based on HKV (2005)
Locations = [cellstr(nc_varget(OPT.ncfile, 'stationName'))...
    cellstr(nc_varget(OPT.ncfile, 'stationAbbreviation'))];
% replace names
Locations{strcmp('Eurogeul', Locations)} = 'Hoek van Holland';
Locations{strcmp('Schiermonnikoog Noord', Locations)} = 'Borkum';

WL = repmat(nc_varget(OPT.ncfile, 'waterlevel')', length(Locations), 1);
WLcell = mat2cell(WL, ones(length(Locations),1), size(WL,2));
Hsig = nc_varget(OPT.ncfile, 'waveheight');
Hsigcell = mat2cell(Hsig, ones(5,1), size(Hsig,2));
Tmmin10 = nc_varget(OPT.ncfile, 'Tm_waveperiod');
Tmmin10cell = mat2cell(Tmmin10, ones(5,1), size(Tmmin10,2));
Tp = nc_varget(OPT.ncfile, 'peak_waveperiod');
Tpcell = mat2cell(Tp, ones(5,1), size(Tp,2));

% remove missing values (nans)
nonnanid = ~isnan(Hsig) & ~isnan(Tmmin10) & ~isnan(Tp);
for i = 1:length(Locations)
    WLcell{i} = WLcell{i}(nonnanid(i,:));
    Hsigcell{i} = Hsigcell{i}(nonnanid(i,:));
    Tmmin10cell{i} = Tmmin10cell{i}(nonnanid(i,:));
    Tpcell{i} = Tpcell{i}(nonnanid(i,:));
end

% create data structure
data = struct(...
    'WL', WLcell,...
    'Hsig', Hsigcell,...
    'Tmmin10', Tmmin10cell,...
    'Tp', Tpcell);

% interpolate Den Helder conditions based on IJmuiden and Eierlandse gat
Locations = [Locations; {'Den Helder' 'DHE'}];
xstp = [79249 98372 106514]';
ystp = [501800 549340 587985]';
xydist = sqrt(diff(xstp).^2 + diff(ystp).^2);

ELDid = strcmp('ELD', Locations(:,2));
IJMid = strcmp('IJM', Locations(:,2));
DHEid = strcmp('DHE', Locations(:,2));

data(DHEid).WL = intersect(data(1).WL, data(3).WL);
for i = 1:length(data(DHEid).WL)
    id1 = data(ELDid).WL == data(DHEid).WL(i);
    id3 = data(3).WL == data(DHEid).WL(i);
    data(DHEid).Hsig(i) = interp1([0 sum(xydist)], [data(IJMid).Hsig(id3); data(ELDid).Hsig(id1)], xydist(1));
    data(DHEid).Tmmin10(i) = interp1([0 sum(xydist)], [data(IJMid).Tmmin10(id3); data(ELDid).Tmmin10(id1)], xydist(1));
    data(DHEid).Tp(i) = interp1([0 sum(xydist)], [data(IJMid).Tp(id3); data(ELDid).Tp(id1)], xydist(1));
end    

%%
vars = {'WL' 'Hsig' 'Tmmin10' 'Tp'};
outvars = {};
for i = 1:length(vars)
    if any(strcmp(vars{i}, varargin))
        outvars(end+1) = vars(i); %#ok<AGROW>
        varargin(strcmp(vars{i}, varargin)) = [];
    end
end

%%
id = zeros(size(varargin));
for i = 1:length(varargin)
    idLocation = strcmpi(Locations(:,1), varargin{i}) |...
        strcmpi(Locations(:,2), varargin{i});
    if any(idLocation)
        id(i) = find(idLocation);
    else
        fprintf('Warning: location "%s" not found.\n', varargin{i})
    end
end

%%
varargout = {};
for j = 1:length(id)
    for i = 1:length(outvars)
        varargout{end+1} = data(id(j)).(outvars{i});
    end
end

if nargout < length(varargout)
    varargout = {varargout};
end