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_SaTh_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)&(difSaTh_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