Index: mttroot/mtt/bin/trans/m/dm2sr.m
==================================================================
--- mttroot/mtt/bin/trans/m/dm2sr.m
+++ mttroot/mtt/bin/trans/m/dm2sr.m
@@ -1,6 +1,6 @@
-function sr = dm2sr(A,B,C,D,E,T);
+function sr = dm2sr(A,B,C,D,E,T,u0,x0);
 % sr = dm2sr(A,B,C,D,E,T);
 % Descriptor matrix to impulse response.
 % NB At the moment - this assumes that E is unity .....
 % A,B,C,D,E - descriptor matrices
 % T vector of time points
@@ -8,34 +8,45 @@
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % %% Version control history
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % %% $Id$
 % %% $Log$
+% %% Revision 1.3  1996/08/11 19:33:24  peter
+% %% Replaced exp by expm - whoops!
+% %%
 % %% Revision 1.2  1996/08/11 10:37:40  peter
 % %% Corrected mistake in step-response calculation.
 % %%
 % %% Revision 1.1  1996/08/11 09:42:40  peter
 % %% Initial revision
 % %%
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+[Ny,Nu] = size(D);
+[Ny,Nx] = size(C);
+
+if nargin<7
+  u0 = zeros(Nu,1);
+  u0(1) = 1;
+end;
+
+if nargin<8
+  x0 = zeros(Nx,1);
+end;
 
 [N,M] = size(T);
 if M>N
   T = T';
   N = M;
 end;
 
-[Ny,Nu] = size(D);
-NN=Ny*Nu;
+one = eye(Nx);
 
-[N_y,N_x] = size(C);
-one = eye(N_x);
-
-sr = zeros(N,NN);
+sr = zeros(N,Ny);
 i = 0;
 for t = T'
   i=i+1;
-  SR = C*( A\(expm(A*t)-one) )*B + D;
-  sr(i,:) = reshape(SR, 1,NN);
+  expAt = expm(A*t);
+  SR = C*( ( A\(expAt-one) )*B*u0 + expAt*x0) + D*u0;
+  sr(i,:) =SR';
 end;