% 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

% force it to pass through origin
q=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; .011; .021; .048; .079; .126; .185; .264; .358; .472; .622; .733; .893];
f=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 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 bidratic comparametric eq'n, we fit to quadratic model y=ax^2+bx+c
V2=[q.*q q ones(size(q))]; % vandermonde matrix for quadratic fit
%V3=[q.^3 q.^2 q.^1 q.^0]; % vandermonde matrix for cubic fit
%V4=[q.^4 q.^3 q.^2 q.^1 q.^0]; % vandermonde matrix for quartic fit

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

coeffs2=pinv(V2)*f;
a=coeffs2(1)
b=coeffs2(2)
c=coeffs2(3)

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

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

f2=a*q.^2+b*q.^1+c*q.^0; % quadratic fit to f
% plot(q,f, q,f2) % plots on separate postscript pages
%plot(q,[f f2]) % plots on same page
plot(q,[f f2 f f2],["-";"-";"+";"+"])