function [phic_resp,phic_ref,Freq,Damp,kterms] = rpfit(G,nlines,n1,n2,wi,nmodes,exterms,respDOF,refDOF)
%
%  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 SWD960920
%
%  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.)

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

nfrf = nresp * nref;

j = sqrt(-1);

npts = n2 - n1 + 1;

%
%  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( (p-1)*nlines + n1 : (p-1)*nlines + n2 ,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,wi,m,n);

%
%  Convert from Polynomial to pole-residue form
%

res = zeros(nfrf,nmodes);
kterms = zeros(nfrf,exterms);

for ir = 1:nfrf,

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

   res(ir,:) = 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));

%
%  Compute modal participation factors from Residues
%

mpf = zeros(nfrf,nmodes);  % Modal participation factors

for im = 1:nmodes,

    mpf(:,im) = res(:,im) * (2 * j * imag(poles(im*2-1)))/ poles(im*2-1)^2 * 9.81; % CONVERSION FROM G's to m/s^2

end

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

%
%  Compute complex mode shapes (both response and reference) from modal participation factors
%	(using singular value decomposition of MPF matrix for each mode)
%
%   ASSUME FRFS ARE FROM ACCELS
%
for p = 1:nmodes,

   mpfmat = zeros(nresp,nref-1);
	
   for iref = 1:nref,

       for iresp = 1:nresp,

         dploc = find(Gindex(:,1) == iresp & Gindex(:,2) == iref);

	 mpfmat(iresp,iref) = mpf(dploc,p);

       end

   end

   [u,s,v] = svd(mpfmat);
  
   scfct = sqrt(s(1,1));

   phic_resp(:,p) = u(:,1)*scfct;  
   phic_ref(:,p) = v(:,1)*scfct;

   end

%
%   Account for relative sensor orientations
%

   for iref = 1:nref,

	phic_ref(iref,:) = phic_ref(iref,:) * sign(refDOF(iref,2));

   end

   for iresp = 1:nresp,

	phic_resp(iresp,:) = phic_resp(iresp,:) * sign(respDOF(iresp,2));

   end

return
