%
%   mdof_frf_curvefit_phase1.m  ver 1.1  January 22, 2014
%
function[fr,dr,AXr,errormax]=...
    mdof_frf_curvefit_phase1(num,NM,maxa,maxf,c_r,c_i,c_r_rec,c_i_rec,...
                             fr,dr,AXr,MFRF_R,MFRF_I,tlabel,yl,errormax,...
                             nh,fest,ntrials,f,Hsum,i1,i2,frf_r,frf_i)
%
    tpi=2*pi;
    k=1;
    for ijk=1:ntrials
%      
       nv=NM(k);
       k=k+1;
       if(k>num)
           k=1;
       end
%
       for i=1:num
           c_r(:,i)=c_r_rec(:,i);
           c_i(:,i)=c_i_rec(:,i);
       end
%
       if(ijk<round(0.2*ntrials))
           ft=fest(nv)*(0.98+0.04*rand());
           damp=0.15*rand();           
       else    
          if(ijk<20 || rand()<0.6)
                ft=fest(nv)*(0.95+0.10*rand());
                damp=0.2*rand();
          else 
              if(rand()<0.5)
                ft=fr(nv)*(0.98+0.04*rand());
                damp=dr(nv)*(0.98+0.04*rand());
              else
                ft=fr(nv)*(0.99+0.02*rand());
                damp=dr(nv)*(0.99+0.02*rand());
              end  
          end
       end  
%
       if(rand()<0.03)
           ft=f(i1)+( f(i2)-f(i1)  )*rand();
       end
%
       omegan=tpi*ft;      
%
       error=0;
%
       for j=i1:i2
          omega=tpi*f(j);
          rho=omega/omegan;
          c_r(j,nv)=MFRF_R(omega,rho,damp);
          c_i(j,nv)=MFRF_I(omega,rho,damp);
       end
%
       Ht=c_r(:,nv)+(1i)*c_i(:,nv);
%
       if(ijk<20 || rand()<0.4)
            AX=maxa/max(abs(Ht))*(0.1+0.9*rand());
       else
            if(rand()<0.5)
                AX=AXr(nv)*(0.95+0.1*rand());       
            else    
                AX=AXr(nv)*(0.99+0.02*rand());
            end
       end    
%
       c_r(:,nv)=AX*c_r(:,nv);
       c_i(:,nv)=AX*c_i(:,nv);
%
       for j=i1:i2
          sum_cr=sum(c_r(j,:));
          sum_ci=sum(c_i(j,:));
%          
          bbb=abs(frf_r(j)-sum_cr)+abs(frf_i(j)-sum_ci);
          error=error+bbb;
       end
%      
       if(error<errormax)
           errormax=error;
           fr(nv)=ft;
           dr(nv)=damp;      
           AXr(nv)=AX;
           c_r_rec(:,nv)=c_r(:,nv);
           c_i_rec(:,nv)=c_i(:,nv);
%%
           out1=sprintf('%d  %8.2e ',ijk,errormax);
           disp(out1);
           for i=1:num
               out1=sprintf('   %8.3e  %8.2f  %8.3f   ',AXr(i),fr(i),dr(i));
               disp(out1);
           end 
%%
           for ij=i1:i2
               Hsum(ij)=sum(c_r_rec(ij,:))+(1i)*sum(c_i_rec(ij,:));
           end    
%%
           fig_num=2;
           figure(fig_num);
           subplot(2,1,1);
           plot(f(i1:i2),frf_r(i1:i2),f(i1:i2),real(Hsum(i1:i2)));
           grid on;
           legend('Measured Data','Curve-fit',1);
           xlabel('Frequency (Hz)');
           ylabel(yl);
           out1=sprintf('%s Real Frequency Response Function',tlabel);
           title(out1);
%
           subplot(2,1,2);
           plot(f(i1:i2),frf_i(i1:i2),f(i1:i2),imag(Hsum(i1:i2)));
           grid on;
           legend('Measured Data','Curve-fit',1);
           xlabel('Frequency (Hz)');
           ylabel(yl);
           out1=sprintf('%s Imaginary Frequency Response Function',tlabel);
           title(out1); 
%%
%%
       end
    end                                                             