function [Sa,drift,runs,NCSa,transout]=traceIDA(maxrpr,iniSa,incSa,...
            Th_drift,Th_Sa,mode,hf,eSa,edr,blackbox,varargin)
% 
% IDA curve tracing routine. 
% This is exactly the same as traceDPO4.m with the only difference
% that it also passes the number of runs executed to the analyzing routine
% that it calls for every dynamic analysis
%-------------------------------------------------------------------------
% INPUT
% maxrpr  : Maximum runs per record. This sets how many runs the routine
%           is allowed to use to trace the DPO
% iniSa   : This is the initial SA  step to be used.
% incSa   : This is the step increment for each iteration. If it is 0
%           then we essentially keep the same constant-length step
%           during Hunt-up
% Th_drift: The maximum drift (usually 10%) up to which you want to 
%          trace the curve.
% Th_Sa   : The maximum Sa value allowed to be run (guards against
%           excessively unrealistic runs being performed for 
%           non-collapsing structures
% mode    : mode(1) => Tracing mode, 
%           0=economy(use as few runs as possible), 
%           1=semi (use a bit more)
%           2=detailed (use all)
%           mode(2) => Stopping-rule mode
%           0=maxrpr is strict :
%             if you expend all "maxrpr" runs and still haven't
%             satisfied any (other) stopping rule, then STOP.
%           1=maxrpr is not strict :
%             the HUNT-UP routine is allowed to use as many runs as
%             needed to meet at least one other stopping
%             rule. Still NSECT and the rest phases will not be performed
%             CAUTION: This may lead the program to infinite loops,
%             especially if the structure can't collapse and there
%             are no Th_Sa, Th_drift limits.
% NCmode  : 0=stop when you find first Non-Convergence (NOT IMPLEMENTED)
%           n=stop after n-steps non-converging (guards against possible
% hf      : 3x1 vector of fractions of maxSA (curve height) used as error 
%           control. For hf(1) a usual value is 5..6 (1/5..1/6), for
%           hf(2)=0.5*hf(1) and for hf(3)=hf(1)*2 are typical values.
%           if they are left 0 then the default values hf(1)=5, hf(2)=2.5
%           and hf(3)=10 are supplied automatically.
%           The higher the fraction is the more runs and the better 
%           the accuracy: hf(1) controls accuracy for tail correction,
%           hf(2) controls body filling (useful for semi mode only) and
%           hf(3) controls if Th_drift is exceeded. 
% eSa,edr : If an elastic run has been performed seperately then the user
%           may provide it thru these variables. If eSa is provided but
%           edr is zero, then the suggested run will be performed.
%           If esa is 0 then a run will be performed at eSa=0.01g.
%           This should be elastic for most structures that I can think
%           of!
% blackbox : the name of the function to use for performing each run
% varargin : this catches all additional input required by the blackbox 
%            routine
% OUTPUT 
% Sa     : Sa values obtained by the routine (g)
% drift  : drift values obtained by the routine
% runs   : total runs used up by the routine(<=maxrpr)
% NCSa   : minimum NC point Sa (useful for error control)
% transout : a cell array 1xlength(Sa), same length as Sa,drift, containing
%            extra output supplied by the black-box routine.
%-------------------------------------------------------------------------
% It can perform quite bad, if the DPO is weird. Definately you can't
% pick all possible features if you aren't supposed to use many points!
% Basic features: Even if you provide a very low maxrpr, it will still 
% run a minimum number of runs required to get some approximation of
% the shape. So very low maxrpr's aren't good. I suggest 8-12, depending
% on the expected capacity of the structure (high capacity-->taller curve
% so more points needed.
% HINT : If you use the "detailed" mode (i.e you want to expend all your
% runs), it makes sense to allow for higher hfrac=10..20, as this
% will improve accuracy a lot close to the collapse limit, while a 
% lower hfrac, would still leave less accuracy there and most of the
% remaining runs would be spend lower, filling in the gaps!
% CAUTION : don''t use ppval to evaluate natural splines
% It isn''t able to work with them. USE ppual or fnval! 
%-------------------------------------------------------------------
% Now you are able to prescribe an elastic run or prescribe how to
% perform it. Also all height fractions are now set externally.
%--------------------------------------------------------------------
% Blackbox: The blackboxx analysis routine must have the following 
% structure [drift,NCflag,transout]=blackbox(Sa,...any variables you want)
% see dyn_push3() and PushCurve7() for examples.
%-------------------------------------------------------------------
% NOTE: Whenever transout(i) is an empty cell, this means that
%      this point didn''t directly come from a run, but was
%      computed instead. Examples of these are
%        a) the NCSA point (we can''t compute transout there)
%        b) the (0,0) point (transout is all 0 there)
% (Last Updated: May 14th, 2001)
% Copyright by D.Vamvatsikos


%// if hf(2) is 0 then supply the default value.
if (hf(2)==0), hf(2)=hf(1)*0.5; end
%// This one is useful to control how many runs will be performed
%// if the interpolation is insufficient to calculate a point @Th_drift..
%// if hf(2) is 0 then supply the default value
if (hf(3)==0), hf(3)=hf(1)*2; end

% default is "maxrpr is strict"
if length(mode)==1, mode(2)=0; end


% some erroneous input control!
if (Th_Sa<iniSa)|(Th_Sa<eSa)
  error('Not allowed to have Th_Sa < max(iniSa,eSa)'); 
end;

% initialize
runs=0;
Sa(1)=0;
drift(1)=0;
NC(1)=0;
transout(1)={cell(0)};
i=2;

% run first at very low Sa, guaranteed to produce elastic behavior.
if (eSa==0)  
  %user hasn't supplied any elastic run info. Use defaults
  Sa(i)=0.01;
  % use feval to call the blackbox routine. Note
  % (a) transout is a cell array but we can address cells as transout()
  % (b) vararrgin{:} produces a comma seperated list of all additional
  % arguments
  [drift(i),NC(i),transout(i)]=feval(blackbox,Sa(i),runs,varargin{:});
  runs=runs+1; i=i+1;
else
  if (edr==0)   
    % user has supplied a value of Sa to do an elastic run
    Sa(i)=eSa; 
    [drift(i),NC(i),transout(i)]=feval(blackbox,Sa(i),runs,varargin{:});
    runs=runs+1; i=i+1;
  else  
    % user has supplied the elastic run. No run performed 
    Sa(i)=eSa;
    drift(i)=edr;transout(i)={cell(0)};
    i=i+1;
  end
end



%// HUNT UP. Stop when you bracket the DPO
%// e.g. stop when you hit NC for the first time
%// There is a  maxrpr limit here. It will not use as many runs as needed
%// to get the DPO shape.

% POSSIBLE HUNT-UP IMPROVEMENT: After you exceed Th_drift or reach NC
% do one more step to see if you exceed it again (or get NC again) as a
% verification that you have indeed found THE flatline.


% Initialize
% if incSa==0, then we do constant-stepping (just like most people do)
% NOTE: Our first "hunt-up" run will be eSa+iniSa. This may cause some
% problems if you supply the elastic run "foutch-style", i.e. eSa is
% not truly an SA with elastic behavior (Foutch performs his elastic
% runs at an arbitrary high SA but makes sure the model is an elastic
% one (simply makes the yield strength huge).. This should be explicitly
% avoided at ALL times!! The same model is to be used for all runs as
% the postprocessing is not able to anticipate that and havoc will follow!

dSa=iniSa;
Sa(i)=Sa(i-1)+dSa;

HUTflag=0;
while ~HUTflag
  [drift(i),NC(i),transout(i)]=feval(blackbox,Sa(i),runs,varargin{:});
  runs=runs+1;i=i+1;
  % I prefer a condition like "isinf(drift(i))" but some routines
  % return a NaN, so this may fail.
  
  % HUTflag is the Hunt-Up Termination flag, and reports which rule has been
  % invoked to stop the hunt-up.
  if (NC(i-1)==1), HUTflag=1; 
  elseif (drift(i-1)>Th_drift), HUTflag=2; 
  elseif (Sa(i-1)==Th_Sa), HUTflag=3; 
  elseif (runs==maxrpr)&(mode(2)~=1), HUTflag=4;    
  else  
    dSa=dSa+incSa;
    Sa(i)=Sa(i-1)+dSa;
    % if you exceed Th_Sa, set new run at Th_Sa exactly. Next loop
    % will actually exit after this run.
    if Sa(i)>Th_Sa, Sa(i)=Th_Sa; end;
  end  
end

%keyboard;

% N-Secting. This converges closer to the collapse capacity point.
% If it is non-convergence that rules, then carefull tri-secting is
% performed to get closer to the flatline, otherwise bisecting is done
% to get closer to Th_drift
if (HUTflag==1)|(HUTflag==2)
  % Now try N-SECTING to locate better point. Seperate the two points
  %that bracket the "capacity" into vectors of Sa, drift
  diSa=[Sa(i-2),Sa(i-1)];
  didr=[drift(i-2),drift(i-1)];
  difSa=diSa(2)-diSa(1);
  % so HUTflag==1 produces Nsect=3, HUTflag==2 produces Nsect=2
  % hence we get trisecting and bisecting respectively.
  Nsect = 4-HUTflag ;
  while 1 
    if ((Nsect==3)&(difSa<diSa(1)/hf(1))) |...
       ((Nsect==2)&(difSa<diSa(1)/hf(3))) |(runs==maxrpr)
       break
    end   
    
    Sa(i)=diSa(1)+(diSa(2)-diSa(1))/Nsect;
    
    [drift(i),NC(i),transout(i)]=feval(blackbox,Sa(i),runs,varargin{:});
    if (NC(i)==1)|(drift(i)>Th_drift)    
      % we found a lower "NC" point
      if NC(i)==1, Nsect=3; else Nsect=2; end
      diSa=[diSa(1),Sa(i)];
      didr=[didr(1),drift(i)]; 
    else
      % a higher "C" (converging)  point was found
      diSa=[Sa(i),diSa(2)];
      didr=[drift(i),didr(2)];
    end
    runs=runs+1;  i=i+1;
    difSa=diSa(2)-diSa(1);
  end
  % reorder, as this n-secting tends to introduce a wealth of unsorted
  % runs! 
  [Sa,drift,NC,transout]=reorder(Sa,drift,NC,transout);
end  
  

% I used to have a loop here to force a better fit around the Th_drift
% but I think this is a waste of runs. The simple "gap-filling"
% procedure described later serves the same purpose. It used to be that
% Th_drift was a vital piece of info for capacity, but I believe this
% to be less useful now. The only reason one may want to keep it is for
% purposes of using the Th_drift as the absolute measure of capacity
% (due to some fracturing maybe). A good bracketing of Th_drift really
% improves our chances of finding good capacity estimates.
% Even better, I have th_drift as a collapse-indicator now, so it is
% included in the N-secting routine above!


% fill-in the gaps  (semi and detailed modes)
% If an intermediate NC point is found, do not try to N-sect around it.
% Just keep it and let the post-processor worry about this. The
% gap-filling quality criterion is good enough to look around if the
% gaps are wide enough. Otherwise why densify our grid there instead of
% checking other wider gaps for possible collapses?

% using NC to get the converging runs is better that isinf or
% isnan, as some routines use NaN and others use Inf to code the
% drift at non-convergence.

uppermost_Sa=min(Sa(find(NC)));
% if no NC points use the run that exceeded Th_drift (or Th_Sa) as
% the uppermost. Note that we do not try to reduce the gap between
% the uppermost_Sa and the Converging run that is closest to
% it. This is taken care of in the N-secting part. Of course if you
% specify maxrpr=100 and hf(1)=2, there is a good chance that you
% get very good resolution everywhere but close to capacity. But
% that is your problem, you geve in bad parameters. I am not going
% to try and modify this to be more lenient to bad input!
% Gap-filling is strictly within the converging runs (not NC or
% those exceeding Th_drift,Th_sa).
if isempty(uppermost_Sa), uppermost_Sa=Sa(end); end
% this process is based upon locating
% where the max difference lies between CONVERGING runs.
% Actually at start of the gap-filling, all NC points are at the
% upper Sa's, so this is fine. During the gap-filling, an 
% intermediate NC may be detected, so if we use something like
% [maxd,j]=max(diff(Sa(find(~NC))));
% then we will end up performing the same run (that produced the
% intermediate NC) over and over. To counter that, just exclude the
% highest NC's from the gap-filling by use of the "uppermost_Sa".
% this way you will be able to keep up your gap-filling around any
% intermediate NC!

if (mode(1)~=0)
  while runs<maxrpr
    [maxd,j]=max(diff(Sa(find(Sa<uppermost_Sa))));
    if (mode(1)==1)&(maxd<uppermost_Sa/hf(2)) 
      % semi-mode "quality-condition" has been achieved
      break; 
    end   
    % put a point midways
    Sa(i)=0.5*(Sa(j)+Sa(j+1));
    [drift(i),NC(i),transout(i)]=feval(blackbox,Sa(i),runs,varargin{:});
    runs=runs+1;
    % reorder SA,drift,transout by inserting the new run at j-th place.
    indx=[1:j,i,(j+1):(i-1)];
    i=i+1;
    [Sa,drift,NC,transout]=reorder(Sa,drift,NC,transout,indx);
  end
end  


NCSa=min(Sa(find(NC)));
% if no drift is NaN (i.e., NC occured) then NCSa=[], it becomes empty
% instead of NaN as we used to in older versions of traceDPOX.m
% Notice that this forces the use of cell array NCSa, instead of
% simple array NCSa. Instead of fixing everything everywhere else,
% jsut change it to NaN here too! No reason why not to!
if isempty(NCSa), NCSa=NaN; end

% now make them column vectors
Sa=Sa';
drift=drift';

%keyboard


function [Sa,drift,NC,transout,indx]=reorder(Sa,drift,NC,transout,indx);
%___________________________________________________________________
% Simply reorders the three given arrays and cell arrays according to
% increasing Sa. If the "indx" is given, then it is used, otherwise the
% sorting is directly done instead.
%___________________________________________________________________

if nargin<5
  [Sa,indx]=sort(Sa);
else
  Sa=Sa(indx);
end
drift=drift(indx);
NC=NC(indx);
transout=transout(indx);