function md=steadystate(md);
%STEADYSTATE - compute the velocity and temperature field of a model in steady state.
%
%   Usage:
%      md=steadystate(md)
%

	%timing
	t1=clock;

	numanalyses=7;
	analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticStokesAnalysisEnum;DiagnosticHutterAnalysisEnum;SlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum];
	solution_type=SteadyStateSolutionEnum;

	displaystring(md.verbose,'%s',['create finite element model']);
	femmodel=NewFemModel(md,solution_type,analyses,numanalyses);

	%retrieve parameters
	verbose=femmodel.parameters.Verbose;
	qmu_analysis=femmodel.parameters.QmuAnalysis;
	control_analysis=femmodel.parameters.ControlAnalysis;
	
	%compute solution
	if ~qmu_analysis,
		if ~control_analysis,
			
			displaystring(verbose,'%s',['call computational core']);
			femmodel=steadystate_core(femmodel);

		else,
			
			displaystring(verbose,'%s',['call computational core']);
			femmodel=control_core(femmodel);

		end

		displaystring(verbose,'%s',['write results']);
		md.results.steadystate=OutputResults(femmodel.elements, femmodel.nodes , femmodel.vertices , femmodel.loads , femmodel.materials, femmodel.parameters, femmodel.results);

	else
		%launch dakota driver for diagnostic core solution
		Qmu(femmodel);
	end

	%stop timing
	t2=clock;
	displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);