Index: /issm/trunk-jpl/examples/LcurveAnalysis/DomainOutline.exp
===================================================================
--- /issm/trunk-jpl/examples/LcurveAnalysis/DomainOutline.exp	(revision 23752)
+++ /issm/trunk-jpl/examples/LcurveAnalysis/DomainOutline.exp	(revision 23752)
@@ -0,0 +1,10 @@
+## Name:DomainOutline
+## Icon:0
+# Points Count  Value
+5 1.000000
+# X pos Y pos
+0 0
+1000000 0
+1000000 1000000
+0 1000000
+0 0
Index: /issm/trunk-jpl/examples/LcurveAnalysis/Front.exp
===================================================================
--- /issm/trunk-jpl/examples/LcurveAnalysis/Front.exp	(revision 23752)
+++ /issm/trunk-jpl/examples/LcurveAnalysis/Front.exp	(revision 23752)
@@ -0,0 +1,10 @@
+## Name:icefront
+## Icon:0
+# Points Count  Value
+5 1.
+# X pos Y pos
+-1000 900000
+-1000 1100000
+1100000 1100000
+1100000 900000
+-1000 900000
Index: /issm/trunk-jpl/examples/LcurveAnalysis/Square.par
===================================================================
--- /issm/trunk-jpl/examples/LcurveAnalysis/Square.par	(revision 23752)
+++ /issm/trunk-jpl/examples/LcurveAnalysis/Square.par	(revision 23752)
@@ -0,0 +1,24 @@
+%Start defining model parameters here
+
+disp('      creating thickness');
+hmin = 300;
+hmax = 1000;
+ymin = min(md.mesh.y);
+ymax = max(md.mesh.y);
+md.geometry.thickness = hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin);
+md.geometry.base      = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
+md.geometry.surface   = md.geometry.base+md.geometry.thickness;
+
+disp('      creating drag');
+md.friction.coefficient=200*ones(md.mesh.numberofvertices,1);
+md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.;
+md.friction.p = ones(md.mesh.numberofelements,1);
+md.friction.q = ones(md.mesh.numberofelements,1);
+
+disp('      creating flow law paramter');
+md.materials.rheology_B=1.8*10^8*ones(md.mesh.numberofvertices,1);
+md.materials.rheology_B(find(md.mesh.x<md.mesh.y))=1.4*10^8;
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+disp('      creating boundary conditions');
+md=SetIceShelfBC(md,'Front.exp');
Index: /issm/trunk-jpl/examples/LcurveAnalysis/runme.m
===================================================================
--- /issm/trunk-jpl/examples/LcurveAnalysis/runme.m	(revision 23752)
+++ /issm/trunk-jpl/examples/LcurveAnalysis/runme.m	(revision 23752)
@@ -0,0 +1,126 @@
+step=[1];
+if any(step==1)
+	% Generate observations
+	md = model;
+	md = triangle(md,'DomainOutline.exp',100000);
+	md = setmask(md,'all','');
+	md = parameterize(md,'Square.par');
+	md = setflowequation(md,'SSA','all');
+	md.cluster = generic('np',2);
+	md = solve(md,'Stressbalance');
+	plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'title','"True" B',...
+		'data',md.results.StressbalanceSolution.Vel,'title','"observed velocities"')
+	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=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,'Stressbalance');
+	plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'title','B first guess',...
+		'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities')
+	save model2 md
+end
+if any(step==3)
+	% Perform L-curve analysis for ice rigidity inversion
+	loadmodel('model2.mat');
+
+	% Set up inversion parameters
+	maxsteps = 20;
+	md.inversion.iscontrol = 1;
+	md.inversion.control_parameters = {'MaterialsRheologyBbar'};
+	md.inversion.maxsteps = maxsteps;
+	md.inversion.cost_functions = [101 502];
+	md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,length(md.inversion.cost_functions));
+	md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1);
+	md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1);
+   md.verbose = verbose('solution',false,'control',true);
+	
+	% Starting L-curve analysis:
+	%
+	% J = Jo + alpha*R.
+	% J: total cost function to be minimized.
+	% Jo: sum of the objective cost function(s) (ex.: 101, or 101+102, or 101+102+103).
+	% R: regularization term (ex.: 502).
+	% alpha: weight of the regularization term.
+	%
+	% L-curve analysis is a method to find the best value for alpha.
+	% Basicaly, it loops over different values of alpha. A plot can be generated for each
+	% respective value of Jo and R (R versus Jo).
+	%
+	min_alpha   = 1.e-20;
+	max_alpha   = 1.e-11;
+	nstep_alpha = 30;
+	log_step    = (log10(max_alpha)-log10(min_alpha))/nstep_alpha;
+	log_alphas  = [log10(min_alpha):log_step:log10(max_alpha)];
+	alphas      = 10.^log_alphas;
+	J           = zeros(length(alphas),length(md.inversion.cost_functions)+1);
+	% Loop over the alphas
+	for i=1:length(alphas),
+		disp('------------------------------------------------------------');
+		disp(['      alpha iteration: ' int2str(i) '/' int2str(length(alphas)) ', alpha value: ' num2str(alphas(i))]);
+		disp('------------------------------------------------------------');
+		md.inversion.cost_functions_coefficients(:,end) = alphas(i);
+		md = solve(md,'Stressbalance');
+		J(i,:) = md.results.StressbalanceSolution.J(end,:); % J comes in [Jo, alphaR, J]. In this example: [101, alpha*502, 101+alpha*502]
+	end
+
+	% Plot the L-curve (log-log)
+	Jo = zeros(length(alphas),1);
+	for i=1:size(J,2)-2,
+		Jo = Jo + J(:,i); % sum of the cost functions (no regularization term). In this example, only 101
+	end
+	R  = J(:,end-1)./alphas(:); % only the regularization term
+	
+	% Tip:
+	% A rescale in the axes may be useful to visualize the L-curve.
+	%
+	% Remember: J = Jo + alpha*R
+	%
+	% Apply a linear transformation on the original axis (Jo, R): 
+   %
+	% |   1       alpha | | Jo  |   | Jo + alpha*R |   |    J    |
+	% |                 | |     | = |              | = |         |
+	% | 1/alpha     1   | |  R  |   | Jo/alpha + R |   | J/alpha |
+	%
+	% Then, use:
+	% Jo2 = J(:,end);
+	% R2  = J(:,end)./alphas(:);
+	% loglog(Jo2,R2,... 
+	%
+	loglog(Jo,R,'-s','Color',[.3 .8 .4],'MarkerSize',6,'MarkerFaceColor','m','MarkerEdgeColor','k','LineWidth',2)
+	voffset=0.1*R;
+	hoffset=0.1*Jo;
+	text(Jo+hoffset,R+voffset,[repmat('\alpha = ',length(alphas),1) num2str(alphas(:),'%2.0e')],...
+            'FontSize',10,'HorizontalAlignment','left','VerticalAlignment','Middle')
+	xlabel('$\mathrm{log}(\mathcal{J}_0$)','Interpreter','latex')
+	ylabel('$\mathrm{log}(\mathcal{R})$','Interpreter','latex')
+end
+if any(step==4)
+	%invert for ice rigidity
+	loadmodel('model2.mat');
+
+	%Set up inversion parameters
+	maxsteps = 20;
+	md.inversion.iscontrol = 1;
+	md.inversion.control_parameters = {'MaterialsRheologyBbar'};
+	md.inversion.maxsteps = maxsteps;
+	md.inversion.cost_functions = [101 502];
+	md.inversion.cost_functions_coefficients      = ones(md.mesh.numberofvertices,1);
+	md.inversion.cost_functions_coefficients(:,2) = 4.e-17*ones(md.mesh.numberofvertices,1); % here you can use the best value found for alpha
+	md.inversion.min_parameters    = cuffey(273)*ones(md.mesh.numberofvertices,1);
+	md.inversion.max_parameters    = cuffey(200)*ones(md.mesh.numberofvertices,1);
+
+	%Go solve!
+	md.verbose=verbose(0);
+	md=solve(md,'Stressbalance');
+	plotmodel(md,'axis#all','tight','data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'title','inferred B',...
+		'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities')
+end