%
%  Name:   runRP.m
%
%  Usage:  runRP.m
%
%  This m-file does a rational polynomial fit to FRF's 
%
%  Version SWD970710
%
%  Author(s): Phillip Cornwell, cornwell@lanl.gov
%	modified by Scott W. Doebling, doebling@lanl.gov
%
%  Institution: Los Alamos National Laboratory,
%               ESA-EA
%
%  Copyright (c) 1996-1997, The Regents of the University
%       of California. (See top level of code for
%       full copyright statement.)


if exist('Gfit') ~= 1, 

    disp('ERROR: MUST DEFINE THE FREQUENCY BINS TO FIT!')

else
%
% clear variables
%
clear Freqnew Dampnew icountnew hmat

%
% get needed quantities
%
j = sqrt(-1);
flow=str2num(get(onHndlVector(8),'string'));
fhigh=str2num(get(onHndlVector(11),'string'));
nmodfit=str2num(get(onHndlVector(15),'string'));
exterms=str2num(get(onHndlVector(17),'string'));
dppair = str2num(get(onHndlVector(24),'string'));

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

n1=floor((flow-fmin)/df)+1;
n2=floor((fhigh-fmin)/df)+1;


%
%  Set variables
%

wifit = freqfit'*2*pi;

%
%  Run Rational Polynomial Curve Fit
%

disp('Calculating Rational Polynomial Curve Fit. Please be patient...')


if istat == 1,
   
   %  Do statistical analysis
   
   %    Initialize storage variables
      
   Freqallnew = zeros(nmodfit,nstatiter);
   Dampallnew = zeros(nmodfit,nstatiter);
   Phicrespallnew = zeros(size(respDOF,1),nstatiter,nmodfit);
   Phicrefallnew = zeros(size(refDOF,1),nstatiter,nmodfit);
   
   if istatmeth == 1,
   
      %  Monte Carlo Analysis here

   
      if exist('navg1') ~=1 , 
         navg1 = input('Number of Averages?'); 
      end

      [mag1sig,phs1sig] = genstat(G1,Gxy1,Coh1,navg1);
   
      mag1sigfit = zeros(nfitbin*nresp,nref);
      phs1sigfit = zeros(nfitbin*nresp,nref);
   
      for i = 1:nresp,
   
         mag1sigfit((i-1)*nfitbin+1:i*nfitbin,:) = mag1sig((i-1)*nlines+fitpts,:);
         phs1sigfit((i-1)*nfitbin+1:i*nfitbin,:) = phs1sig((i-1)*nlines+fitpts,:);
   
      end
   
      for k = 1:nstatiter,

         fprintf(1,'Monte Carlo Run #%d of %d\n',k,nstatiter)
      
         if k==1,

            Gmc = Gfit;

         else

            randn('state',sum(100*clock));
            Gmc_mag = abs(Gfit) + mag1sigfit .* randn(size(mag1sigfit));
            Gmc_phs = angle(Gfit) + phs1sigfit .* randn(size(phs1sigfit));

            Gmc = Gmc_mag .* exp(j * Gmc_phs);
        
         end
      
         [phic_resp_new,phic_ref_new,Freqnew,Dampnew,kterms] = ...
            rpfit(Gmc,nfitbin,1,nfitbin,wifit,nmodfit,exterms,respDOF,refDOF);
      
         Freqallnew(:,k) = Freqnew;
         Dampallnew(:,k) = Dampnew;
         Phicrespallnew(:,k,:) = phic_resp_new;
         Phicrefallnew(:,k,:) = phic_ref_new;
         
      end
      %
      % Display Modal Analysis Statistics
      %   
   
      mcdisp  
   
   elseif istatmeth == 2,
   
      %
      %  BOOTSTRAP Analysis
      %
      
      %  Define the number of bootstrap ensembles
      
      nboot = size(G,3);
     
      %  Start the loop
   
      for k = 1:nstatiter,
         
         fprintf(1,'Bootstrap Run #%d of %d\n',k,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(Gfit(:,:,rrand),3);
         
         %  Run the rp fit
         
          [phic_resp_new,phic_ref_new,Freqnew,Dampnew,kterms] = ...
            rpfit(Gboot,nfitbin,1,nfitbin,wifit,nmodfit,exterms,respDOF,refDOF);
         
         %  Add the results of this run to "allnew" sets

         Freqallnew(:,k) = Freqnew;
         Dampallnew(:,k) = Dampnew;
         Phicrespallnew(:,k,:) = phic_resp_new;
         Phicrefallnew(:,k,:) = phic_ref_new;
         
         %  End the loop
         
      end
      
      
      %  Display the results
      
      fprintf(1,'\n\nBootstrap Results:\n')
      
      for i = 1:size(Freqallnew,1),
         fprintf(1,'\nNew Mode #%d of %d\n',i,size(Freqallnew,1))
         fprintf(1,'\n%s\n','Run #  Freq(Hz)     Damp(%)   ')
         for j = 1:nstatiter,
            fprintf(1,'%-7d%10.3e%13.3e\n',j,Freqallnew(i,j),Dampallnew(i,j))
         end
      end
      
      if input('Accept these values (1=y)?'),
         disp('New values accepted')
         if exist('Freq')~=1,Freq=[];end;
         if exist('Damp')~=1,Damp=[];end;
         Freq = [Freq; mean(Freqallnew,2)];
         Damp = [Damp; mean(Dampallnew,2)];
         if exist('phic_resp')~=1,
            phic_resp = [];
            phic_ref = [];
         end
         for i = 1:size(Phicrespallnew,3),
            phic_resp = [phic_resp , mean(Phicrespallnew(:,:,i).'').''];
            phic_ref  = [phic_ref  , mean(Phicrefallnew(:,:,i).'').''];
         end
         if exist('Freqall')~=1,
            Freqall=[];
            Dampall=[];
         end
         Freqall = [Freqall; Freqallnew];
         Dampall = [Dampall; Dampallnew];
         if exist('Phicrespall')~=1, 
            Phicrespall = []; 
            Phicrefall = []; 
            nmodtemp=0;
         else,  
            nmodtemp = size(Phicrespall,3); 
         end
         nmodnew = size(Phicrespallnew,3);
         Phicrespall(:,:,[nmodtemp+1 : nmodtemp + nmodnew]) = Phicrespallnew;
         Phicrefall(:,:,[nmodtemp+1 : nmodtemp + nmodnew]) = Phicrefallnew;
      else
         disp('New values rejected')
      end
      
   end
   
else
   
   %
   %  Regular Rational Polynomial fit analysis
   %
   
   [phic_resp_new,phic_ref_new,Freqnew,Dampnew,kterms] = ...
         rpfit(Gfit(:,:,1),nfitbin,1,nfitbin,wifit,nmodfit,exterms,respDOF,refDOF);
      
   %
   % Calculate the frequencies and mode shapes
   %
   if exist('Freq')~=1,Freq=[];end;
   if exist('Damp')~=1,Damp=[];end;
   if exist('phic_resp')~=1,phic_resp=[];end;
   if exist('phic_ref')~=1,phic_ref=[];end;
   if exist('icountfull')~=1,icountfull=[];end;

   ifull=length(Freq);
   icountnew=[ifull+1:1:ifull+nmodfit]';

   textfull=RPtext(icountfull,Freq,Damp);
   set(hedit1,'string',textfull);

   textnew=RPtext(icountnew,Freqnew,Dampnew);
   set(hedit2,'string',textnew);

   set(figRPresults,'visible','on')

end

disp('Finished Rational Polynomial Curve Fit.')

end
