function [magout,phout] = ggen(model,magtol,phtol)
%%%FIX function [mag,phase] = ggen(model,a,w,ns,magtol,phtol)
%GGEN   Generate the SIDF input-output model by direct simulation
%       ggen('model',a,w,ns,magtol,phtol) returns the magnitude and phase
%       of the SIDF input-output model G(jw,a) corresponding to the
%       nonlinear system (set of first-order ODEs) in the file model.m
%       created for simulation via the matlab numerical integrator ode45.
%
%	CALL:
%       [mg,ph] = ggen('model', a, w, ns, magtol, phtol)
%
%	INPUT:
%	model   - String containing name of user-supplied nonlinear
%                   system model suitable for simulation via ODE45
%                   that includes a sinusoidal input and two extra
%                   additional states for generating Fourier integrals
%       a       - amplitude of the sinusoidal input (must be global)
%       w       - frequency of the sinusoidal input (must be global)
%       ns      - number of states in the ODE model, including the
%                   Fourier integrals
%       magtol  - tolerance on magnitude of G(jw,a) (convergence test)
%       phtol   - tolerance on magnitude of G(jw,a) (convergence test)
%
%	OUTPUT:
%	magout  - Returned SIDF I/O model magnitude, |G(jw,a)|, in dB
%	phout   - Returned SIDF I/O model phase, < G(jw,a), in degrees
%
%%%%%%% JH Taylor - University of New Brunswick - 6 July 2002

%Initialize
load temp_gwsim % obtain a, w, ns provided by sidf_io
dpr = 180/pi;
k = 0; T = 2*pi/w;
t0 = 0; tf = T; tspan = [ t0 tf ];
x0 = zeros(ns,1);
[t,x] = ode45(model,tspan,x0);
[nrows,ncols] = size(x);
xf = x(nrows,:);
magast = w/(pi*a)*abs(xf(ncols-1)+j*xf(ncols));
phaseast = atan2(xf(ncols),xf(ncols-1));
while (k >= 0)
	k = k+1;
	if rem(k,10)==0 
	  disp(['Info: # cycles = ',num2str(k),' - probably high frequency']);
	end;
	t0 = k*T; tf = t0+T;
	x0 = xf;
	x0(ncols-1) = 0;
	x0(ncols) = 0;
	[t,x] = ode45(model,tspan,x0);
	nrows = size(x,1);
	xf = x(nrows,:);
	mag = w/(pi*a)*abs(xf(ncols-1)+j*xf(ncols));
	phase = atan2(xf(ncols),xf(ncols-1));
	magdiff = abs(20*log10(mag/magast));
	phasediff = (180/pi)*abs(phase-phaseast);
	if ((magdiff >= magtol) | (phasediff >= phtol))
		magast = mag;
		phaseast = phase;
	else
		k = -1;
	end
end;
magout = mag;
phout = phase .* dpr;  %% change radians to  degrees
