function runrp
%
%  Name:   runrp.m
%
%  Usage:  runrp.m
%
%  This m-file does a rational polynomial fit to FRF's 
%
%  Version SWD971126
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global DISPLAY MODES DATA iaxes inds

%
% get needed quantities
%

j = sqrt(-1);

%
%  Determine the part of the FRF in the freq band
%

n1 = floor(DISPLAY{iaxes}.cursfreq{1}/DISPLAY{iaxes}.dx) + 1 ;
n2 = floor(DISPLAY{iaxes}.cursfreq{2}/DISPLAY{iaxes}.dx) + 1 ;

%  Set up the variables

switch MODES.ispecfunfit
case 'FRF'
   G =  DATA{1}.G(:,:,:,DISPLAY{iaxes}.iavg{1});
case 'CPS'
   G = DATA{1}.Gxy(:,:,:,DISPLAY{iaxes}.iavg{1});
end

if ~isfield(MODES,'nmodes'),
   MODES.nmodes = 0;
end

modevec = [MODES.nmodes+1:(MODES.nmodes+MODES.nmodfit)];

%
%  Run Rational Polynomial Curve Fit
%

switch MODES.istatmeth,
case 'None',  % NO STATISTICAL ANALYSIS   
   
   [res,freq,damp,kterms] = rpfit(G,n1,n2,2*pi*DISPLAY{iaxes}.xvec.',MODES.nmodfit,MODES.nkterms);

   for im = 1:length(modevec),
   
      MODES.residue{modevec(im)} = res(:,:,im).' / (2 * pi * freq(im))^2 ;  % ASSUMES ACCELERATION FRF
      MODES.Freq{modevec(im)} = freq(im);
      MODES.Damp{modevec(im)} = damp(im);
      MODES.kterms{modevec(im)} = kterms(im);
      MODES.timestamp{modevec(im)} = now;
      MODES.methodstr{modevec(im)} = 'Rational Polynomial';

   end
   
   if MODES.nmodes == 0,
      MODES.respDOF = DATA{1}.respDOF;
      MODES.refDOF = DATA{1}.refDOF;
   end
   
   MODES.nmodes = MODES.nmodes + MODES.nmodfit;
   
   residues2modes

case 'Monte Carlo',
   
   %  Initialize modal variables
   
   for im = 1:length(modevec),
      MODES.Freqall{modevec(im)} = zeros(MODES.nstatiter,1);
      MODES.Dampall{modevec(im)} = zeros(MODES.nstatiter,1);
      MODES.ktermsall{modevec(im)} = zeros(MODES.nstatiter,MODES.nkterms);
      MODES.residueall{modevec(im)} = zeros(MODES.nstatiter,size(DATA{1}.refDOF,1),size(DATA{1}.respDOF,1));
   end
   
   if ~isfield(DATA{1},'mag1sig'),
      [DATA{1}.mag1sig,DATA{1}.phs1sig,DATA{1}.coh1sig,DATA{1}.cps1sig] = genstat(DATA{1}.G(:,:,:,DISPLAY{iaxes}.iavg{inds}),DATA{1}.Gxy(:,:,:,DISPLAY{iaxes}.iavg{inds}),DATA{1}.Coh(:,:,:,DISPLAY{iaxes}.iavg{inds}),DATA{1}.navg);
   end
   
   for k = 1:MODES.nstatiter,

      fprintf(1,'Monte Carlo Run #%d of %d\n',k,MODES.nstatiter)
      
      switch k,
      case 1,
         Gmc = G;
      otherwise
         randn('state',sum(100*clock));
         Gmc_mag = abs(G) + DATA{1}.mag1sig .* randn(size(DATA{1}.mag1sig));
         Gmc_phs = angle(G) + DATA{1}.phs1sig .* randn(size(DATA{1}.phs1sig));
         
         Gmc = Gmc_mag .* exp(j * Gmc_phs);     
      end
      
      [res,freq,damp,kterms] = rpfit(Gmc,n1,n2,2*pi*DISPLAY{iaxes}.xvec.',MODES.nmodfit,MODES.nkterms);

      for im = 1:length(modevec),
   
         MODES.residueall{modevec(im)}(k,:,:) = res(:,:,im).' / (2 * pi * freq(im))^2 ;  % ASSUMES ACCELERATION FRF
         MODES.Freqall{modevec(im)}(k,:) = freq(im);
         MODES.Dampall{modevec(im)}(k,:) = damp(im);
         MODES.ktermsall{modevec(im)}(k,:) = kterms(im);
         
      end
        
   end
   
   for im = 1:length(modevec),
      MODES.timestamp{modevec(im)} = now;
      MODES.methodstr{modevec(im)} = 'Rational Polynomial';      
   end
   
   if MODES.nmodes == 0,
      MODES.respDOF = DATA{1}.respDOF;
      MODES.refDOF = DATA{1}.refDOF;
   end
   
   MODES.nmodes = MODES.nmodes + MODES.nmodfit;
   
   % Compute Mean and 95% CI for Freq, Damp, and residue
   
   for im = 1:length(modevec),
   
      MODES.Freq{modevec(im)} = mean(MODES.Freqall{modevec(im)});
      MODES.Damp{modevec(im)} = mean(MODES.Dampall{modevec(im)});
      
      MODES.Freq95CI{modevec(im)} = 2*std(MODES.Freqall{modevec(im)});
      MODES.Damp95CI{modevec(im)} = 2*std(MODES.Dampall{modevec(im)});
           
      MODES.residue{modevec(im)} =  shiftdim(mean(MODES.residueall{modevec(im)},1),1);    
      MODES.residue95CI{modevec(im)} =  2*shiftdim(std(real(MODES.residueall{modevec(im)}),0,1),1) + ...
         sqrt(-1)*2*shiftdim(std(imag(MODES.residueall{modevec(im)}),0,1),1);
      
   end
   
   % Compute modal vectors, etc.
   
   residues2modes
   
case 'Bootstrap',
   
   switch MODES.ispecfunfit
   case 'FRF'
      G =  DATA{1}.G;
   case 'CPS'
      G = DATA{1}.Gxy;
   end

   %  Define the number of bootstrap ensembles
   
   nboot = size(G,4);
   
   if nboot == 1,
      errordlg('ERROR: This dataset does not contain multiple ensembles - no bootstrap analysis possible')
      return
   end

   %  Initialize modal variables
   
   for im = 1:length(modevec),
      MODES.Freqall{modevec(im)} = zeros(MODES.nstatiter,1);
      MODES.Dampall{modevec(im)} = zeros(MODES.nstatiter,1);
      MODES.ktermsall{modevec(im)} = zeros(MODES.nstatiter,MODES.nkterms);
      MODES.residueall{modevec(im)} = zeros(MODES.nstatiter,size(DATA{1}.refDOF,1),size(DATA{1}.respDOF,1));
   end
   

   %
   %  BOOTSTRAP Analysis
   %
   
   
   %  Start the loop
      
   for k = 1:MODES.nstatiter,
         
      fprintf(1,'Bootstrap Run #%d of %d\n',k,MODES.nstatiter)
         
      %  Select which ensembles to use in this bootstrap run
         
      rand('state',sum(100*clock));
      rrand = floor(rand(nboot,1)*nboot+1);
      
      %  Form the averaged FRF for this run
         
      Gboot = mean(G(:,:,:,rrand),4);
         
      %  Run the rp fit
      
      [res,freq,damp,kterms] = rpfit(Gboot,n1,n2,2*pi*DISPLAY{iaxes}.xvec.',MODES.nmodfit,MODES.nkterms);

      for im = 1:length(modevec),
   
         MODES.residueall{modevec(im)}(k,:,:) = res(:,:,im).' / (2 * pi * freq(im))^2 ;  % ASSUMES ACCELERATION FRF
         MODES.Freqall{modevec(im)}(k,:) = freq(im);
         MODES.Dampall{modevec(im)}(k,:) = damp(im);
         MODES.ktermsall{modevec(im)}(k,:) = kterms(im);
         
      end
      
      %  End the loop
         
   end
   
   for im = 1:length(modevec),
      MODES.timestamp{modevec(im)} = now;
      MODES.methodstr{modevec(im)} = 'Rational Polynomial';      
   end
   
   if MODES.nmodes == 0,
      MODES.respDOF = DATA{1}.respDOF;
      MODES.refDOF = DATA{1}.refDOF;
   end
   
   MODES.nmodes = MODES.nmodes + MODES.nmodfit;
   
   % Compute Mean and 95% CI for Freq, Damp, and residue
   
   for im = 1:length(modevec),
   
      MODES.Freq{modevec(im)} = median(MODES.Freqall{modevec(im)});
      MODES.Damp{modevec(im)} = median(MODES.Dampall{modevec(im)});
      
      MODES.Freq95CI{modevec(im)} = boot95ci(MODES.Freqall{modevec(im)});
      MODES.Damp95CI{modevec(im)} = boot95ci(MODES.Dampall{modevec(im)});      
      
      MODES.residue{modevec(im)} =  shiftdim(median(MODES.residueall{modevec(im)},1),1);
      
      for iref = 1:size(MODES.residue{modevec(im)},1),
         for iresp = 1:size(MODES.residue{modevec(im)},2),
      
            MODES.residue95CI{modevec(im)}(iref,iresp) =  boot95ci(real(MODES.residueall{modevec(im)}(:,iref,iresp))) + ...
               sqrt(-1)*boot95ci(imag(MODES.residueall{modevec(im)}(:,iref,iresp))) ;
            
         end
      end
      
   end
   
   % Compute modal vectors, etc.
   
   residues2modes
   

   
end

return


