echo off
%	tdofss_modal_time_ode45.m	 state-space modal form transient analysis
%   of tdof model, proportional damping, modal contribution plotting
%	Plots tdof transient responses for overall and individual modal 
%	contributions.  Calls the function files below, which define the 
%	state space system and individual modes.

	clf;

%	run tdofss_eig.m to provide eigenvalues and eigenvectors

	tdofss_eig;

	global  a_ss a1_ss a2_ss a3_ss b b1 b2 b3 u

%	note, this is the point where we would start if we had eigenvalue results from ANSYS,
%	using the eigenvalues and eigenvectors to define state space equations in 
%	principal coordinates

%	define damping ratio to be used for proportional damping in the state space equation
%	in principal coordinates
	
	zeta = input('input zeta, 0.02 = 2% of critical damping (default) ...  ');

   	if (isempty(zeta))
      	zeta = 0.02;
	else
   	end

% 	setup 6x6 state-space system matrix for all three modes in principal
%	coordinates, a_ss

    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];

%	setup three 2x2 state-space matrices, one for each individual mode

	a1_ss = a_ss(1:2,1:2);

	a2_ss = a_ss(3:4,3:4);

	a3_ss = a_ss(5:6,5:6);

%	transform the 3x1 force vector in physical coordinates to principal coordinates and
%	then insert the principal forces in the appropriate rows in the state-space
%	6x1 input matrix, padding with zeros as appropriate

	F = [1 0 -2]';

	Fp = xn'*F;

%	expand the force vectors in principal coordinates from 3x1 to 6x1, padding with zeros

	b = [0 Fp(1) 0 Fp(2) 0 Fp(3)]';	%  principal forces applied to all masses

	b1 = b(1:2);

	b2 = b(3:4);

	b3 = b(5:6);

%	the output matrix c is setup in one step, to allow the "bode" command to
%	output the desired physical coordinates directly without having to go 
%	through any intermediate steps.

% 	setup the output matrix for displacement transfer functions, each row
%   represents the position	outputs of mass 1, mass 2 and mass 3
%	velocities not included, so c is only 3x6 instead of 6x6

  	c = [xn(1,1) 0 xn(1,2) 0 xn(1,3) 0
		 xn(2,1) 0 xn(2,2) 0 xn(2,3) 0
		 xn(3,1) 0 xn(3,2) 0 xn(3,3) 0];

	c1 = c(:,1:2);

	c2 = c(:,3:4);

	c3 = c(:,5:6);

% 	define direct transmission matrix d

	d = 0;

%	transient response using the ode45 command

	u = 1;

	ttotal = input('Input total time for Simulation, default = 10 sec, ...  ');

   	if (isempty(ttotal))
     	ttotal = 10;
	else
   	end

	tspan = [0 ttotal];

%	calculate the initial conditions in principal coordinates using the inverse of the
%	normalized modal matrix

	x0phys = [0 -1 1]';			%	initial condition position, physical coord

	x0dphys = [-1 2 -2]';		%	initial condition velocity, physical coord

	x0 = inv(xn)*x0phys;

	x0d = inv(xn)*x0dphys;

%	create the initial condition state vector

	x0ss = [x0(1) x0d(1) x0(2) x0d(2) x0(3) x0d(3)];

	x0ss1 = x0ss(1:2);

	x0ss2 = x0ss(3:4);

	x0ss3 = x0ss(5:6);

%	use the ode45 non-stiff differential equation solver

	options = [];							%  no options specified

%  total response, principal coord, states are modes of vibration

  	[t,x] = ode45('tdofssmodalfun',tspan,x0ss,options);

%  mode 1 response, principal coord	

  	[t1,x1] = ode45('tdofssmodal1fun',tspan,x0ss1,options);

%  mode 2 response, principal coord

  	[t2,x2] = ode45('tdofssmodal2fun',tspan,x0ss2,options);

%  mode 3 response, principal coord

  	[t3,x3] = ode45('tdofssmodal3fun',tspan,x0ss3,options);

%  total response, physical coord

	z_ode = c*x';

%  mode 1 response, physical coord

	z_ode1 = c1*x1';

%  mode 2 response, physical coord

	z_ode2 = c2*x2';

%  mode 3 response, physical coord

	z_ode3 = c3*x3';

%	plot displacements in principal coordinates

	subplot(1,1,1);

	plot(t1,x1(:,1),'k+-',t2,x2(:,1),'kx-',t3,x3(:,1),'k-')
	title('Displacements in Principal Coordinate System, ode45')
	xlabel('Time, sec')
	ylabel('Displacements')
	legend('zp1','zp2','zp3',3)
	axis([0 10 -35 10])
	grid

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

	axis([0 1 -2 2]);

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

%	plot displacements in physical coordinates

	plot(t,z_ode(1,:),'k+-',t,z_ode(2,:),'kx-',t,z_ode(3,:),'k-')
	title('Displacements in Physical Coordinate System, ode45')
	xlabel('Time, sec')
	ylabel('Displacements')
	legend('z1','z2','z3',3)
	grid

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

%	load previous closed-form solutions for tplot, z1, z2, z3 if zeta = 0

	if  zeta == 0

	load tdof_modal_time_z1z2z3;

	plot(t,z_ode(1,:),'k-',t,z_ode(2,:),'k-',t,z_ode(3,:),'k-',tplot,z1,'k.-',tplot,z2,'k.-',tplot,z3,'k.-')
	title('Displacements in Physical Coordinate System from ode45 (ode) and Closed Form (cf)')
	xlabel('Time, sec')
	ylabel('Vibration Displacements')
	legend('ode dof 1','ode dof 2','ode dof 3','cf dof 1','cf dof 2','cf dof 3')
	grid

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

	else
	end
	
%	plot the modal contributions to the motion of masses 1, 2 and 3

	plot(t1,z_ode1(1,:),'k+-',t2,z_ode2(1,:),'kx-',t3,z_ode3(1,:),'k-')
	title('Displacement of dof 1 for Modes 1, 2 and 3, ode45')
	xlabel('Time, sec')
	ylabel('Displacements')
	legend('Mode 1','Mode 2','Mode 3')
	grid

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

	plot(t1,z_ode1(2,:),'k+-',t2,z_ode2(2,:),'kx-',t3,z_ode3(2,:),'k-')
	title('Displacement of dof 2 for Modes 1, 2 and 3, ode45')
	xlabel('Time, sec')
	ylabel('Displacements')
	legend('Mode 1','Mode 2','Mode 3')
	grid

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

	plot(t1,z_ode1(3,:),'k+-',t2,z_ode2(3,:),'kx-',t3,z_ode3(3,:),'k-')
	title('Displacement of dof 3 for Modes 1, 2 and 3, ode45')
	xlabel('Time, sec')
	ylabel('Displacements')
	legend('Mode 1','Mode 2','Mode 3')
	grid