## Copyright (C) 2006 Doug Stewart
## 
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

## ret=ilt(sys [,Rtype])
##  This function calculates the inverse Laplace transform of simple laplace space signals.
##  sys is a system type function ( struct )
##  Rtype is the respnse type that you are requesting
##       if Rtype =
##                  ' s '  = step response
##                  ' i '  = impulse response
##              The default is step response.
##   ret is a string containing the time space info  F(t)  
## If you want the result when the input is a sin or cos function then
## apply this before you make the sys and then ask for in impulse responce.
##
##  Example
##   numerator = 3                        3
##   denominator =[1 4 3]         (s+1)(s+3)
##   transferFunction =tf(numerator,denominator)
##  The equasion of the step responce in time space is   ilt(transferFunction,'s')
##                        the result is                       F(t)  =  + 1 -1.5 e^(-t)+ 0.5 e^(-3t)
##   and the equasion of the impulse responce in time space is ilt(transferFunction,'i')
##                         the result is                       F(t)  =  + 1.5 e^(-t) -1.5 e^(-3t)
##   and the inverse Laplace transform of the transferFunction is  ilt(transferFunction,'i')
##                         the result is                       F(t)  =  + 1.5 e^(-t) -1.5 e^(-3t)
##  If you just want the inverse Laplace transform of some signal then put the signal into a sys type 
##  variable and do    ilt(sys,'i')
##
## The method used is to calculate the partial fraction expansion using the residue funtion
## and then form up the Time space answer.
## This has been tested with real roots non repeating , real roots repeating, complex roots and
## imaginary roots.
## This function uses the residue function and thus the acuracy of the values  is 
##  what ever "residue" returns.

## Author: Doug Stewart
## 
## 2006-12-11 
## * Initial revision
## rev 2  Nov 2007
## added impulse and fixed some other problems.

 function [ ret ] = ilt (sys,Rtype)

if(nargin<1 | nargin>2)
 usage("[ ret ] = ilt (sys,Rtype)     Type    help ilt     for more info.")
endif 

[n d]=sys2tf(sys);
 
 if(nargin==2)
      if(Rtype=='i' | Rtype=='I')
      constantterm=false;
      else
      d=conv([1 0],d);
      endif
 else
     d=conv([1 0],d);
endif


ss=cell(30,1);
q=1;
if(length(n)==length(d))
    [b n]=deconv(n,d);
      ss(q,1)=cellstr(sprintf("impulse  "));          
     q=q+1;
endif

[r p k e1]=residue(n,d);

# we must sort then so that the complex conjugate pairs are together.

 [e1 kk]=sortcp(e1);
 p=p(kk);
 r=r(kk);
%[p kk]=sortcp(p)
%r=r(kk)


# try and get rid of roundoff errors

for k =1 :length(r)
if(!isreal(p(k)))
if(abs(real(p(k))) > 1000*abs(imag(p(k))));
p(k)=real(p(k));
    
 endif
 endif
endfor

k=length(r);

while ( k>0)                

## the poles are real
if (isreal(p(k)))
   if (r(k)>=0)
        ss(q,1)= ' +';  ## insert a "+" for positive and nothing for negative
	q=q+1;
   endif
   if(abs(p(k))<100*eps )  %% pole is at the origin
                ss(q,1)=cellstr(sprintf(" %g ",r(k)));  
               
  else
                 ss(q,1)=cellstr(sprintf(" %g *",r(k)));  
	
   endif
  
   q=q+1;
  
## handle the repeating  real poles  
   if(e1(k)>1)
        if(e1(k)==2)
	   ss(q,1)=cellstr(sprintf(" t  .*" ))   ;  ## if only one t
	   q=q+1; 
	else
             ss(q,1)=cellstr(sprintf(" t.^%g  * 1/ %g   .*",e1(k)-1,factorial(e1(k)-1)))   ;   ## more than 1 t
           q=q+1; 
	endif   
   endif  
   if(p(k)!=0)
        if(abs(abs(p(k))-1)<10*eps)    % if exactly 1
	    if(p(k)>0)
	        ss(q,1)=cellstr(sprintf(" e.^(t)  "  ));
	    else
	       ss(q,1)=cellstr(sprintf(" e.^(-t)  "  ));
	    endif
	else
	     ss(q,1)=cellstr(sprintf(" e.^(%g*t)  " ,p(k)) );% if not 1
        endif
      q=q+1;	
   endif
   k=k-1;
   
