function  [z, out_data, res_vec, varargout] = Newton_Raphson_solver(z, f_and_jac, xxdotu_z, ...
  model, params, g, v, stap, geop, trap, verbosity)
% -- NEWTON SOLVER --
% starting from initial guess z, iteratively look for solution z
% for which f(z)=0
% the solver is blind to the underlying physics contained in f_and_jac
% each step requires evaluation of f(z) and df_dz(z)
% OUTPUT: the parameter vector z solving f(x(z), xdot(z), u(z)) = 0
% INPUTS: * function handles f_and_jac, x_z, xdot_z, u_z provide the
%         parametrization for the case of interest

% parameters Newton solver
tau_max = 1; red_tau = 0.5; tau_min = 1e-3; nr_iter = 150; res_tol = 1e-8;
% check_jac = true;

res_acc = Inf;
inewt = 1;
res_vec = zeros(1, nr_iter);
while inewt <= nr_iter
  % -----------EVALUATION f AND JAC
  [f, df_dz, stap, geop, trap, inner_jacs] = f_and_jac(z, stap, geop, trap, inewt-1);
%   if check_jac
%     zpert = .01*randn(size(z));
%     ndelta = 100;
%     delta = logspace(-10,1,ndelta);
%     resdiff = zeros(ndelta,1);
%     for i = 1:ndelta
%       fpert = f_and_jac(z+delta(i)*zpert, stap, geop, trap, inewt-1);
%       diffnum = (fpert - f)/delta(i);
%       diffana = df_dz*zpert;
%       resdiff(i) = norm(diffnum - diffana);
%     end
%     loglog(delta, resdiff,'*')
%     pause()
%   end
  % -----------EVALUATION RESIDUE
  res = norm(f);
  res_vec(inewt) = res;
  if res < res_tol
    if verbosity > 0
      fprintf('State equations solved in %u iterations\n', inewt);
    end
    out_data = get_out_data(xxdotu_z, z, g, v, stap, geop, trap, model, params);
    res_vec = res_vec(res_vec ~= 0);
    % varargout containing the jacs of interest for optimization
    % routine
    if nargout == 4
      jacs.df_dz = df_dz;
      % jacs.df_dx = inner_jacs.df_dx;
      % jacs.df_dxdot = inner_jacs.df_dxdot;
      jacs.df_du = inner_jacs.df_du;
      varargout{1} = jacs;
    end
    return
  elseif res < res_acc
    if verbosity > 2
      plot_iteration(xxdotu_z, z, g, v, model);
    end
    tau = tau_max;
    z_acc = z; res_acc = res;
    % -----------NEWTON DESCENT DIRECTION
    [L,U,P] = lu(df_dz);
    z_dir = -U\(L\(P*f));
  else
    res_vec(inewt) = res_acc;
    tau = red_tau*tau;
    z = z_acc;
    if verbosity > 1
      fprintf('it %u: increasing residue; tau reduced to value %.4f\n', inewt, tau);
    end
    if tau < tau_min
      error('RAPTOR:IncreasingRes',...
        'tau (%4.2e) is lower than the minimum value (%4.2e), no solution found (increasing residue)', tau, tau_min);
    end
  end
  z = z + tau*z_dir;
  unphys_flag = z_validity_flag(xxdotu_z, z, g, v, model);
  while unphys_flag ~= -1
    if verbosity > 1
      fprintf('it %u: state not valid; tau reduced to value %.4f\n', inewt, red_tau*tau);
    end
    tau = red_tau*tau;
    z = z_acc + tau*z_dir;
    if tau < tau_min
      error('RAPTOR:UnphysicalState',...
        'tau (%4.2e) is lower than the minimum value (%4.2e), no solution found (unphysical state, flag %.2d)', tau, tau_min, unphys_flag);
    end
    unphys_flag = z_validity_flag(xxdotu_z, z, g, v, model);
  end
  inewt = inewt+1;
  % maximum nr of iterations?
  if inewt == nr_iter + 1
    error('RAPTOR:MaxIter',...
      'maximum number of iterations (%.2d) reached, no solution can be found', nr_iter);
  end
end
end

function flag = z_validity_flag(xxdotu_z, z, g, v, model)
[x, ~, ~, ~, ~, ~] = xxdotu_z(z);
flag = check_state_validity(x, g, v, model);
end

function [] = plot_iteration(xxdotu_z, z, g, v, model)
[x, ~, xdot, ~, ~, ~] = xxdotu_z(z);
RAPTOR_plot_iteration(x, xdot, g, v, model)
pause(1)
end

function out_data = get_out_data(xxdotu_z, z, g, v, stap, geop, trap, model, params)
[x, ~, xdot, ~, u, ~] = xxdotu_z(z);
out_data = RAPTOR_out(x, g, v, xdot, u, 1, stap, geop, trap, model, params);
end