% Power law fit to White formula on 1 to 10^3
% Could split int otwo regions and do fit to each one.
clear
x=logspace(1,3);
%y=24./x+6./(1+sqrt(x))+0.4; %While
%y=24./x+4./x.^(1/3);        % Putnam
y=24./x.*(1+0.15.*x.^.687);  %Schiller-Naumann
[m,b]=polyfit(log(x),log(y),1);
figure(1)
plot(log(x),log(y),'o',log(x),m(1)*log(x)+m(2))
figure(2)
plot(x,y,'ok','MarkerSize',3)
hold on
plot(x,x.^m(1)*exp(m(2)),'k-')
%plot(x,12./sqrt(x),'r')
%axis([10,1000,0,10])
 xlabel('Re')
 ylabel('C_D')
hold off
display([num2str(exp(m(2))) ' x^' num2str(m(1))])

x=linspace(10,1000,100);
%y=24./x+6./(1+sqrt(x))+0.4; %While
%y=24./x+4./x.^(1/3);        % Putnam
y=24./x.*(1+0.15.*x.^.687);  %Schiller-Naumann
f = @(x,p) p(1).*x.^p(2);
errf = @(p,x,y) sum((y(:)-f(x(:),p)).^2);
p0 = [mean(y) (max(y)-min(y)) (max(x) - min(x))/2];
P = fminsearch(errf,p0,[],x,y); 
x2 = linspace(10,1000,100);
figure(3)
%plot(x,y,'ko',x2,f(x2,P),'b-',x2,12./sqrt(x2),'r')
 plot(x,y,'ok','MarkerSize',3) 
 hold on
 plot(x2,f(x2,P),'k-');
 xlabel('Re')
 ylabel('C_D')
 hold off
display([num2str(P(1)) ' x^' num2str(P(2))])