function md=control(md)

	%timing
	t1=clock;
	
	%Build all models requested for control simulation
	md.analysis_type='control';
	m=CreateFemModel(md); 

	% figure out number of dof: just for information purposes.
	md.dof=m.nodesets.fsize; %biggest dof number

	%initialize control parameters, gradients and observations
	u_g_obs=m.parameters.u_g_obs;
	p_g=m.parameters.p_g;
	grad_g=zeros(m.nodesets.gsize,1);
		
	%set optimization options.
	options=ControlOptions(m.parameters);

	for n=1:m.parameters.nsteps,
		
		%set options
		options=optimset(options,'MaxFunEvals',m.parameters.maxiter(n));

		disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m.parameters.nsteps)]));

		%initialize inputs, ie m.nparameters on which we invert.
		eval(['inputs.' m.parameters.control_type '=p_g;']);
		inputs.fit=m.parameters.fit(n);

		disp('      computing gradJ...');
		c(n).grad_g=GradJCompute(m,inputs,u_g_obs);
		disp('      done.');

		disp('      normalizing directions...');
		if n>=2,
			c(n).grad_g=Orth(c(n).grad_g,c(n-1).grad_g);
		else
			c(n).grad_g=Orth(c(n).grad_g,{});
		end
		disp('      done.');
		
		%visualize direction.
		if m.parameters.plot
			plot_direction;
		end
		
		disp('      optimizing along gradient direction...'); 
		[search_scalar c(n).J]=fminbnd('objectivefunctionC',-1,1,options,m,p_g,u_g_obs,c(n).grad_g,n);
		disp('      done.');

		disp('      updating parameter using optimized search scalar...');
		p_g=p_g+search_scalar*m.parameters.optscal(n)*c(n).grad_g;
		disp('      done.');

		disp('      constraining the new distribution...');    
		p_g=ControlConstrain(p_g,m.parameters);
		disp('      done.');

		%visualize direction.
		if m.parameters.plot,
			plot_newdistribution;
		end

		%some temporary saving 
		if(mod(n,5)==0),
			solution=controlfinalsol(c,m,p_g);
			save temporary_control_results solution
		end
		disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);

	end

	%Create final solution
	disp('      preparing final velocity solution...');
	solution=controlfinalsol(c,m,p_g);
	disp('      done.');
	
	disp('      load results ontol model...');
	md=loadcontrolfinalsol(md,solution);
	disp('      done.');
	
	
	t2=clock;
	disp(['      overall time spent on control code ' num2str(etime(t2,t1)) ' seconds'])
