step=3;
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);
    
    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.results.StressbalanceSolution.Vel,'figure',2);  
	save model1 md
end
if step==2
	%Modify rheology, now constant
	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);
	plotmodel(md,'data',md.friction.coefficient,'caxis',[0 100],'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;
	md.inversion.control_parameters = {'FrictionCoefficient'};
	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);

	%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.Vel,'figure',2);  
end
