function retval= OPLSQR_cof(xi,yi,m)
%%ortogonal polinomial least square method
%% calculates coefficients of least square equations
%%Referance : A yirst Course in Numerical Analysis
%% Anthony Ralston,Philip Rabinowitz, Mc Graw Hill ISBN 0-07-051158-6 
%% m degree of polynomial xi yi input data
n=length(xi);
mp2=m+2;
mp1=m+1;
p=zeros(mp2,n);
gamma=zeros(mp1,1);
beta=zeros(mp1,1);
omega=zeros(mp1,1);
alpha=zeros(mp1,1);
b=zeros(mp1,1);
wi=zeros(n,1);
a=zeros(3,mp1);
for i=1:n 
   p(2,i)=1.0;  %%p0
   p(1,i)=0.0;  %%p-1
   wi(i)=1.0;
endfor
gamma(1,1)=0.0;
for i=1:n   
  gamma(1,1)+=wi(i);
endfor
beta(1,1)=0.0;
for j=1:m+1  
  omega(j,1)=0;
  for i=1:n  
    omega(j,1)+=wi(i,1)*yi(i)*p(j+1,i);
  endfor %% i=1:n 
  b(j,1)=omega(j,1)/gamma(j,1);
if  j!=(m+1) 
    alpha(j+1,1)=0.0;
    for i=1:n  
       alpha(j+1,1)+=wi(i)*xi(i)*p(j+1,i)*p(j+1,i)/gamma(j,1);
    endfor %% i=1:n 
    for i=1:n  
      p(j+2,i)=(xi(i)-alpha(j+1,1))*p(j+1,i)-beta(j,1)*p(j,i);
    endfor %% i=1:n 
      gamma(j+1,1)=0;
    for i=1:n  
      gamma(j+1,1)+=wi(i,1)*p(j+2,i)*p(j+2,i);
    endfor %% i=1:n  
    beta(j+1,1)=gamma(j+1,1)/gamma(j,1);
endif
endfor  %% j=1:m+1
for j=1:m+1 
  a(1,j)=b(j,1);
  a(2,j)=alpha(j,1);
  a(3,j)=beta(j,1);
endfor 
retval=a;
endfunction
