## orthogonal  polynomial least square 
## A First Course in Numerical Analysis, 
## Anthony Ralston,Philip Rabinowitz, Mc Graw Hill ISBN 0-07-051158-6” 
## Author: M. Turhan ÇOBAN
## Created: 2021-08-22

function retval = orthogonal_least(xi,fi,m)
    #ortogonal polinomial least square method
    #Referance : A First Course in Numerical Analysis
    #Anthony Ralston,Philip Rabinowitz, Mc Graw Hill ISBN 0-07-051158-6 
    #m degree of polynomial xi fi input data
    # int i,j,k;
    n=length(xi);
    mp2=m+2;
    mp1=m+1;
    p=zeros(mp2,n);
    gamma=zeros(mp1);
    beta=zeros(mp1);
    omega=zeros(mp1);
    alpha=zeros(mp1);
    ym=zeros(n);
    b=zeros(mp1);
    wi=ones(n);
    a=zeros(3,mp1);
    for i=1:n
        p(2,i)=1.0;
        p(1,i)=0.0;
        wi(i)=1.0;
    end #for
    gamma(1)=0.0;
    for i=1:n 
        gamma(1)=gamma(1)+wi(i);
    end #for
    beta(1)=0.0;
    for j=1:mp1
        omega(j)=0;
        for i=1:n 
            omega(j)=omega(j)+wi(i)*fi(i)*p(j+1,i);
        end # for i=1:n
        b(j)=omega(j)/gamma(j);
        if j!= (m+1)
            alpha(j+1)=0.0;
            for i=1:n 
                alpha(j+1)+=wi(i)*xi(i)*p(j+1,i)*p(j+1,i)/gamma(j);
            end # for i=1:n
            for i=1:n 
                p(j+2,i)=(xi(i)-alpha(j+1))*p(j+1,i)-beta(j)*p(j,i);
            end # for i=1:n
            gamma(j+1)=0.0;
            for i=1:n 
                gamma(j+1)+=wi(i)*p(j+2,i)*p(j+2,i);
            end # for i=1:n
            beta(j+1)=gamma(j+1)/gamma(j);
        end #if j != m
     end #for j=1:(m+1) 
    for j=1:(m+1) 
       a(1,j)=b(j);
       a(2,j)=alpha(j);
       a(3,j)=beta(j);
    end #j=1:(m+1) 
    retval =a;
endfunction
