Index: /issm/trunk/src/m/classes/public/parametercontroloptimization.m
===================================================================
--- /issm/trunk/src/m/classes/public/parametercontroloptimization.m	(revision 3133)
+++ /issm/trunk/src/m/classes/public/parametercontroloptimization.m	(revision 3134)
@@ -10,5 +10,5 @@
 %
 %   Example:
-%      md=parametercontroloptimization(md,'nsteps',30)
+%      md=parametercontroloptimization(md,'nsteps',6)
 
 %process options
@@ -21,61 +21,52 @@
 %variable
 cmax_max=1000;
+cm_maxs=linspace(cmax_max/10,cmax_max,3);
 optscal_max=1000;
+optscals=linspace(optscal_max/10,optscal_max,3);
+fits=[0 2];
 %Kept constant
+md2.verbose=0;
 md2.cm_min=1;
-md2.nsteps=5;
+md2.eps_cm=NaN;
 md2.cm_noisedmp=0;
 md2.cm_maxdmp_slope=0;
 md2.cm_mindmp_slope=0;
-md2.eps_cm=NaN;
+md2.nsteps=getfieldvalue(options,'nsteps',5);
+md2.control_type=getfieldvalue(options,'md2.control_type','drag');
 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 parameters
+best=0;
+count=1;
+for fit=fits,
+	md2.fit=fit*ones(md2.nsteps,1);
+	for cm_max=cm_maxs;
+		md2.cm_max=cm_max;
+		for optscal=optscals;
+			md2.optscal=optscal*ones(md2.nsteps,1);
 
-%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
+			disp(sprintf('\n   Step %i/%i, fit=%i, cm_max=%g, optscal=%g\n',count,length(cm_maxs)*length(optscals)*length(fits),fit,cm_max,optscal))
+			md2=solve(md2,'analysis_type','diagnostic');
+			count=count+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
+			%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(sprintf('   --> current best'))
+				fit_final    =md2.fit;
+				cmmax_final  =md2.cm_max;
+				optscal_final=md2.optscal;
+				best=rel;
 			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(:);
+md.fit    =fit_final;
+md.cm_max =cmmax_final;
+md.optscal=optscal_final;