function femmodel=control_core(femmodel)
%CONTROL_CORE - compute the core inversion
%
%   Usage:
%      femmodel=control_core(femmodel);
%

	%recover parameters common to all solutions
	num_controls=femmodel.parameters.InversionNumControlParameters;
	num_responses=femmodel.parameters.InversionNumCostFunctions;
	responses=femmodel.parameters.InversionCostFunctions;
	control_type=femmodel.parameters.InversionControlParameters;
	solution_type=femmodel.parameters.SolutionType;
	nsteps=femmodel.parameters.InversionNsteps;
	maxiter=femmodel.parameters.InversionMaxiterPerStep;
	step_threshold=femmodel.parameters.InversionStepThreshold;
	tol_cm=femmodel.parameters.InversionCostFunctionThreshold;
	gradient_only=femmodel.parameters.InversionGradientOnly;
	dim=femmodel.parameters.MeshDimension;
	isstokes=femmodel.parameters.FlowequationIsstokes;

	%Initialise options with maxiter
	options.MaxIter=femmodel.parameters.InversionMaxiterPerStep;

	%Initialize misfits with a vector of zeros
	J=zeros(nsteps,1);
	search_scalar=1;

	%Get core from solution type
	solutioncore=CorePointerFromSolutionEnum(femmodel.parameters,solution_type);
	adjointcore=AdjointCorePointerFromSolutionEnum(solution_type);

	%Launch once a complete solution to set up all inputs
	issmprintf(VerboseControl,'   preparing initial solution');
	if isstokes,
		eval(['femmodel=' solutioncore '(femmodel);']);
	end

	for n=1:nsteps,

		issmprintf(VerboseControl(),['control method step ' num2str(n) '/' num2str(femmodel.parameters.InversionNsteps)]);
		femmodel.parameters.StepResponses=responses(n,:);

		%In steady state inversion, compute new temperature field now
		if (solution_type==SteadystateSolutionEnum)
			femmodel=steadystate_core(femmodel);
		end

		issmprintf(VerboseControl,'   compute adjoint state');
		eval(['femmodel=' adjointcore '(femmodel);']);

		femmodel=gradient_core(femmodel,n-1,search_scalar==0);

		%Return gradient if asked
		if gradient_only,
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,GradientEnum);
			return;
		end

		issmprintf(VerboseControl,'   optimizing along gradient direction');
		[search_scalar J(n)]=ControlOptimization('objectivefunctionC',0,1,options,femmodel,n,femmodel.parameters);

		issmprintf(VerboseControl,'   updating parameter using optimized search scalar');
		[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar,1);
		issmprintf(VerboseControl,['   value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]);

		%Has convergence been reached?
		converged=controlconvergence(J(n),tol_cm);
		if converged,
			break;
		end

	end

	%generate output
	issmprintf(VerboseControl,'   preparing final velocity solution');
	femmodel.parameters.InversionIscontrol=0;
	eval(['femmodel=' solutioncore '(femmodel);']);

	%Some results not computed by diagnostic or steadystate
	for i=1:num_controls,
		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type(i));
	end
	femmodel.results=AddExternalResult(femmodel.results,JEnum,J);

end %end function
