function femmodel=gradient_core(femmodel,varargin),
%GRADIENT_CORE - Brief compute inverse method gradient direction
% 
%   Usage:
%       femmodel=gradient_core(femmodel,varargin);
% 
%   Examples:
%      femmodel=gradient_core(femmodel);
%      femmodel=gradient_core(femmodel,step,search_scalar);

if nargin==3,
	step=varargin{1};
	orthogonalize=varargin{2};
elseif nargin==1
	step=0;
	orthogonalize=false;
else
	help gradient_core
	error('bad usage');
end

issmprintf(VerboseControl,['   compute cost function gradient']);
[grad norm_list]=Gradj(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);

if orthogonalize,
	issmprintf(VerboseControl,'%s',['   orthogonalization']);
	old_gradient=ControlInputGetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters);
	new_gradient=Orth(grad,old_gradient);
else
	new_gradient=grad;
end

%Check that gradient is clean
norm_grad=norm(new_gradient,inf);
if(norm_grad<=0),     error(['||∂J/∂α||∞ = 0   gradient norm is zero']); end
if(isnan(norm_grad)), error(['||∂J/∂α||∞ = NaN gradient norm is NaN' ]); end

%plug back into inputs:
[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputSetGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,  femmodel.parameters,new_gradient);

%Scale all gradients
[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=ControlInputScaleGradient(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,  femmodel.parameters,norm_list,step);
