function [fval,bv]=zjc_fitMultiSeg3(c,x,y,plotornot)
%	Fit spectrum with classical Brune's model.
% Input:
%	a  - [c1,c2], position of fc0 & fmax0;
% Output:
% 	bv - slops of 3 segment line.
% =============================================================
% JC Zheng. Apr. 9, 2016 @ 325# Zhiyuan Building, CCU Taiwan.
% The slope of 2nd lineseg is restrict as 0.

if ~exist('plotornot','var')
	plotornot=0;
end
	bv=nan(3,1);

	ind1=1:c(1);
	ind2=c(1)+1:c(2);
	ind3=c(2)+1:numel(x);
	
	x1=x(ind1);
	x2=x(ind2);
	x3=x(ind3);
	
	y1=y(ind1);
	y2=y(ind2);
	y3=y(ind3);
	
	[brob1,stats1] = robustfit(x1,y1);
	 brob2(1)=brob1(1)+brob1(2)*x1(end);
	 brob2(2)=0;
	 stats2.resid=y2-brob2(1);
% ----------------------------------------------------------------------
	% No.3 --> ///constrain a fitted curve through specified points 
	x0=x2(end);	%y0=y2(end);
	y0=brob2(1);
	n = 1;							% Degree of polynomial to fit
							% 'C' is the Vandermonde matrix for 'x'
	V(:,n+1) = ones(length(x3),1,class(x3));
	for j = n:-1:1
		V(:,j) = x3.*V(:,j+1);
	end
	C = V;
	
	d = y3;				% 'd' is the vector of target values, 'y'.
	% There are no inequality constraints in this case, i.e., 
	A = [];
	b = [];
	% We use linear equality constraints to force the curve to hit the required point. 
	% In this case, 'Aeq' is the Vandermoonde matrix for 'x0'
	% and 'beq' is the value the curve should take at that point
	Aeq = x0.^(n:-1:0);
	beq = y0;
% 	Options = optimoptions('fmincon');
% 	Options.MaxIter = 100;
	Options = optimset('Display','off','MaxIter',100,'TolX',0.001);
	[p,resnorm,residual] = lsqlin( C, d, A, b, Aeq, beq, [], [], [], Options);
	%yfit = polyval( p, x3 );
% ----------------------------------------------------------------------
	brob3(1)=polyval( p, 0 );
	brob3(2)=p(1);

	bv(1) =  brob1(2);
	bv(2) =  brob2(2);
	bv(3) =  brob3(2);

	if plotornot==1
		yfit = [brob1(1)+brob1(2).*x1;...
				brob2(1)+brob2(2).*x2;...
				brob3(1)+brob3(2).*x3]; 
		figure; hold on;
		plot(x,y);
		plot([x1;x2;x3],yfit,'g','LineWidth',2);
		text(mean(x1),mean(yfit(ind1)),sprintf('%.2f',bv(1)));
		text(mean(x2),mean(yfit(ind2)),sprintf('%.2f',bv(2)));
		text(mean(x3),mean(yfit(ind3)),sprintf('%.2f',bv(3)));
		strtmp=sprintf('fc0=%.2f, fmax0=%.2f',10^x(c(1)),10^x(c(2)));
		title(strtmp);
	end
	
	fval = sum(abs(stats1.resid)) + sum(abs(stats2.resid)) + sum(abs(residual)) ; 
% 	fval = stats1.robust_s + stats2.robust_s + stats3.robust_s;
% 	fval = [brob1(1)+brob1(2)*x1;brob2(1)+brob2(2)*x2;brob3(1)+brob3(2)*x3];
% 	fprintf('%.3e\n',fval);
end