function [mag,phase] = sidf_io(sim_model,nstates,avec,wvec,magtol,phtol,progress)
%SIDF_IO Generate the SIDF input/output (I/O) model by direct simulation
%        sidf_io('model',ns,avec,wvec,magtol,phtol) returns the magnitude
%        and phase of the SIDF input/output model G(jw,a) corresponding
%        to the nonlinear system described in the file model.m created
%        for simulation via the matlab numerical integrator ode45.
%
%        CALL:
%        [mg,ph] = sidf_io('sim_model', ns, a, w, magtol, phtol)
%
%        INPUT:
%        sim_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
%                    NOTE: this model must read input amplitude and
%                          frequency from the file temp_gwsim
%        nstates   - number of states in the ODE model, including the
%                    Fourier integrals
%        avec      - amplitude(s) of the sinusoidal input
%        wvec      - frequency vector for the sinusoidal input
%        magtol    - tolerance on magnitude of G(jw,a) (convergence test)
%                      optional - by default, magtol = 2 dB
%        phtol     - tolerance on phase of G(jw,a) (convergence test)
%                      optional - by default, phtol = 5 degrees
%        progress  - if ~= 0 the amplitude and frequency are echoed
%
%        OUTPUT:
%        mag       - Returned SIDF I/O model magnitude, |G(jw,a)|, in dB
%        phase     - Returned SIDF I/O model phase, < G(jw,a), in degrees
%
%
%%%%%%% JH Taylor - University of New Brunswick - 9 July 2002

%% This is a "shell" for ggen, which evaluates single instances of G(jw,a)
%%

if nargin > 7, error('too many inputs, try help sidf_io'); end;
if nargin < 7, progress = 0; end; % you don't want the amplitude and frequency
if nargin < 6, phtol = 5; end; % phase tolerance (deg)
if nargin < 5, magtol = 2; end; % magnitude tolerance (dB)
if nargin < 4, error('too few inputs, try help sidf_io'); end;
ns = nstates;
mdl = sim_model;
%
% amplitude loop
for jj=1:length(avec)
   a = avec(jj);
   if progress, disp(['amplitude = ',num2str(a)]); end;
   % frequency loop
   for ii=1:length(wvec)
      w = wvec(ii);
      if progress, disp(['  frequency = ',num2str(w)]); end;
      save temp_gwsim a w ns  %% pass run definition to ggen and the model
      [mag(ii,jj),phase(ii,jj)] = ggen(mdl,magtol,phtol);
   end % frequency loop
end % amplitude loop
