echo off
%	tdof_prop_damped.m		proportionally damped tdof model
%	Calculates poles and zeros of proportionally damped tdof system.  
%	Plots initial condition responses for modes 2 and 3 in physical 
%	and principal coordinate systems.

%	solve for the eigenvalues of the undamped system model

	tdofss_eig;

	subplot(1,1,1);

%	now normalize the undamped eigenvectors with respect to the position of
%   mass 1, which will be set to 1.0 - for plotting of undamped Argand diagram

	for  cnt = 1:length(omegad)

		xmon1(:,cnt) = xmo(:,cnt)/xmo(1,cnt);

	end

	xmon1

%	input proportional damping for equations in principal coordinate system

   	zeta = input('input value for zeta, default = 0.02, 2% of critical ...  ');
	
	if  (isempty(zeta))
		zeta = 0.02;
	else
	end

% 	setup proportionally damped state-space system matrix in principal coordinates

    a_ss = [  0    1    0      0        0       0
         	 		  0    0    0      0        0       0
             		  0    0    0      1        0       0
	         		  0    0  -w2^2 -2*zeta*w2  0       0
             		  0    0    0      0        0       1
             		  0    0    0      0      -w3^2 -2*zeta*w3];

%	solve for the eigenvalues of the system matrix with proportional damping

	[xmp,omegap] = eig(a_ss);

%	take the diagonal elements of the generalized eigenvalue matrix omegap

	omegapd = diag(omegap);

%	now reorder the eigenvalues and eigenvectors from low to high frequency, 
%	keeping track of how the eigenvalues are ordered in reorder the 
%	eigenvectors to match, using indexhz

	[omegaporder,indexhz] = sort(abs(imag(omegapd)));

	for  cnt = 1:length(omegapd)

		omegapo(cnt,1) = omegapd(indexhz(cnt));	%  reorder eigenvalues

		xmpo(:,cnt) = xmp(:,indexhz(cnt));		%  reorder eigenvector columns

	end

%	now calculate the magnitude and phase angle of each of the eigenvector
%	entries

	for  row = 1:length(omegapd)

		for  col = 1:length(omegapd)

			xmpomag(row,col) = abs(xmpo(row,col));

			xmpoang(row,col) = (180/pi)*angle(xmpo(row,col));

		end

	end

	omegapo

	xmpo

	xmpomag

	xmpoang

%	calculate the percentage of critical damping for each mode

	zeta1 = 0
	
	theta2 = atan(real(omegapo(3))/imag(omegapo(3)));
	zeta2 = abs(sin(theta2))

	theta3 = atan(real(omegapo(5))/imag(omegapo(5)));
	zeta3 = abs(sin(theta3))

	plot(omegap,'k*')
	grid on
	axis([-3 1 -2 2])
	axis('square')
	title('Proportionally Damped Eigenvalues')
	xlabel('real')
	ylabel('imaginary')
	text(real(omegapo(3))-1,imag(omegapo(3))+0.1,['zeta = ',num2str(zeta2)])	
	text(real(omegapo(5))-1,imag(omegapo(5))+0.1,['zeta = ',num2str(zeta3)])	

	disp('execution paused to display figure, "enter" to continue'); pause

%	calculate the motions of the three masses for all three modes - damped case

	t = 0:.12:15;

	sigma11 = real(omegapo(1));		%  sigma for first eigenvalue for mode 1
	omegap11 = imag(omegapo(1));   	%  omegap for first eigenvalue for mode 1

	sigma12 = real(omegapo(2));		%  sigma for second eigenvalue for mode 1
	omegap12 = imag(omegapo(2));   	%  omegap for second eigenvalue for mode 1

	sigma21 = real(omegapo(3));		%  sigma for first eigenvalue for mode 2
	omegap21 = imag(omegapo(3));   	%  omegap for first eigenvalue for mode 2

	sigma22 = real(omegapo(4));		%  sigma for second eigenvalue for mode 2
	omegap22 = imag(omegapo(4));   	%  omegap for second eigenvalue for mode 2

	sigma31 = real(omegapo(5));		%  sigma for first eigenvalue for mode 3
	omegap31 = imag(omegapo(5));   	%  omegap for first eigenvalue for mode 3

	sigma32 = real(omegapo(6));		%  sigma for second eigenvalue for mode 3
	omegap32 = imag(omegapo(6));  	%  omegap for second eigenvalue for mode 3

%	displacements of mode 1 in principal coordinates

	zp111 = exp(sigma11*t).*(exp(i*omegap11*t)*xmpo(1,1));		%  mass 1
	zp112 = exp(sigma12*t).*(exp(i*omegap12*t)*xmpo(1,2));		%  mass 1

%	displacements of mode 2 in principal coordinates

	zp221 = exp(sigma21*t).*(exp(i*omegap21*t)*xmpo(3,3));		%  mass 2
	zp222 = exp(sigma22*t).*(exp(i*omegap22*t)*xmpo(3,4));		%  mass 2

%	displacements of mode 3 in principal coordinates

	zp331 = exp(sigma31*t).*(exp(i*omegap31*t)*xmpo(5,5));		%  mass 3
	zp332 = exp(sigma32*t).*(exp(i*omegap32*t)*xmpo(5,6));		%  mass 3

%	calculate the displacements of each mass for mode 2
%	define matrix of displacements vs time for each eigenvector

	z221 = [zeros(1,length(t))
			zp221
			zeros(1,length(t))];

	z222 = [zeros(1,length(t))
			zp222
			zeros(1,length(t))];

%	back-transform from principal to physical coordinates

	zmode21 = xn*z221;

	zmode22 = xn*z222;

	z1mode21 = zmode21(1,:);

	z2mode21 = zmode21(2,:);

	z3mode21 = zmode21(3,:);

	z1mode22 = zmode22(1,:);

	z2mode22 = zmode22(2,:);

	z3mode22 = zmode22(3,:);

