%
%  plot_magnitude_phase_two.m  ver 1.1  January 7, 2014
%
%  by Tom Irvine
%
disp(' This script plots the respective magnitude & phase of two ');
disp(' complex frequency response functions.');
clear THM;
clear THM1;
clear THM2;
clear ff;
clear FRF_p1;
clear FRF_m1;
%
close all;
fig_num=1;
%
disp(' ');
disp(' The input array must have two columns:  freq(Hz) complex amplitude ');
disp(' ');
THM1=input(' Enter complex array name 1:  ');
label1=input(' Enter curve name for legend ','s');
%
disp(' ');
THM2=input(' Enter complex array name 2:  ');
label2=input(' Enter curve name for legend ','s');
%
disp(' ');
fmin=input(' Enter minimum plot frequency (Hz)  ');
fmax=input(' Enter maximum plot frequency (Hz)  ');
%
disp(' ');
t_string=input(' Enter plot title: ','s');
%
sz=size(THM1);
%
n=sz(1);
ff1=zeros(n,1);
FRF_m1=zeros(n,1);
FRF_p1=zeros(n,1);
%
for i=1:n
    ff1(i)=THM1(i,1);
    FRF_m1(i)=abs(THM1(i,2));
        aa=real(THM1(i,2));
        bb=imag(THM1(i,2));
%
    FRF_p1(i)=(180/pi)*atan2(bb,aa);
end 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
sz=size(THM2);
n=sz(1);
ff2=zeros(n,1);
FRF_m2=zeros(n,1);
FRF_p2=zeros(n,1);
%
for i=1:n
    ff2(i)=THM2(i,1);
    FRF_m2(i)=abs(THM2(i,2));
        aa=real(THM2(i,2));
        bb=imag(THM2(i,2));
%
    FRF_p2(i)=(180/pi)*atan2(bb,aa);
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(fig_num);
fig_num=fig_num+1;
%
subplot(3,1,1);
plot(ff1,FRF_p1,ff2,FRF_p2);
title(t_string);
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_p1)<=0. && max(FRF_p2)<=0.)
%
axis([fmin,fmax,-180,0]);
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
                    'YScale','lin','ytick',[-180,-90,0]);
end  
%
if(min(FRF_p1)>=-90. && max(FRF_p1)<=90. && min(FRF_p2)>=-90. && max(FRF_p2)<=90.)
%
axis([fmin,fmax,-90,90]);
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
                    'YScale','lin','ytick',[-90,0,90]);
end 
%
if(min(FRF_p1)>=0. && max(FRF_p2)>=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(ff1,FRF_m1,ff2,FRF_m2);
grid on;
xlabel('Frequency(Hz)');
ylabel('Magnitude');
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log',...
         'YScale','log');
%
legend (label1,label2,1);
maxF=max(FRF_m1);
minF=min(FRF_m1);
if(max(FRF_m2)>max(FRF_m1))
    maxF=max(FRF_m2);    
end   
if(min(FRF_m2)<min(FRF_m1))
    minF=min(FRF_m2);    
end 
%
ymax= 10^(ceil(log10(maxF)+0.1));
ymin= 10^(floor(log10(minF)-0.1));
%
if(ymin<(ymax/1.0e+05))
    ymin=(ymax/1.0e+05);
end    
%
axis([fmin,fmax,ymin,ymax]);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%