% create a portion of a Vandermonde matrix for least squares approximation
% (first n columns of a vander matrix where n-1 is the order of polynomial

%     q         f(q)*255
%    0.011000   23.757634
%    0.021000   44.856354
%    0.048000   70.447514
%    0.079000   90.480663
%    0.126000  116.176796
%    0.185000  140.685083
%    0.264000  171.232044
%    0.358000  193.215470
%    0.472000  216.505525
%    0.622000  233.182320
%    0.733000  239.433702
%    0.893000  246.718232

q=[.011; .021; .048; .079; .126; .185; .264; .358; .472; .622; .733; .893];
f=[23.757634;44.856354;70.447514;90.480663;116.176796;140.685083;171.232044;193.215470;216.505525;233.182320;239.433702;246.718232]/255;

%V=[q.*q q ones(size(q))]; % vandermonde matrix

% for squadratic comparametric eq'n, we need inverse (reverse roles x and y):
V=[f.*f f ones(size(f))]; % vandermonde matrix

%compute the least squares approximation
%coeffs = (inv(V'*V)*V')*f;

%coeffs=pinv(V)*f;
coeffs=pinv(V)*q;
a=coeffs(1)
b=coeffs(2)
c=coeffs(3)

q_out=a*f.^2+b*f.^1+c*f.^0;
%%%plot(a*f.^2+b*f.^1+c*f.^0, f,   q, f)
%%%plot(q_out, f,   q, f)

system("rm bacon.ps") % remove else it appends???
gset term postscript portrait
gset output "bacon.ps"
gset title "INV QUADRATIC FIT TO f(q)"
%gset size square
gset size 1,0.5

grid
xlabel("quantigraphic unit, q")
ylabel("response function, f(q)")

%plot([q q_out],f,["-";"-";"+";"+"])
plot([q q_out],f)