function a=powerfit(x,y,p)
%function a=powerfit(x,y,p)
%Least square fit of data using linear combination of powers
%The function will plot results
%Input data arrays x and y of dimensions n*1
%Input "power" array p of dimension m*1, the values could be any real numbers, including zero
%The fit is:
%   yfit=a(1)*x^p(1) + a(2)*x^p(2) + ... +a(m)*x^p(m)
%Note that if any of the x value is zero, the computation will not work since it is a forced
%fit through zero


q=p;
if size(q,2)>1
q=q';
end

m=length(q);
a=zeros(m,1);

for i=1:m
z(:,i)=x.^q(i);
end

a=(z'*z)\(z'*y);

x1=min(x);
x2=max(x);
xout=[x1:(x2-x1)/99:x2]';

z=zeros(length(xout),m);
for i=1:m
z(:,i)=xout.^q(i);
end

yout=z*a;

plot(xout,yout,'-',x,y,'o');