endif
if(k>1)   ## must have so we don't count out of range                                                                   
  ## now handle the conplex poles
  
    
  if (!isreal(p(k)))
    a=(real(r(k)));
    b=(imag(r(k)));
    c=(real(p(k)));
    d=abs(imag(p(k)));
     
         
    if(abs(a)<=10*eps)  ## if a=0  therefore  pure imaginary residue
   
        if (b>=0)
          ss(q,1)=' +';
	  q=q+1;
        endif
        ss(q,1)=cellstr(sprintf(" %g  *",2*b) );
        q=q+1;
  
        if(e1(k)>1)
           if(e1(k)==2)
	      ss(q,1)=cellstr(sprintf(" t .*" ))   ;  ## if only one t
	      q=q+1; 
	   else
		ss(q,1)=cellstr(sprintf(" t.^%g  * 1/ %g   .*",e1(k)-1,factorial(e1(k)-1)))   ;   ## more than 1 t
          
              q=q+1; 
	   endif   
        endif    
    
     
	 if(abs(abs(c)-1)<10*eps)  ## if c=1
	        if(c>=0)             ## is c is positive 
		      ss(q,1)=cellstr(sprintf("  e.^(t)   .*sin(%g*t)",d) );
		 else                  ## if c is negative put a  "-"
		      ss(q,1)=cellstr(sprintf("  e.^(-t)   .*sin(%g*t)",d)) ;
		endif
	 else         ## if c != 1
		 
		if (c==0)
			ss(q,1)=cellstr(sprintf("  sin(%g*t)",d)) ;
		else
			ss(q,1)=cellstr(sprintf("  e.^(%g *t)   .*sin(%g*t)",c,d)) ;
		endif
         endif
   else
        if (a>=0)
           ss(q,1)=' +';
	   q=q+1;
        endif 
    
       ss(q,1)=cellstr(sprintf(" %g *",2*a) );
       q=q+1;
  
       if(e1(k)>1)
          if(e1(k)==2)
	     ss(q,1)=cellstr(sprintf(" t  .*" ))   ;  ## if only one t
	     q=q+1; 
	   else
              ss(q,1)=cellstr(sprintf(" t.^%g  * 1/ %g   .*",e1(k)-1,factorial(e1(k)-1)))   ;   ## more than 1 t
              q=q+1; 
	   endif   
        endif  
    
    
    
 
         if(abs(c)<=10*eps)  ##  if c==0
	
            if(b/a<0)
                 ss(q,1)=cellstr(sprintf("  ( cos(%g*t) + %g *sin(%g*t))",d,(b/a),d)) ;
            else
                if(b/a==0)
		ss(q,1)=cellstr(sprintf("   cos(%g*t) ",d)) ; 
	     
		else
		ss(q,1)=cellstr(sprintf("  ( cos(%g*t) + %g *sin(%g*t))",d,(b/a),d)) ; 
	        endif
            endif
         else
            if(abs(abs(c)-1)<10*eps)  ## if c=1
	        if(c>=0) 
	             if(b/a<0) 
                        ss(q,1)=cellstr(sprintf(" e.^(t) .*( cos(%g*t)    %g *sin(%g*t))",d,(b/a),d)) ;
                     else
                        ss(q,1)=cellstr(sprintf(" e.^(t) .*( cos(%g*t) + %g *sin(%g*t))",d,(b/a),d)) ;
                     endif
	        else
	             if(b/a<0)
	                ss(q,1)=cellstr(sprintf(" e.^(-t) .*( cos(%g*t)    %g *sin(%g*t))",d,(b/a),d)) ;
	             else
	                ss(q,1)=cellstr(sprintf(" e.^(-t) .*( cos(%g*t) + %g *sin(%g*t))",d,(b/a),d)) ;
	             endif
	        endif
            else
                if(b/a<0)
                  ss(q,1)=cellstr(sprintf(" e.^(%g *t) .* ( cos(%g*t)   %g *sin(%g*t))",c,d,(b/a),d) );      
                else
                  ss(q,1)=cellstr(sprintf(" e.^(%g *t) .*( cos(%g*t) + %g *sin(%g*t))",c,d,(b/a),d) );
	        endif
	     endif      
         endif
     endif   
     q=q+1;
     k=k-2;
  endif
endif
endwhile 


##ss
%%w='f(t)  =  ';
w='yy= ';
for h=1:q-1
  w=strcat(w,strvcat( ss(h,1)));
endfor
w=strcat(w,";");
ret=w;

endfunction