% 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

myInf = 4; % it doesn't seem to like Inf
%q=[0; .011; .021; .048; .079; .126; .185; .264; .358; .472; .622; .733; .893; 2*myInf; 4*myInf; 8*myInf; 16*myInf; 32*myInf; 64*myInf; 128*myInf; 256*myInf; 512*myInf; 1024*myInf; 2048*myInf; 4096*myInf; 8192*myInf; 16384*myInf; 32768*myInf; 65536*myInf];
%f=[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;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]/255;

q=[0; .011; .021; .048; .079; .126; .185; .264; .358; .472; .622; .733; .893; 2*myInf; 4*myInf];
f=[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;1;1]/255;

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

% for squadratic comparametric eq'n, we need inverse (reverse roles x and y):
V=[q sqrt(q)]; % vandermonde matrix WITH NO c IN THE MODEL

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

%coeffs=pinv(V)*f;
coeffs=pinv(V)*f;
a=coeffs(1)
b=coeffs(2)
%%%c=coeffs(3)
c=0; % was forced to zero

q=q(1:13); % get rid of the infinity at the end
f=f(1:13); % get rid of the infinity at the end

f_out = a*q+b*sqrt(q)+c*q.^0;
system("rm bacon.ps") % remove else it appends???
gset term postscript portrait
gset output "bacon.ps"
gset title "SQRTIC FIT TO f(q)"
%gset size square
gset size 1,0.5
%axis([0 .9 0 1])
grid
xlabel("quantigraphic unit, q")
ylabel("response function, f(q)")

plot(q, [f f_out],"-@")

% now find the comparametric equation to which the above is a solution
% modified squadraticlike system f(q) 
%eqfq =f==a*q+b*q^(1/2)+c
%eqgfq=g==a*(k*q)+b*(k*q)^(1/2)+c
%qinv =Solve[eqfq, q]
%kqinv=Solve[eqgfq,q]
%qinv1=qinv[[1]][[1]][[2]]
%qinv2=qinv[[2]][[1]][[2]]
%kqinv1=kqinv[[1]][[1]][[2]]
%kqinv2=kqinv[[2]][[1]][[2]]
%eqgf=qinv1==kqinv1
%gsolution=Solve[eqgf,g]
%g1=gsolution[[1]][[1]][[2]]

% test comparagram: or use the other m file,  figsqrticcomparagram.m
%g=(4*c + 2*b**2*k/a - 4*c*k + 4*f*k + 2*b*sqrt(b**2 - 4*a*c + 4*a*f)*k/a + 2*b*sqrt(k)*sqrt(2*b**2 - 4*a*c + 4*a*f + 2*b*sqrt(b**2 - 4*a*c + 4*a*f) - b**2*k + 4*a*c*k - 4*a*f*k + (b**2 - 4*a*c + 4*a*f)*k)/a)/4;
%k=2;
%fmax=0.967522;
%plot((a*(q/k)+b*sqrt(q/k)+c*(q/k).^0)/fmax,(a*(q)+b*sqrt(q)+c*(q).^0)/fmax);