disp(' ');
disp(' joint_receptance.m  ver 1.0  January 12, 2014 '); 
disp(' ');
disp(' by Tom Irvine ');
disp('  ');
%
disp(' This script calculates the joint receptance for two subsystems ');
disp(' which are to be connected at a common dof. ');
%
clear H;
clear HSA;
clear HSB;
clear FRF_p;
clear a;
%
close all;
%
fig_num=1;
%
disp(' ');
disp(' Each receptance function must have the same dimension and ');
disp(' have the same frequency column. ');
disp(' ');
%
disp(' Each receptance function must have two columns:  ');
disp('   freq (Hz) & complex ');
disp('  ');
%
HSA=input(' Enter the receptance for subsystem A:  ');
%
HSB=input(' Enter the receptance for subsystem B:  ');
%
disp(' ');
disp(' Select units:  1=English  2=metric ');
iu=input(' ');
%
H=HSA;
%
sza=size(HSA);
szb=size(HSB);
%
if(sza(1)==szb(1))
%
    n=sza(1);
    FRF_p=zeros(n,1);
%
    for i=1:n
        H(i,2)=HSA(i,2)*HSB(i,2)/( HSA(i,2) + HSB(i,2)  );
        FRF_p(i)=atan2(imag(H(i,2)),real(H(i,2)))*180/pi;
    end
%
%
    ff=H(:,1);
    disp(' ');
    fmin=input(' Enter minimum plot frequency (Hz) ');
    fmax=input(' Enter maximum plot frequency (Hz) ');
    ymax=max(abs(H(:,2)));
    ymin=min(abs(H(:,2)));
    ymax=10^(ceil(log10(ymax)));
    ymin=10^(floor(log10(ymin)));    
%   
%
    if(ymin<ymax*0.0001)
        ymin=ymax*0.0001;
    end
%
    figure(fig_num);
    fig_num=fig_num+1;
%
    subplot(3,1,1);
    plot(ff,FRF_p);
    out1=sprintf('Joint Receptance Function Magnitude & Phase ');
    title(out1);   
    grid on;
    ylabel('Phase (deg)');
    axis([fmin,fmax,-180,180]);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
                    'YScale','lin','ytick',[-180,-90,0,90,180]);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    if(max(FRF_p)<=0.)
%
        axis([fmin,fmax,-180,0]);
        set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
                    'YScale','lin','ytick',[-180,-90,0]);
    end  
%
    if(min(FRF_p)>=-90. && max(FRF_p)<90.)
%
        axis([fmin,fmax,-90,90]);
        set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
                    'YScale','lin','ytick',[-90,0,90]);
    end 
%
    if(min(FRF_p)>=0.)
%
        axis([fmin,fmax,0,180]);
        set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
                    'YScale','lin','ytick',[0,90,180]);
    end 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    subplot(3,1,[2 3]);
    plot(freq,PPP);
    grid on;
    a=[ff abs(H(:,2))]; 
    [C,I]=max(a);
    xmax=a(I(2),2);
    max_freq=a(I(2),1);    
%
    if(iu==1)
            ylabel('Disp/Force (in/lbf)');
            out1=sprintf('\n maximum = %8.4g (in/lbf) at %8.4g Hz \n',xmax,max_freq);  
    else
            ylabel('Disp/Force (m/N) ');
            out1=sprintf('\n maximum = %8.4g  (m/N) at %8.4g Hz \n',xmax,max_freq);
    end 
%
    xlabel('Frequency(Hz)');
    axis([fmin,fmax,ymin,ymax]);
    set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
    disp(out1);    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    disp('  ');
    disp(' The joint receptance function is:  H  - freq(Hz) complex');
    disp('  ');
else
    disp(' ');
    disp(' Error:  subsystems have different sizes. ');
end