Index: /issm/trunk/src/m/classes/public/parametercontroloptimization.m
===================================================================
--- /issm/trunk/src/m/classes/public/parametercontroloptimization.m	(revision 3129)
+++ /issm/trunk/src/m/classes/public/parametercontroloptimization.m	(revision 3129)
@@ -0,0 +1,81 @@
+function md=parametercontroloptimization(md,varargin),
+%PARAMETERCONTROLOPTIMIZATION - parameterization for control method on drag
+%
+%   It is possible to specify the number of steps, values for the
+%   minimum and maximum values of the drag, specify a noise dampening, the 
+%   kind of fit to use or the the optscal.
+%   
+%   Usage:
+%       md=parametercontroloptimization(md,varargin)
+%
+%   Example:
+%      md=parametercontroloptimization(md,'nsteps',30)
+
+%process options
+options=pairoptions(varargin{:});
+
+%Copy model
+md2=md;
+
+%Hard coded parameters
+%variable
+cmax_max=1000;
+optscal_max=1000;
+%Kept constant
+md2.cm_min=1;
+md2.nsteps=5;
+md2.cm_noisedmp=0;
+md2.cm_maxdmp_slope=0;
+md2.cm_mindmp_slope=0;
+md2.eps_cm=NaN;
+md2.maxiter=10*ones(md2.nsteps,1);
+md2.cm_jump=0.99*ones(md2.nsteps,1);
+md2.verbose=0;
+md2.control_analysis=1;
+md2.control_type=getfieldvalue(options,'md2.control_type','drag');
+
+%Get nsteps
+nsteps=getfieldvalue(options,'nsteps',50);
+if mod(nsteps,5)~=0,
+	error('''nsteps'' option should be a multiple of 5')
+end
+
+%loop over the set of 5 steps
+fit_final=[];
+cmmax_final=[];
+optscal_final=[];
+for i=1:nsteps/5,
+	disp(['Step ' num2str(i) '/' num2str(floor(nsteps/5))]);
+	best=0;
+	for fit=[0 2],
+		md2.fit=fit*ones(md2.nsteps,1);
+		for cm_max=linspace(cmax_max/10,cmax_max,3),
+			md2.cm_max=cm_max;
+			for optscal=linspace(optscal_max/10,optscal_max,3),
+				md2.optscal=optscal*ones(md2.nsteps,1);
+
+				%current run
+				md2=solve(md2,'analysis_type','diagnostic');
+
+				%Check misfit difference
+				rel=(md2.results.diagnostic.J(1)-md2.results.diagnostic.J(end))/md2.results.diagnostic.J(1);
+				disp(['   ΔJ/J=' num2str(rel*100) '% for fit=' num2str(fit) ',cm_max=' num2str(cm_max) ', and optscal=' num2str(optscal) ]);
+				if (rel>best),
+					disp(['   --> best for this step']);
+					fit_final(:,i)    =md2.fit;
+					cmmax_final(:,i)  =md2.cm_max;
+					optscal_final(:,i)=md2.optscal;
+					drag=md2.results.diagnostic.parameter;
+					best=rel;
+				end
+			end
+		end
+	end
+	%Plug best drag;
+	md2.drag=drag;
+end
+
+%load final parameters onto initial model
+md.fit    =fit_final(:);
+md.cm_max =cmmax_final(:);
+md.optscal=optscal_final(:);