function G=trfcont(th,nnu,w)
%TRFCONT Computes the continuous time transfer function of a model.
%	A recommended newer variant to use is TH2FF(THD2THC(TH))
%	G = trfcont(TH)
%
%	TH: A matrix defining a model, as described in HELP THETA
%
%	G is returned as the transfer function of a continous time counter-
%	part of the model TH. This functions is of the standard frequency
%	function format (see HELP FREQFUNC). Note that no automatic compu-
%	tation of the noise spectrum is offered.
%	
%	If the model TH has several inputs, G will be returned as the transfer
%	functions of selected inputs # j1 j2 .. jk by 
%	G = trfcont(TH,[j1 j2 ... jk])  [default is all inputs]
%	The functions are computed at 128 logathmically spaced frequency-values
%	going about a decade above the Nyquist frequency. The functions can be
%	computed at arbitrary frequencies w (a row vector, generated e.g. by 
%	LOGSPACE) by G = trfcont(TH,ku,w).
%	The transfer function can be plotted by BODEPLOT. bodeplot(trfcont(TH))
%	is a possible construction.
%	See also TRF and TRFSD.

%	L. Ljung 7-7-87
%	Copyright (c) 1987 by the MathWorks, Inc.
%	All Rights Reserved.
 
T=th(1,2); nu=th(1,3);

% *** Set up default values ***
wdef=logspace(log10(pi/T/100),log10(10*pi/T),128);
if nargin<3, w=wdef;,end
if nargin<2 , nnu=1:th(1,3);end
if length(w)==1, if w<0, w=wdef;end,end
if length(nnu)==1, if nnu<0, nnu=1:th(1,3);end,end

%  *** Transform to continuous time (remove delays)***

[bc,fc]=contin(th,nnu,1);

nk=th(1,2*nu+7:3*nu+6);   % The delays

[nrb,ncb]=size(bc);[nrf,ncf]=size(fc);
nm=max(ncb,ncf);
   i=sqrt(-1);


% *** Compute the transfer function(s) GC=bc/fc ***
nrc=length(w)+1;
sc=1;
for k=nnu
    G(1,sc:sc+2)=[100+k,k,20+k];	
    G(2:nrc,sc)=w';
    GC=polyval(bc(k,:),i*w)./polyval(fc(k,:),i*w);
    G(2:nrc,sc+1)=abs(GC');
    G(2:nrc,sc+2)=(phase(GC)'-w'*T*max(0,(nk(k)-1)))*180/pi; % Add the delays
    sc=sc+3;
end

