disp(' ');
disp(' transfer_from_modes.m   ver 1.7   December 17, 2011 ');
disp(' by Tom Irvine ');
%
disp(' ');
disp(' This program calculates a transfer function for each   ');
disp(' degree-of-freedom in a system based on the mode shapes, ');
disp(' natural frequencies, and damping ratios. ');
%
clear freq;
clear fnv;
clear dampv;
clear QE;
clear H;
clear HM;
clear HP;
clear max;
clear PHA;
clear PPP;
tpi=2*pi;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select input method ');
disp('  1=mass & stiffness matrices ');
disp('  2=natural frequencies and mass-normalized eigenvectors ');
%
imethod=input(' ');
%
disp(' ');
disp(' Select output metric ');
disp('  1=displacement/force ');
disp('  2=velocity/force ');
disp('  3=acceleration/force ');
%
iam=input(' ');
%
if(imethod==1)
    disp(' ');
        mass = input(' Enter the mass matrix name:  ');
        mass
    disp(' ');
    disp(' Divide mass by 386? ');
    disp(' 1=yes 2=no ');
    idm=input(' ');
    if(idm==1)
        mass=mass/386;
        mass
    end
    m=mass;
    disp(' ');
    stiffness = input(' Enter the stiffness matrix name:  ');
    stiffness
%
    [ModeShapes,Eigenvalues]=eig(stiffness,mass);
%
    omega = sqrt(Eigenvalues);
    num=max(size(mass));
%
    for(i=1:num)
        fnv(i)=omega(i,i)/tpi;
    end
%
    disp(' ');
    disp('  Natural Frequencies (Hz) ');
    for(i=1:num)
       out1=sprintf('  %8.4g ',fnv(i));
       disp(out1);
    end
%
    MST=ModeShapes';
    temp=zeros(num,num);
    temp=m*ModeShapes;
    QTMQ=MST*temp;
%   
    clear scale;
    for(i=1:num)
        scale(i)=1./sqrt(QTMQ(i,i));
        ModeShapes(:,i) = ModeShapes(:,i)*scale(i);    
    end
%
    disp(' ');
    disp('  Modes Shapes (column format)');
    QE=ModeShapes
else
    disp(' ');
    disp(' Enter the name of the mass-normalized eigenvector matrix ');
    QE=input(' ');
    num=max(size(QE));
%
    disp(' ');
    disp(' Enter the natural frequency (Hz) vector or matrix ');
    clear fff;
    fff=input(' ');
    sz=size(fff);
    if(sz(1)==sz(2))
        for(i=1:sz(1))
            fnv(i)=fff(i,i);
        end
    else
        fnv=fff;
    end
end
%
disp(' ');
disp(' Select damping input method ');
disp(' 1=constant for all modes   2=vector ');
dm=input(' ');
if(dm==1)
   disp(' ');
   disp(' Enter the modal damping ratio ');
   mdr=input(' ');
   dampv=mdr*ones(max(size(fnv)),1);
else
   disp(' ');
   disp(' Enter the modal damping ratio vector ');
   dampv=input(' ');
end
%
disp(' ');
disp(' Enter the frequency step ');
df=input(' ');
%
disp(' ');
disp(' Enter the minimum frequency for plots');
minf=input(' ');
%
disp(' ');
disp(' Enter the maximum frequency for plots');
maxf=input(' ');
%
if(maxf<=minf)
    maxf=100*minf;
end
%
nf=round(maxf/df);
clear omega;
%
for(i=1:nf)
    freq(i)=i*df;
    omega(i)=2*pi*freq(i);
    omega2(i)=(omega(i))^2;
end
%
j=sqrt(-1);
%
clear omn;
omn=tpi*fnv;
omn2=omn.*omn;
H=zeros(num,num,nf);
%
while(1)
disp(' ');
i=input(' Enter first  dof of interest ');
k=input(' Enter second dof of interest ');
for(s=1:nf) % frequency loop
%
%%    for(i=1:num)
%%        for(k=1:num)
            for(r=1:num)
                rho=freq(s)/fnv(r);
                den=1-rho^2+j*2*dampv(r)*rho;
