function check_jacobian_error_results(err,tol,varargin)
% evaluates the results of the numerical Jacobian check for each profile
% and catches error if necessary

verbosity = 0;
if numel(varargin)>=1 && ~isempty(varargin{1})
    verbosity = varargin{1};
    assert(isnumeric(verbosity),'verbosity must be numeric, not %s',class(verbosity));
end

% initialize MException
profnames = fieldnames(err); % lisf of all profiles to check
if verbosity
    fprintf(' \n\n **** Jacobian residuals ****** \n\n')
end
for iprof = 1:numel(profnames)
    % for each profile, check the errors per channel for x, xdot and u
    myprof = profnames{iprof};
    varnames = fieldnames(err.(myprof));
    for ivar = 1:numel(varnames)
        myvar = varnames{ivar};
        try
            myerr = err.(myprof).(myvar);
            check_error_per_channel(myerr,tol);
            if verbosity
                print_err_info(myerr,myprof,myvar)
            end
        catch ME
            if ~exist('MEstore','var')
                % MEstore stores all causes of failure
                msg = sprintf('Failed Jacobian test');
                MEstore = MException('JacTest:failed',msg);
            end
            causemsg = sprintf('failed for d%s/d%s',myprof,myvar);
            causeException = MException('JacTest:JacobianFail',causemsg);
            thisME = addCause(ME,causeException); % add cause to this exception
            MEstore = addCause(MEstore,thisME); % append to MEstore structure
        end
    end
end

if exist('MEstore','var') % if something went wrong somewhere, rethrow the exception containing all the info
    throw(MEstore)
end
return

function check_error_per_channel(myErrorForAllChans,tol)
if isstruct(myErrorForAllChans) % error is a structure 
  % for one channel (psi, te, etc) check each error
  channames = fieldnames(myErrorForAllChans);
  for ichan = 1:numel(channames)
    mychan = channames{ichan};
    myerr = myErrorForAllChans.(mychan);
    try
      check_error_value(myerr,tol)
    catch ME
      msg = sprintf('Problem while evaluating derivative component for channel %s',mychan);
      causeException = MException('JacTest:valueCheckFailed',msg);
      ME = addCause(ME,causeException); % add details of cause
      rethrow(ME);
    end
  end
else
   % one value only
   myerr = myErrorForAllChans;
    try
      check_error_value(myerr,tol)
    catch ME
      msg = sprintf('Problem while evaluating derivative component');
      causeException = MException('JacTest:valueCheckFailed',msg);
      ME = addCause(ME,causeException); % add details of cause
      rethrow(ME);
    end
end
  
return

function check_error_value(myerr,tol)
% compare to tolerance or check for other problems.
%throws error if there is a problem.
if isnan(myerr)
    % NaN indicates that Jacobian was not defined (in check_jacobians)
    error('JacTest:JacobianNotDefined','Jacobian not defined while numerical difference exists');
else
    % otherwise check numerical tolerance
    if myerr== Inf
        error('JacTest:numzero','numerical perturbation is zero while analytical perturbation is not zero');
    elseif myerr == -Inf
        error('JacTest:analzero','analytical perturbation is zero while numerical perturbation is not zero');
    elseif myerr>tol
        error('JacTest:exceedTol','Exceeded error tolerance error=%.2e, tol=%.2e',myerr,tol);
    end
end
return

function str = stringOfChannelSettings(channelSettings)
% convert channel settings stucture to string so it can be displayed

channames = fieldnames(channelSettings);
str = '\nSettings per equation:'; % init string
for ichan = 1:numel(channames) % loop over channel names and put the in the string
    mychan = channames{ichan};
    setting = channelSettings.(mychan);
    str = sprintf('%s\n %5s:%7s',str,mychan,setting);
end

return

function str = all_causes(ME)
% get all causes and make a nice string to display

str = '';

if ~isempty(ME.cause)
    for ii=1:numel(ME.cause)
        str = sprintf('%s\n%s',str,ME.message);
        str = sprintf('%s\n%s',str,all_causes(ME.cause{ii})); %recursive call (!)
    end
else
    str = ME.message;
end
return

function print_err_info(myerr,myprof,myvar)

if isstruct(myerr) % error is structure (for every channel)
  components = fieldnames(myerr);
  for icomp = 1:numel(components)
    %%
    mycomp = components{icomp};
    thiserr = myerr.(mycomp);
    derstr = sprintf('d%s/%s(%s)',myprof,myvar,mycomp);
    fprintf('%17s = %-10.3g',derstr,thiserr)
  end
else % single error only
  thiserr = myerr;
  derstr = sprintf('d%s/%s',myprof,myvar);
  fprintf('%14s = %-10.3g\n',derstr,thiserr)
end
fprintf('\n');
%%


return
