step=1;
if step==1
	%Generate observation
	md = model;
	md = triangle(md,'DomainOutline.exp',100000);

	%Hands on: STEP 2
	md = setmask(md,'','');

	md = parameterize(md,'Square.par');
	md = setflowequation(md,'SSA','all');
	md.cluster = generic('np',2);

	%Hands on: STEP 1
	md.geometry.base    = md.geometry.base+100;
	md.geometry.surface = md.geometry.surface+100;

	%Hands on: STEP 3
	md.materials.rheology_B(:) = 1.8*10^8;

	%Hands on: STEP 4
	md.friction.coefficient(:)=50;
	md.friction.coefficient(find(md.mesh.x<600000 & md.mesh.x>400000))=10;

	md = solve(md,StressbalanceSolutionEnum);

	plotmodel(md,'data',md.friction.coefficient,'figure',1);
	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);

	save model1 md
end
if step==2
	loadmodel('model1.mat');

	%Hands on: STEP 1
	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);

	plotmodel(md,'data',md.friction.coefficient,'figure',1);
	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);

	save model2 md
end
if step==3
	%invert for ice rigidity
	loadmodel('model2.mat');

	%Set up inversion parameters
	maxsteps = 20;
	md.inversion.iscontrol = 1;

	%Hands on: STEP 1
	md.inversion.control_parameters = {'FrictionCoefficient'};
	md.inversion.maxsteps = maxsteps;

	%Hands on: STEP 2
	md.inversion.cost_functions = 101;
	md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1);

	%Hands on: STEP 3
	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,'caxis',[10 50]);
	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
end
