step=[1 2 3];
if any(step==1)
	%Generate observation
	md=model;
	md=triangle(md,'DomainOutline.exp',100000);
	md=setmask(md,'all','');
	md=parameterize(md,'Square.par');
	md=setflowequation(md,'macayeal','all');
	md.cluster=generic('name',oshostname,'np',2);
	md=solve(md,DiagnosticSolutionEnum);
	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
	plotmodel(md,'data',md.results.DiagnosticSolution.Vel,'figure',2);  
	save model1 md
end
if any(step==2)
	%Modify rheology, now constant
	loadmodel('model1.mat');
	md.materials.rheology_B(:)=1.8*10^8;

	%results of previous run are taken as observations
	md.inversion.vx_obs=md.results.DiagnosticSolution.Vx;
	md.inversion.vy_obs=md.results.DiagnosticSolution.Vy;
	md.inversion.vel_obs=md.results.DiagnosticSolution.Vel;

	md=solve(md,DiagnosticSolutionEnum);
	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
	plotmodel(md,'data',md.results.DiagnosticSolution.Vel,'figure',2);  
	save model2 md
end
if any(step==3)
	%invert for rheology
	loadmodel('model2.mat');

	%Set up inversion parameters
	nsteps=40;
	md.inversion.iscontrol=1;
	md.inversion.control_parameters={'MaterialsRheologyBbar'};
	md.inversion.nsteps=nsteps;
	md.inversion.cost_functions=101*ones(nsteps,1);
	md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1);
	md.inversion.maxiter_per_step=10*ones(nsteps,1);
	md.inversion.step_threshold=.8*ones(nsteps,1);
	md.inversion.gradient_scaling=10^7*ones(nsteps,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=solve(md,DiagnosticSolutionEnum);
	plotmodel(md,'data',md.results.DiagnosticSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'figure',1);
	plotmodel(md,'data',md.results.DiagnosticSolution.Vel,'figure',2);  
end