%	calculate the displacements of each mass for mode 3
%	define matrix of displacements vs time for each eigenvector

	z331 = [zeros(1,length(t))
			zeros(1,length(t))
			zp331];

	z332 = [zeros(1,length(t))
			zeros(1,length(t))
			zp332];

	zmode31 = xn*z331;

	zmode32 = xn*z332;

	z1mode31 = zmode31(1,:);

	z2mode31 = zmode31(2,:);

	z3mode31 = zmode31(3,:);

	z1mode32 = zmode32(1,:);

	z2mode32 = zmode32(2,:);

	z3mode32 = zmode32(3,:);

%	plot principal displacements of mode 2

	plot(t,real(zp221),'k-',t,real(zp222),'k+-',t,imag(zp221),'k.-',t,imag(zp222),'ko-')
	title('principal real and imag disp for mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	plot physical disp of masses for mode 2

	plot(t,real(z1mode21),'k-',t,real(z1mode22),'k+-',t,imag(z1mode21),'k.-',t,imag(z1mode22),'ko-')
	title('physical real and imag disp for mass 1, mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z2mode21),'k-',t,real(z2mode22),'k+-',t,imag(z2mode21),'k.-',t,imag(z2mode22),'ko-')
	title('physical real and imag disp for mass 2, mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z3mode21),'k-',t,real(z3mode22),'k+-',t,imag(z3mode21),'k.-',t,imag(z3mode22),'ko-')
	title('physical real and imag disp for mass 3, mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z1mode21+z1mode22),'k-',t,real(z2mode21+z2mode22),'k+-',t,real(z3mode21+z3mode22),'k.-')
	title('physical disp z1, z2, z3 mode 2')
	legend('mass 1','mass 2','mass 3')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	plot subplots for notes

%	plot principal disp of mode 2

	subplot(3,2,1)
	plot(t,real(zp221),'k-',t,real(zp222),'k+-',t,imag(zp221),'k.-',t,imag(zp222),'ko-')
	title('principal disp for mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -1 1])
	grid on

%	plot physical disp of masses for mode 2

	subplot(3,2,3)
	plot(t,real(z1mode21),'k-',t,real(z1mode22),'k+-',t,imag(z1mode21),'k.-',t,imag(z1mode22),'ko-')
	title('physical real and imag disp for mass 1, mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	subplot(3,2,4)
	plot(t,real(z2mode21),'k-',t,real(z2mode22),'k+-',t,imag(z2mode21),'k.-',t,imag(z2mode22),'ko-')
	title('physical real and imag disp for mass 2, mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	subplot(3,2,5)
	plot(t,real(z3mode21),'k-',t,real(z3mode22),'k+-',t,imag(z3mode21),'k.-',t,imag(z3mode22),'ko-')
	title('physical real and imag disp for mass 3, mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	subplot(3,2,6)
	plot(t,real(z1mode21+z1mode22),'k+-',t,real(z2mode21+z2mode22),'k.-',t,real(z3mode21+z3mode22),'ko-')
	title('physical disp for z1, z2, z3 mode 2')
	legend('mass 1','mass 2','mass 3')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	subplot(1,1,1)

%	plot principal disp of mode 3	

	plot(t,real(zp331),'k-',t,real(zp332),'k+-',t,imag(zp331),'k.-',t,imag(zp332),'ko-')
	title('principal disp for mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	plot physical disp of masses for mode 3

	plot(t,real(z1mode31),'k-',t,real(z1mode32),'k+-',t,imag(z1mode31),'k.-',t,imag(z1mode32),'ko-')
	title('physical real and imag disp for mass 1, mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z2mode31),'k-',t,real(z2mode32),'k+-',t,imag(z2mode31),'k.-',t,imag(z2mode32),'ko-')
	title('physical real and imag disp for mass 2, mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z3mode31),'k-',t,real(z3mode32),'k+-',t,imag(z3mode31),'k.-',t,imag(z3mode32),'ko-')
	title('physical real and imag disp for mass 3, mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z1mode31+z1mode32),'k-',t,real(z2mode31+z2mode32),'k+-',t,real(z3mode31+z3mode32),'k.-')
	title('physical disp for z1, z2, z3 mode 3')
	legend('mass 1','mass 2','mass 3')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	plot subplots for notes

%	plot principal disp of mode 3	

	subplot(3,2,1)
	plot(t,real(zp331),'k-',t,real(zp332),'k+-',t,imag(zp331),'k.-',t,imag(zp332),'ko-')
	title('principal disp for mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

%	plot physical disp of masses for mode 3

	subplot(3,2,3)
	plot(t,real(z1mode31),'k-',t,real(z1mode32),'k+-',t,imag(z1mode31),'k.-',t,imag(z1mode32),'ko-')
	title('physical real and imag disp for mass 1, mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	subplot(3,2,4)
	plot(t,real(z2mode31),'k-',t,real(z2mode32),'k+-',t,imag(z2mode31),'k.-',t,imag(z2mode32),'ko-')
	title('physical real and imag disp for mass 2, mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	subplot(3,2,5)
	plot(t,real(z3mode31),'k-',t,real(z3mode32),'k+-',t,imag(z3mode31),'k.-',t,imag(z3mode32),'ko-')
	title('physical real and imag disp for mass 3, mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -0.5 0.5])
	grid on

	subplot(3,2,6)
	plot(t,real(z1mode31+z1mode32),'k+-',t,real(z2mode31+z2mode32),'k.-',t,real(z3mode31+z3mode32),'ko-')
	title('physical disp for z1, z2, z3 mode 3')
	legend('mass 1','mass 2','mass 3')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause