clear

% Exponential fit using fminsearch
% h =  linspace(0,40,1000);
 h = 0:4:40
 D = density2(h*1000);
 
% Define exponential function
% f = @(x,p) p(1) + p(2)*exp(-x./p(3));
f = @(x,p) p(1)*exp(x.*p(2));

% define the error function. 
 errf = @(p,x,y) sum((y(:)-f(x(:),p)).^2);

% an initial guess of the exponential parameters
 p0 = [mean(D) (max(D)-min(D)) (max(h) - min(h))/2];

% search for solution
 P = fminsearch(errf,p0,[],h,D); 

% plot the result
 h2 = linspace(0,40,100);
 figure(1)
 plot(h,D,'ko',h2,f(h2,P),'b-')
  
% Compare to linearization approach
 LD=log(D);
 X=h'; Y=LD';
 B = [ones(length(X),1) X] \ Y;
 Y2=B(2)*h2+B(1);
 figure(2)
 semilogy(h2,exp(Y2),'b') 
 hold
 semilogy(h,D,'ko')
 hold
 
 figure(3)
 plot(h,density2(h*1000),'ko')
 hold
 plot(h2,exp(Y2),'k')
 plot(h2,f(h2,P),'k--')
 hold
 xlabel('h (km)')
 ylabel('\rho (kg/m^3)')
 
 
 
 
 % Output functions
 display('Function:')
 %display([num2str(P(1)) ' + ' num2str(P(2)) ' exp(-x/' num2str(P(3)) ')'])
 % Nonlinear
 display([num2str(P(1)) ' exp(' num2str(P(2)) '*x)'])
 % Linearization
 display([num2str(exp(B(1))) ' exp(' num2str(B(2)) '*x)'])