disp(' ');
disp(' powertrans_from_modes.m   ver 1.0   January 28, 2010 ');
disp(' by Tom Irvine ');
%
disp(' ');
disp(' This program calculates a power transmissibility function (G^2/G^2) ');
disp(' for a given pair of degree-of-freedoms in a system based on the     ');
disp(' mode shapes, natural frequencies, and damping ratios.               ');
disp(' The transmissibility function is referenced to a force input node.  ');
%
clear freq;
clear fnv;
clear dampv;
clear QE;
clear H;
clear HM;
clear max;
clear TM;
clear PT;
clear PTM;
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=3;
%
if(imethod==1)
    disp(' ');
        massm = input(' Enter the mass matrix name:  ');
        massm
    disp(' ');
    disp(' Divide mass by 386? ');
    disp(' 1=yes 2=no ');
    idm=input(' ');
    if(idm==1)
        massm=massm/386;
        massm
    end
    m=massm;
    disp(' ');
    stiffness = input(' Enter the stiffness matrix name:  ');
    stiffness
%
    [ModeShapes,Eigenvalues]=eig(stiffness,massm);
%
    omega = sqrt(Eigenvalues);
    num=max(size(massm));
%
    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 massm-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(' Enter the modal damping ratio vector ');
dampv=input(' ');
%
disp(' ');
disp(' Enter the frequency step ');
df=input(' ');
%
disp(' ');
disp(' Enter the maximum frequency ');
maxf=input(' ');
%
disp(' ');
disp(' Enter the minimum frequency for plots');
minf=input(' ');
%
nf=round(maxf/df);
for(i=1:nf)
    freq(i)=i*df;
end
%
j=sqrt(-1);
%
clear omn;
omn=tpi*fnv;
H=zeros(num,num,nf);
%
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)/(omn(r)^2);
%        
                if(iam==2)
                    term=term*tpi*freq(s);                   
                end
                if(iam==3)
                    term=term*(tpi*freq(s))^2;                    
                end 
                H(i,k,s)=H(i,k,s)+term;
            end
        end
    end
end
%
HM=abs(H);
%
%% 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
%
disp(' ');
disp(' Enter the force input node number ');
force_node=input(' ');
%
disp(' ');
disp(' Enter the base node dof number ');
base_node=input(' ');
%
disp(' ');
disp(' Enter the response node dof number ');
response_node=input(' ');
%
for(s=1:nf) % frequency loop  
    TM(s)=HM(response_node,force_node,s)/HM(base_node,force_node,s);
    PT(s)=(TM(s))^2;
end
%
PTM=[freq' PT'];
figure(1);
plot(freq,PT);
out1=sprintf('Power Transmissibility (Mass %d/Mass %d) with Force at Mass %d ',response_node,base_node,force_node);
title(out1);
xlabel('Frequency(Hz)');
ylabel('Trans(G^2/G^2)');
grid on;
ymax=max(PT);
ymin=min(PT);
      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.001)
          ymin=ymax*0.001;
      end
axis([minf,maxf,ymin,ymax]);
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
%
disp(' ');
disp(' The power transmissibility matix name is PTM.');
%%
choice=input(' Export as txt file?  1=yes  2=no  ' );
disp(' ')
%
if choice == 1 
       [writefname, writepname] = uiputfile('*','Save data as');
	   writepfname = fullfile(writepname, writefname);
	   writedata = [PTM(:,1) PTM(:,2) ];
	   fid = fopen(writepfname,'w');
	   fprintf(fid,'  %g  %g\n',writedata');
	   fclose(fid);
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
disp(' ');
disp(' Multiply the power transmissibility function by an input PSD?');
disp('  1=yes  2=no ');
ium=input(' ');
if(ium==1)
%
    yone=double(PTM(:,2));
    xone=double(PTM(:,1));
%
    out1=sprintf(' Select file input method for the PSD at node dof %d',base_node);
    disp(out1);
    disp('   1=external ASCII file ');
    disp('   2=file preloaded into Matlab ');
    file_choice = input('');
%
    if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        THM=THM';
    else
        THM = input(' Enter the matrix name:  ');
    end
%
    clear ytwo;
    clear xtwo;
    clear a1;
    clear f1;
    ytwo=double(THM(:,2));
    xtwo=double(THM(:,1));
%
    f1=xone;
    a1=yone;
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Q2 = real_mult_intlog(xone,ytwo,df);
    clear xi;
    clear yi;
    clear xL;
    clear yL;
    clear xiL;
    clear yiL;
%
    xi=min(xone):df:max(xone); 
%
    xL=log10(xtwo);
    yL=log10(ytwo);
    xiL=log10(xi);
%
    disp('ref 1');
    yiL = interp1(xL,yL,xiL);
%
    for(i=1:max(size(xi)))
        xi(i)=10^xiL(i);
        yi(i)=10^yiL(i);
    end
%   
    disp('ref 2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear a2;
    clear f2;
    ijk=1;
    for(i=1:max(size(yi)))
       if(yi(i)>=-1.0e-30 && yi(i)<=1.0e+30)
          f2(ijk)=xi(i); 
          a2(ijk)=yi(i);
          ijk=ijk+1;
       end
    end
%
    df2=df/2;
%
    clear ab;
    clear ff;
%
    disp('ref 3');
%
    ijk=1;
    for(i=1:max(size(f1)))
        for(j=1:max(size(f2)))
            if(abs(f1(i)-f2(j))<df2)
                ab(ijk)=a1(i)*a2(j);
                ff(ijk)=(f1(i)+f2(j))/2.;
                ijk=ijk+1;
                break;
            end
        end
    end    
%
    area=df*sum(ab);
    grms=sqrt(area);
%
    area=df*sum(a2);
    grms_in=sqrt(area);    
%
    disp(' See figure(2) ');
%
    figure(2);
    plot(f2,a2,'r',ff,ab,'b');
    out1=sprintf('PSD  dof %d  overall level=%6.3g GRMS',response_node,grms);
    title(out1);
    out1=sprintf(' Input %6.3g GRMS',grms_in);  
    out2=sprintf('Response %6.3g GRMS',grms);   
    legend (out1,out2,2);      
    ylabel('Accel(G^2/Hz)');
    xlabel('Frequency(Hz)');
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
    grid('on');
%
    disp(' ');
    choice=input(' Export output PSD as txt file?  1=yes  2=no  ' );
    disp(' ')
%
    if choice == 1 
       [writefname, writepname] = uiputfile('*','Save data as');
	   writepfname = fullfile(writepname, writefname);
	   writedata = [ff' ab' ];
	   fid = fopen(writepfname,'w');
	   fprintf(fid,'  %g  %g\n',writedata');
	   fclose(fid);
    end
end