%Which steps to be performed
steps=[1] ;

%Run Steps

% {{{ Mesh Generation #1
if any(steps==1)

	md.miscellaneous.name='PIG.Mesh_generation';

	%Mesh parameters
	domain =['./DomainOutline.exp'];
	hmax=40000;    % maximum element size of the final mesh
	hmin=5000;     % minimum element size of the final mesh
	hinit=10000;   % element size for the initial mesh
	gradation=1.7; % maximum size ratio between two neighboring elements
	err=8;         % maximum error between interpolated and control field

	% Generate an initial uniform mesh (resolution = hinit m)
	md=bamg(model,'domain',domain,'hmax',hinit,'MaxCornerAngle',1);

	%ploting
	plotmodel(md,'data','mesh')

	% Load Velocities
	% http://nsidc.org/data/nsidc-0484.html
	nsidc_vel='../Data/Antarctica_ice_velocity.nc'; 	

	% Get necessary data to build up the velocity grid
	xmin = ncreadatt(nsidc_vel,'/','xmin');
	xmin = strtrim(xmin);  % this is a string, and we need to recover the double value
	xmin = xmin(1:end-2);  % get rid of the unit
	xmin = str2num(xmin);  % convert to double
	
	ymax = ncreadatt(nsidc_vel,'/','ymax'); 
	ymax = strtrim(ymax);  
	ymax = ymax(1:end-2);  
	ymax = str2num(ymax); 
	
	nx = ncreadatt(nsidc_vel,'/','nx');
	ny = ncreadatt(nsidc_vel,'/','ny');
	
	spacing = ncreadatt(nsidc_vel,'/','spacing'); 
	spacing = strtrim(spacing);
	spacing = spacing(1:end-2);  
	spacing = str2num(spacing); 
	
	% Get velocities (Note: You can use ncdisp('file') to see an ncdump)
	vx = double(ncread(nsidc_vel,'vx'));
	vy = double(ncread(nsidc_vel,'vy'));
	
	x=xmin+(0:1:nx)'*spacing; 
	x=double(x);
	y=(ymax-ny*spacing)+(0:1:ny)'*spacing; 
	y=double(y);

	% Interpolate velocities onto coarse mesh
	vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
	vy_obs=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0);
	vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
	clear vx vy x y;

	% Adapt the mesh to minimize error in velocity interpolation
	md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',vel_obs,'err',err);
	
	%ploting
	plotmodel(md,'data','mesh')

	% Convert x,y coordinates (Polar stereo) to lat/lon
	[md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,-1);
	
	% Save model
	save ./Models/PIG.Mesh_generation md;
end
% }}}

% {{{ Masks #2
if any(steps==2) 

	md = loadmodel('./Models/PIG.Mesh_generation');	

	% Load SeaRISe dataset for Antarctica  
	% http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
	searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
	
	%read thickness mask from SeaRISE
	x1=double(ncread(searise,'x1'));
	y1=double(ncread(searise,'y1'));
	thkmask=double(ncread(searise,'thkmask'));
	
	%interpolate onto our mesh vertices
	groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
	groundedice(groundedice<=0)=-1;
	clear thkmask;

	%fill in the md.mask structure
	md.mask.groundedice_levelset=groundedice; %ice is grounded for mask equal one
	md.mask.ice_levelset=-1*ones(md.mesh.numberofvertices,1);%ice is present when negatvie

	%ploting
	plotmodel(md,'data',md.mask.groundedice_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
	
	% Save model
	save ./Models/PIG.SetMask md;
end
% }}}

% {{{ Parameterization #3
if any(steps==3) 

	md = loadmodel('./Models/PIG.SetMask');
	md = parameterize(md,'./Pig.par');

	% Use a MacAyeal flow model
	md = setflowequation(md,'SSA','all');
	
	% Save model
	save ./Models/PIG.Parameterization md;
end
% }}}

% {{{ Control Method #4
if any(steps==4)

	md = loadmodel('./Models/PIG.Parameterization');

	% Control general
	md.inversion.iscontrol=1;
	md.inversion.maxsteps=20;
	md.inversion.maxiter=40;
	md.inversion.dxmin=0.1;
	md.inversion.gttol=1.0e-4;
	md.verbose=verbose('solution',true,'control',true);

	% Cost functions
	md.inversion.cost_functions=[101 103 501];
	md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
	md.inversion.cost_functions_coefficients(:,1)=1;
	md.inversion.cost_functions_coefficients(:,2)=1;
	md.inversion.cost_functions_coefficients(:,3)=8e-15;

	% Controls
	md.inversion.control_parameters={'FrictionCoefficient'};
	md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
	md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);

	% Additional parameters
	md.stressbalance.restol=0.01;
	md.stressbalance.reltol=0.1;
	md.stressbalance.abstol=NaN;

	% Solve
	md.toolkits=toolkits;
	md.cluster=generic('name',oshostname,'np',2);
	md=solve(md,StressbalanceSolutionEnum);

	% Update model friction fields accordingly
	md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;

	plotmodel(md,'data',md.friction.coefficient)

	% Save model
	save ./Models/PIG.Control_drag md;
end
% }}}

% {{{ Plot #5
if any(steps==5)

	md = loadmodel('./Models/PIG.Control_drag');

	plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',...
		'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,...
		'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,...
		'FontSize#all',12,...
		'data',md.initialization.vel,'title','Observed velocity',...
		'data',md.results.StressbalanceSolution.Vel,'title','Modeled Velocity',...
		'data',md.geometry.base,'title','Bed elevation',...
		'data',md.results.StressbalanceSolution.FrictionCoefficient,'title','Friction Coefficient',...
		'colorbar#all','on','colorbartitle#1-2','[m/yr]',...
		'caxis#1-2',([1.5,4000]),...
		'colorbartitle#3','[m]', 'log#1-2',10);
end
% }}}

% {{{ HO #6
if any(steps==6)

	% Load Model

	% Disable inversion

	% Extrude Mesh

	% Set Flowequation

	% Solve

	% Save Model

end
% }}}

% {{{ Plot #7
if any(steps==7)

	mdHO = loadmodel('./Models/PIG.ModelHO');
	mdSSA = loadmodel('./Models/PIG.Control_drag');

	basal=find(mdHO.mesh.vertexonbase);
	surf=find(mdHO.mesh.vertexonsurface);

	plotmodel(mdHO,'nlines',3,'ncols',2,'unit#all','km','axis#all','equal',...
						'xlim#all',[min(mdHO.mesh.x) max(mdHO.mesh.x)]/10^3,...
						'ylim#all',[min(mdHO.mesh.y) max(mdHO.mesh.y)]/10^3,...
						'FontSize#all',12,...
						'data',mdHO.initialization.vel,'title','Observed velocity',...
						'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
						'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
						'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
						'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
						'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
						'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
						'colorbar#all','on','view#all',2,...
						'colorbartitle#all','[m/yr]',...
						'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
end
% }}}
