steps=[0];

if ismember(1,steps),
	md=mesh(model,'DomainOutline.exp',2000);
	md=geography(md,'','');
	md=parameterize(md,'Square.par');
	md=setelementstype(md,'macayeal','all');

	%Run 1: given basal drag
	centerx=50000;
	centery=50000;
	radius=50000;
	pos=find(sqrt((md.x-centerx).^2+(md.y-centery).^2)<radius); 
	md.drag_coefficient(pos)=md.drag_coefficient(pos)+50*cos(pi/(2*radius)*(sqrt((md.y(pos)-centery).^2+(md.x(pos)-centerx).^2)));

	md.cluster=eval(oshostname);
	md=solve(md,DiagnosticSolutionEnum);
	md.vx_obs=PatchToVec(md.results.DiagnosticSolution.Vx);
	md.vy_obs=PatchToVec(md.results.DiagnosticSolution.Vy);
	md.dummy=md.drag_coefficient;
end

if ismember(2,steps),
	md.drag_coefficient(:)=30;
	md.control_analysis=1;
	md.nsteps=20;
	md.cm_responses=103*ones(md.nsteps,1);
	md.optscal=30*ones(md.nsteps,1);
	md.maxiter=15*ones(md.nsteps,1);
	md.cm_jump=0.99*ones(md.nsteps,1);
	md.cm_noisedmp=0;
	md.weights=1*ones(md.numberofnodes,1);
	md.eps_cm=NaN;
	md=solve(md,DiagnosticSolutionEnum);
	md.results.DiagnosticSolution_NO=md.results.DiagnosticSolution;
end

if ismember(3,steps),
	md.drag_coefficient(:)=30;
	md.control_analysis=1;
	md.nsteps=20;
	md.cm_responses=103*ones(md.nsteps,1);
	md.optscal=30*ones(md.nsteps,1);
	md.maxiter=15*ones(md.nsteps,1);
	md.cm_jump=0.99*ones(md.nsteps,1);
	md.cm_noisedmp=0.6*10^-9; %0.5 too small
	md.weights=1*ones(md.numberofnodes,1);
	md.eps_cm=NaN;
	md=solve(md,DiagnosticSolutionEnum);
	md.results.DiagnosticSolution_YES1=md.results.DiagnosticSolution;
end

if ismember(4,steps),
	md.drag_coefficient(:)=30;
	md.control_analysis=1;
	md.nsteps=20;
	md.cm_responses=103*ones(md.nsteps,1);
	md.optscal=30*ones(md.nsteps,1);
	md.maxiter=15*ones(md.nsteps,1);
	md.cm_jump=0.99*ones(md.nsteps,1);
	md.cm_noisedmp=5*10^-9;
	md.weights=1*ones(md.numberofnodes,1);
	md.eps_cm=NaN;
	md=solve(md,DiagnosticSolutionEnum);
	md.results.DiagnosticSolution_YES2=md.results.DiagnosticSolution;
end

if ismember(5,steps),
	md.drag_coefficient(:)=30;
	md.control_analysis=1;
	md.nsteps=20;
	md.cm_responses=103*ones(md.nsteps,1);
	md.optscal=30*ones(md.nsteps,1);
	md.maxiter=15*ones(md.nsteps,1);
	md.cm_jump=0.99*ones(md.nsteps,1);
	md.cm_noisedmp=8*10^-8;
	md.weights=1*ones(md.numberofnodes,1);
	md.eps_cm=NaN;
	md=solve(md,DiagnosticSolutionEnum);
	md.results.DiagnosticSolution_YES3=md.results.DiagnosticSolution;
end

save model.mdl md

printflag=true;
plotmodel(md,'data',md.results.DiagnosticSolution_NO.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
if printflag, printmodel('SquareTiknonov0','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
plotmodel(md,'data',md.results.DiagnosticSolution_YES1.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
if printflag, printmodel('SquareTiknonov1','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
plotmodel(md,'data',md.results.DiagnosticSolution_YES2.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
if printflag, printmodel('SquareTiknonov2','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
plotmodel(md,'data',md.results.DiagnosticSolution_YES3.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
if printflag, printmodel('SquareTiknonov3','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
plotmodel(md,'data',md.dummy,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
if printflag, printmodel('SquareRealCoeff','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end

if false,
	plotmodel(md,...
		'data',PatchToVec(md.results.DiagnosticSolution_NO.Vy)-md.vy_obs,...
		'data',PatchToVec(md.results.DiagnosticSolution_YES1.Vy)-md.vy_obs,...
		'data',PatchToVec(md.results.DiagnosticSolution_YES2.Vy)-md.vy_obs,...
		'data',PatchToVec(md.results.DiagnosticSolution_YES3.Vy)-md.vy_obs,...
		'caxis#all',[-0.3 0.3]);
end
if false
	plotmodel(md,...
		'data',md.results.DiagnosticSolution_NO.DragCoefficient,...
		'data',md.results.DiagnosticSolution_YES1.DragCoefficient,...
		'data',md.results.DiagnosticSolution_YES2.DragCoefficient,...
		'data',md.results.DiagnosticSolution_YES3.DragCoefficient,...
		'data',md.dummy,...
		'nlines',2,'ncols',3,'caxis#all',[10 90]);
end

%plotmodel(md,'data',md.results.DiagnosticSolution_NO.Vx,'data',md.results.DiagnosticSolution_YES1.Vx,'data',md.results.DiagnosticSolution_YES2.Vx,'data',md.vx_obs);
%plotmodel(md,'data',PatchToVec(md.results.DiagnosticSolution_NO.Vy)-md.vy_obs,'data',PatchToVec(md.results.DiagnosticSolution_YES1.Vy)-md.vy_obs,'data',PatchToVec(md.results.DiagnosticSolution_YES2.Vy)-md.vy_obs);
