function [param,opt_tws,res_tws,NS] = Optimization_iGrav_TWS_annual_function(tws_grav,ws_series,param_spec,varargin)
% Inputs:
%   tws_grav    ... total water storage as derived from gravimeter
%                   (verctor)
%   ws_series   ... water storage time series used for water balance
%                   equation (matrix, columns = time series)
%   param_spec  ... specifying estimation parameters: structure array
%                   param_spec.ini_val ... initial value (=number of
%                                           columns in (ws_series)
%                   param_spec.min_val ... minimum value of est. parameter  
%                   param_spec.max_val ... maximum value of est. parameter  

%% Set initial parameters
ini_val = param_spec.ini_val;
min_lim = param_spec.min_val;
max_lim = param_spec.max_val;

%% Call optimization toolbox
param = fmincon(@objfn,ini_val,[],[],[],[],min_lim,max_lim);

%% Generate outputs
% total water storage = sum of all water storage time series. To
% subtract time series, use negative 'param_spec.ini_val'
opt_tws = 0;
for j = 1:length(param)
    opt_tws = opt_tws + param(j)*ws_series(:,j);
end
% Residuals: Use the original time series!
res_tws = tws_grav - opt_tws;
NS = mm_efficiency_measure(tws_grav,tws_calc,'KGE2009');

%% Objective function
    function err = objfn(param)
        % total water storage = sum of all water storage time series. To
        % subtract time series, use negative 'param_spec.ini_val'
        tws_calc = 0;
        for i = 1:length(param)
            tws_calc = tws_calc + param(i)*ws_series(:,i);
        end
%         tws_calc = tws_calc - tws_calc(1);
        % Compute resulting error
        err = 1 - mm_efficiency_measure(tws_grav,tws_calc,'KGE2009');
        err = err*10000;
    end %objfn

end % main function