echo off
%	tdofss.m	state-space transfer function solution of tdof undamped 
%	model using	state-space matrices directly and the bode command.
%	Calculates and plots the four distinct frequency responses for the 
%	tdof model.

	clf;

	legend off;

	subplot(1,1,1);

	clear all;

%	define the values of masses, springs, dampers and Forces

	m1 = 1;
	m2 = 1;
	m3 = 1;

	c1 = input('input value for c1, default 0, ...  ');
	
	if  (isempty(c1))
		c1 = 0;
	else
	end

	c2 = input('input value for c2, default 0, ...  ');
	
	if  (isempty(c2))
		c2 = 0;
	else
	end

	k1 = 1;
	k2 = 1;

	F1 = 1;
	F2 = 1;
	F3 = 1;

%	define the system matrix, a

	a = [  0      1         0			 0		   0	  0
		-k1/m1 -c1/m1     k1/m1        c1/m1       0	  0	
		   0      0         0            1         0      0
		 k1/m2  c1/m2  -(k1+k2)/m2  -(c1+c2)/m2  k2/m2  c2/m2
		   0      0         0            0         0      1
		   0      0       k2/m3        c2/m3    -k2/m3 -c2/m3];

%	define the input matrix, b, a 6x3 matrix

	b = [  0	  0		 0
		 F1/m1	  0		 0
		   0	  0		 0
		   0	F2/m2	 0
		   0	  0		 0
		   0	  0	   F3/m3];

%	define the output matrix, c, the 6x6 identify matrix

	c = eye(6,6);

%	define the direct transmission matrix

	d = 0;

%	solve for the eigenvalues of the system matrix

	[xm,omega] = eig(a);

%	Define a vector of frequencies to use, radians/sec.  The logspace command uses
%	the log10 value as limits, i.e. -1 is 10^-1 = 0.1 rad/sec, and 1 is
%	10^1 = 10 rad/sec.  The 200 defines 200 frequency points.

	w = logspace(-1,1,200);

%	use the "ss" function to define state space system for three inputs, forces at 
%	masses 1, 2 and 3 and for all 6 states, three displacements and three velocities

	sssys = ss(a,b,c,d);

%	use the bode command with left hand magnitude and phase vector arguments
%	to provide values for further analysis/plotting

%	the mag and phs matrices below will be 6x3x200 in size
%	the appropriate magnitude and phase to plot for each transfer function
%	are called by appropriate indexing

%	first index 1-6:  z1 z1dot z2 z2dot z3 z3dot
%	second index 1-3:  F1  F2  F3
%	third index 1-200:  all frequency points, use ":" 

	[mag,phs] = bode(sssys,w);

	z11mag = mag(1,1,:);
	z11phs = phs(1,1,:);

	z21mag = mag(3,1,:);
	z21phs = phs(3,1,:);

	z31mag = mag(5,1,:);
	z31phs = phs(5,1,:);

	z22mag = mag(3,2,:);
	z22phs = phs(3,2,:);

%	calculate the magnitude in decibels, db

	z11magdb = 20*log10(z11mag);

	z21magdb = 20*log10(z21mag);

	z31magdb = 20*log10(z31mag);

	z22magdb = 20*log10(z22mag);

%	plot the four transfer functions separately, in a 2x2 subplot form

	subplot(2,2,1)
	semilogx(w,z11magdb(1,:),'k-')
	title('state space, z11, z33 db magnitude')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,2)
	semilogx(w,z21magdb(1,:),'k-')
	title('state space, z21, z12, z23, z32 db magnitude')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,3)
	semilogx(w,z31magdb(1,:),'k-')
	title('state space, z31, z13 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

	subplot(2,2,4)
	semilogx(w,z22magdb(1,:),'k-')
	title('state space, z22 db magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude, db')
	axis([.1 10 -150 50])
	grid

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

	subplot(2,2,1)
	semilogx(w,z11phs(1,:),'k-')
	title('state space, z11, z33 phase')
	ylabel('phase, deg')
	%axis([.1 10 -400 -150])
	grid

	subplot(2,2,2)
	semilogx(w,z21phs(1,:),'k-')
	title('state space, z21, z12, z23, z32 phase')
	ylabel('phase, deg')
	%axis([.1 10 -400 -150])
	grid

	subplot(2,2,3)
	semilogx(w,z31phs(1,:),'k-')
	title('state space, z31, z13 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	%axis([.1 10 -400 -150])
	grid

	subplot(2,2,4)
	semilogx(w,z22phs(1,:),'k-')
	title('state space, z22 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	%axis([.1 10 -400 -150])
	grid

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