function h = fm_tstep(flag, convergency, iteration,t)
% FM_TSTEP determine the time step during time domain
%          simulations
%
% H = FM_TSTEP(FLAG,CONV,ITER,T)
%       FLAG: 1 - initialized time step and fix maximum time
%                 step
%             2 - check time step and change it if necessary
%       CONV: 1 - last time step computation converged
%             0 - last time step computation did not converge
%       ITER: number of iterations needed for getting the
%             convergence of the last time step computation
%       T:    actual time
%       H:    new time step
%
%see also FM_INT
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    07-Mar-2004
%Version:   1.1.0
%
%E-mail:    fmilano@thunderbox.uwaterloo.ca
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global DAE Syn Settings Fault

switch flag
 case 1

  % estimate the minimum time step
  if DAE.n == 0
    freq = 1;
  elseif DAE.n == 1,
    As = DAE.Fx - DAE.Fy*(DAE.Jlfv\DAE.Gx);
    freq = abs(As);
  else,
    try
      opts.disp = 0;
      if DAE.n <= 6
        freq = max(abs(eig(full(DAE.Fx))));
      else
        freq = max(abs(eigs(DAE.Fx,nfreq,'LI',opts)));
      end
    catch
      freq = 40;
    end
    if freq > Settings.freq,
      freq = Settings.freq;
    end
  end

  % set the minimum time step
  deltaT = abs(Settings.tf-Settings.t0);
  Tstep = 1/freq;
  Settings.deltatmax = min(5*Tstep,deltaT/100);
  Settings.chunk = 100;
  Settings.deltat = Tstep;
  Settings.deltatmin = min(Tstep/64,Settings.deltatmax/20);
  if Settings.fixt
    if Settings.tstep < 0
      fm_disp('Error: fixed time step is negative or zero',2)
      fm_disp('Automatic time step has been set',1)
      Settings.fixt = 0;
    elseif Settings.tstep < Settings.deltatmin
      fm_disp('Warning: fixed time step is less than estimated minimum time step',2)
      Settings.deltat = Settings.tstep;
    else
      Settings.deltat = Settings.tstep;
    end
  end

 case 2

  % check time step
  switch convergency
   case 1,  % should we change the time step?
    if iteration >= 15,
      Settings.deltat = max(Settings.deltat*0.9,Settings.deltatmin);
    end
    if iteration <= 10,
      Settings.deltat = min(Settings.deltat*1.3,Settings.deltatmax);
    end
    if Settings.fixt,
      Settings.deltat = min(Settings.tstep,Settings.deltat);
    end
   case 0,  % bisecting time step if no convergence
    Settings.deltat = Settings.deltat*0.5;
    if Settings.deltat < Settings.deltatmin;
      Settings.deltat == 0;
    end
  end

  % check fault occurrencies
  if Fault.n
    index_times = find (t == [Fault.con(:,5); Fault.con(:,6)]);
    if ~isempty(index_times)
      Settings.deltat = min(Settings.deltat,0.005);
    end
  end

end

h = Settings.deltat;