function [A,B] = rpcoeffs(h,wi,m,n)
%
%  Name: rpcoeffs
%
%  Usage: [A,B] = rpcoeffs(h,wi,m,n)
%
%  Description: Computes rational polynomial
%		curve-fit coefficients using Forsythe
%		Orthogonal Polynomials -- 
%		method adopted from Richardson and Formenti
%

%  Version SWD960507
%
%  Author(s): Scott W. Doebling, doebling@lanl.gov
%
%  Institution: Los Alamos National Laboratory,
%		ESA-EA
%
%  Copyright (c) 1996, The Regents of the University
%     	of California. (See top level of code for
%	full copyright statement.)

%  List of variables:
%    m = order of numerator polynomial
%    n = order of denomenator polynomial
%    L = number of frequency bins to curve fit
%    h = FRF matrix values at selected frequencies
%    wi = vector of frequencies to use (indexed 1 to L)
%    Pnum, Pden = matrices of complex orthogonal Forsythe Polynomials (one-sided)
%           row i = ith frequency 
%	    column j = (j-1)th order polynomial
%    qnum = weights for numerator Forsythe poly. (m+1 x 1)
%    qden = weights for denomenator Forsythe poly. (n+1 x 1)
%   

%
%  Normalize frequency variables
%  and make FRF curve 2-sided

j = sqrt(-1);
L1 = length(wi);

wmax = wi(L1);

wi = wi/wmax;
wi = [-wi(L1:-1:1);wi];
xi = j*wi;

h = [conj(h(L1:-1:1,:));h];

L = length(wi);
nh = size(h,2);  % Total number of FRFs to Curve-fit

%
%  Generate Orthogonal Forsythe Polynomials
%

% Numerator Polynomials

qnum = ones(L,1); 	% numerator poly. weights

[Pnum,unum,vnum] = forsythe(xi,qnum,m);


% Denomenator Polynomials

qden = zeros(L,1);

for i = 1:L,
    qden(i) = max(svd(h(i,:)));   % Use CMIF as denomenator poly. weights
end

[Pden,uden,vden] = forsythe(xi,qden,n);


%
%  Form matrices for Least-Squares solution
%

P = Pnum;

SOLROW = zeros(L*nh,(m+1)*nh + n);
WCUM = zeros(L*nh,1);

for hc = 1:nh,

    T = zeros(L,n);

    for i = 1:n,

        T(:,i) = Pden(:,i) .* h(:,hc);

    end   

    W = h(:,hc) .* Pden(:,n+1);

    SOLROW((hc-1)*L + 1 : hc*L , (hc-1)*(m+1) + 1 : hc*(m+1)) = P;
    SOLROW((hc-1)*L + 1 : hc*L , (m+1)*nh + 1 : (m+1)*nh + n) = -T;

    WCUM((hc-1)*L + 1 : hc*L,1 ) = W;

end
 
%
%  Solve for the Forsythe polynomial coefficients
%

SOL = SOLROW \ WCUM;

for hc = 1:nh,

  C(:,hc) = SOL((hc-1)*(m+1) + 1: hc*(m+1));

end

D = [SOL((m+1)*nh + 1:(m+1)*nh + n);1];

%
%  Convert Forsythe Poly. coefficients back to normal poly. coeffs.
%

bkm_num = bkm_comp(m,unum,vnum);
bkm_den = bkm_comp(n,uden,vden);

for hc = 1:nh,

%   Numerator

    for k = 0:m,

 	C1(k+1) = 0;

 	for mm = k:m,

	    C1(k+1) = C1(k+1) + C(mm+1,hc) * bkm_num(k+1,mm+1);

	end

 	C1(k+1) = C1(k+1) / wmax^k;
    end

    A(:,hc) = real(C1(length(C1):-1:1))';

end

%   Denomenator

for k = 0:n,

    D1(k+1) = 0;

    for mm = k:n,

	D1(k+1) = D1(k+1) + D(mm+1) * bkm_den(k+1,mm+1);

    end

    D1(k+1) = D1(k+1) / wmax^k;
	
end

B = real(D1(length(D1):-1:1))';

%%
%%   ALL DONE
%%

return
