function [res,Freq,Damp,kterms] = rpfit(G,n1,n2,wi,nmodes,exterms)
%
%  Name: rpfit
%
%  Usage: 
%
%  Description: Performs rational polynomial
%		curve-fit using Forsythe
%		Orthogonal Polynomials -- 
%		method adopted from Richardson and Formenti
%
%  Returns complex modes for both reference and response DOF
%
%  Uses a singular value decomposition to convert from modal participation factors
%   to individual mode shapes

%  Version SWD970821
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This matlab source code was originally     %
% developed as part of "DIAMOND" at          %
% Los Alamos National Laboratory. It may     %
% be copied, modified, and distributed in    %
% any form, provided:                        %
%  a) This notice accompanies the files and  %
%     appears near the top of all source     %
%     code files.                            %
%  b) No payment or commercial services are  %
%     received in exchange for the code.     %
%                                            %
% Original copyright is reserved by the      %
% Regents of the University of California,   %
% in addition to Scott W. Doebling, Phillip  %
% J. Cornwell, Erik G. Straser, and Charles  %
% R. Farrar.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nlines = size(G,1);
nresp = size(G,2);
nref = size(G,3);

nfrf = nresp * nref;

j = sqrt(-1);

npts = n2 - n1 + 1;

wifit = wi(n1:n2);

%
%  Form matrix of FRFs to curve-fit (h)
%

h = zeros(npts,nfrf);

nc = 1;

for k = 1:nref,

    for p = 1:nresp,

       h(:,nc) = G(n1:n2,p,k);
       
       Gindex(nc,:) = [p,k];
       
       nc = nc + 1;

    end
    
end


%
%  Compute rational polynomial coefficients
%

n = nmodes*2;

m = nmodes*2 - 1 + exterms;

[A,B] = rpcoeffs(h,wifit,m,n);

%
%  Convert from Polynomial to pole-residue form
%

res = zeros(max(Gindex(:,1)),max(Gindex(:,2)),nmodes);
kterms = zeros(nfrf,exterms);

for ir = 1:nfrf,

   [ires,temp,ikterms] = residue(A(:,ir),B);

   res(Gindex(ir,1),Gindex(ir,2),:) = ires(1:2:length(ires));
   
   if ~isempty(ikterms),
      kterms(ir,:) = ikterms;
   end

end

poles = roots(B);

%
%  Convert Poles to Frequency and Damping
%

Freq = abs(poles(1:2:nmodes*2))/2/pi;
Damp = -real(poles(1:2:nmodes*2))./ abs(poles(1:2:nmodes*2));

[Freq,ip] = sort(Freq);
Damp = Damp(ip);
res = res(:,:,ip);

return
