Index: /issm/trunk-jpl/examples/Inversion/runme.m
===================================================================
--- /issm/trunk-jpl/examples/Inversion/runme.m	(revision 18270)
+++ /issm/trunk-jpl/examples/Inversion/runme.m	(revision 18271)
@@ -1,22 +1,13 @@
-step=3;
+step=1;
 if step==1
 	%Generate observation
 	md = model;
 	md = triangle(md,'DomainOutline.exp',100000);
-    
-	md = setmask(md,'','');
-    
+	md = setmask(md,'all','');
 	md = parameterize(md,'Square.par');
 	md = setflowequation(md,'SSA','all');
 	md.cluster = generic('np',2);
-    
-    md.geometry.surface=md.geometry.surface+100;
-    md.geometry.base = md.geometry.base+100;
-    md.materials.rheology_B(:)=1.8*10^8;
-    md.friction.coefficient(:)=50;
-    md.friction.coefficient(find(md.mesh.x>400000 & md.mesh.x<600000))=10;
-    
 	md = solve(md,StressbalanceSolutionEnum);
-	plotmodel(md,'data',md.friction.coefficient,'caxis',[0 100],'figure',1);
+	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
 	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
 	save model1 md
@@ -25,5 +16,5 @@
 	%Modify rheology, now constant
 	loadmodel('model1.mat');
-	md.friction.coefficient(:)=50;
+	md.materials.rheology_B(:) = 1.8*10^8;
 
 	%results of previous run are taken as observations
@@ -34,5 +25,5 @@
 
 	md = solve(md,StressbalanceSolutionEnum);
-	plotmodel(md,'data',md.friction.coefficient,'caxis',[0 100],'figure',1);
+	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
 	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
 	save model2 md
@@ -45,17 +36,36 @@
 	maxsteps = 20;
 	md.inversion.iscontrol = 1;
-	md.inversion.control_parameters = {'FrictionCoefficient'};
+	md.inversion.control_parameters = {'MaterialsRheologyBbar'};
 	md.inversion.maxsteps = maxsteps;
-	md.inversion.dxmin=10^-6;
-	md.inversion.cost_functions = [103 501];
-	md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,2);
-	md.inversion.cost_functions_coefficients(:,2)=0;
-	md.inversion.min_parameters    = 10^-5*ones(md.mesh.numberofvertices,1);
-	md.inversion.max_parameters    = 100*ones(md.mesh.numberofvertices,1);
+	md.inversion.cost_functions = 101;
+	md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1);
+	md.inversion.min_parameters    = paterson(273)*ones(md.mesh.numberofvertices,1);
+	md.inversion.max_parameters    = paterson(200)*ones(md.mesh.numberofvertices,1);
 
 	%Go solve!
 	md.verbose=verbose(0);
 	md=solve(md,StressbalanceSolutionEnum);
-	plotmodel(md,'data',md.results.StressbalanceSolution.FrictionCoefficient,'caxis',[0 100],'figure',1);
+	plotmodel(md,'data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'figure',1);
 	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
 end
+if step==4
+	%invert for ice rigidity
+	loadmodel('model2.mat');
+
+	%Set up inversion parameters
+	maxsteps = 20;
+	md.inversion.iscontrol = 1;
+	md.inversion.control_parameters = {'MaterialsRheologyBbar'};
+	md.inversion.maxsteps = maxsteps;
+	md.inversion.cost_functions = [101 502];
+	md.inversion.cost_functions_coefficients      = ones(md.mesh.numberofvertices,1);
+	md.inversion.cost_functions_coefficients(:,2) = 10^-16*ones(md.mesh.numberofvertices,1);
+	md.inversion.min_parameters    = paterson(273)*ones(md.mesh.numberofvertices,1);
+	md.inversion.max_parameters    = paterson(200)*ones(md.mesh.numberofvertices,1);
+
+	%Go solve!
+	md.verbose=verbose(0);
+	md=solve(md,StressbalanceSolutionEnum);
+	plotmodel(md,'data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
+end
Index: /issm/trunk-jpl/examples/Inversion/runme2.m
===================================================================
--- /issm/trunk-jpl/examples/Inversion/runme2.m	(revision 18271)
+++ /issm/trunk-jpl/examples/Inversion/runme2.m	(revision 18271)
@@ -0,0 +1,55 @@
+step=1;
+if step==1
+	%Generate observation
+	md = model;
+	md = triangle(md,'DomainOutline.exp',100000);
+	md = setmask(md,'','');
+	md = parameterize(md,'Square.par');
+	md = setflowequation(md,'SSA','all');
+	md.cluster = generic('np',2);
+
+	%START
+	md.geometry.base    = md.geometry.base+100;
+	md.geometry.surface = md.geometry.surface+100;
+	md.friction.coefficient(:)=50;
+	md.friction.coefficient(find(md.mesh.x<600000 & md.mesh.x>400000))=10;
+	%END
+
+	md = solve(md,StressbalanceSolutionEnum);
+
+	save model1 md
+end
+if step==2
+	loadmodel('model1.mat');
+
+	md.friction.coefficient(:)=50;
+
+	%results of previous run are taken as observations
+	md.inversion=m1qn3inversion();
+	md.inversion.vx_obs  = md.results.StressbalanceSolution.Vx;
+	md.inversion.vy_obs  = md.results.StressbalanceSolution.Vy;
+	md.inversion.vel_obs = md.results.StressbalanceSolution.Vel;
+
+	md = solve(md,StressbalanceSolutionEnum);
+	save model2 md
+end
+if step==3
+	%invert for ice rigidity
+	loadmodel('model2.mat');
+
+	%Set up inversion parameters
+	maxsteps = 20;
+	md.inversion.iscontrol = 1;
+	md.inversion.control_parameters = {'FrictionCoefficient'};
+	md.inversion.maxsteps = maxsteps;
+	md.inversion.cost_functions = 101;
+	md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1);
+	md.inversion.min_parameters    = 1*ones(md.mesh.numberofvertices,1);
+	md.inversion.max_parameters    = 100*ones(md.mesh.numberofvertices,1);
+
+	%Go solve!
+	md.verbose=verbose(0);
+	md=solve(md,StressbalanceSolutionEnum);
+	plotmodel(md,'data',md.results.StressbalanceSolution.FrictionCoefficient,'figure',1);
+	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
+end