%                
                term=(QE(i,r)*QE(k,r)/den)/omn2(r);
%        
                if(iam==2)
                    term=term*omega(s);                   
                end
                if(iam==3)
                    term=term*omega2(s);                    
                end 
                H(i,k,s)=H(i,k,s)+term;
            end
%%        end
%%    end
end
HM=abs(H);
HP=-atan2(imag(H),real(H))*180/pi;
%
disp(' ');
disp(' Enter units: ');
if(iam==1)
    disp(' 1= inches, lbf ');
    disp(' 2= meters, N ');
end
if(iam==2)
    disp(' 1= inches/sec, lbf ');
    disp(' 2= meters/sec, N ');
end
if(iam==3)
    disp(' 1= G, lbf ');
    disp(' 2= meters/sec^2, N ');
end
iu=input(' ');
%
if(iam==3 && iu==1)
   HM=HM/386;
end
%
ijk=0;
%%% for(i=1:num)
%%%     for(k=1:num)
      ijk=ijk+1;
      clear PPP;
      for(ia=1:nf)
        PPP(ia)=HM(i,k,ia);
        PHA(ia)=HP(i,k,ia);
%
        while(PHA(ia)<0)
            PHA(ia)=PHA(ia)+360.;
        end
%
        ymax=max(PPP);
        ymin=min(PPP);
      end
%
      for(ia=1:20)
          if(ymax<0.9*(10^(ia-10)))
              ymax=10^(ia-10);
              break;
          end
      end
      for(ia=1:20)
          if(ymin<10^(ia-10))
              ymin=10^(ia-11);
              break;
          end
      end    
%
      if(ymin<ymax*0.0001)
          ymin=ymax*0.0001;
      end
%      
      figure(2*ijk-1);
%
      sz=size(freq);
      if(sz(2)>sz(1))
          freq=freq';
      end
%
      sz=size(PPP);
      if(sz(2)>sz(1))
          PPP=PPP';
      end 
%
      str=sprintf('H_%d_%d',i,k);
      varname = genvarname(str);
      eval([varname ' = [freq PPP]']);
      out1=sprintf(' Matlab array name:  %s ',varname);
      disp(out1);
%
      plot(freq,PPP);
      grid on;
      xlabel('Frequency(Hz)');
%      
      if(iam==1)
        if(iu==1)
            ylabel('Disp/Force (in/lbf)');
            out1=sprintf('\n maximum = %8.4g (in/lbf) \n',max(PPP));  
        else
            ylabel('Disp/Force (m/N) ');
            out1=sprintf('\n maximum = %8.4g  (m/N) \n',max(PPP));
        end      
      end
%
      if(iam==2)
        if(iu==1)
            ylabel('Vel/Force (in/sec/lbf)');
            out1=sprintf('\n maximum = %8.4g (in/sec/lbf) \n',max(PPP));        
        else
            ylabel('Vel/Force (m/sec/N) ');
            out1=sprintf('\n maximum = %8.4g (m/sec/N) \n',max(PPP));      
        end
      end
%
      if(iam==3)
        if(iu==1)
            ylabel('Accel/Force (G/lbf)');
            out1=sprintf('\n maximum = %8.4g (G/lbf) \n',max(PPP));   
        else
            ylabel('Accel/Force (m/sec^2/N) '); 
            out1=sprintf('\n maximum = %8.4g (m/sec^2/N) \n',max(PPP));        
        end
      end
%
      disp(out1);  
%
      out1=sprintf('Transfer Function Magnitude H%d%d',i,k);
      title(out1);
      axis([minf,maxf,ymin,ymax]);
      set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
      figure(2*ijk);
      plot(freq,PHA);
      out1=sprintf('Transfer Function Phase H%d%d',i,k);
      title(out1);     
      grid on;
      xlabel('Frequency(Hz)');
      ylabel('Phase(deg)');    
%    
%%%     end
%%% end
%
disp(' ');
disp(' Perform calculation for another pair ? ');
disp(' 1=yes  2=no ');
inc=input(' ');
    if(inc==2)
        break;
    end
end