% Script for finding sheetpile depth with LHS/RHS

% Approaches:
% 1) manual
%    - user enters an interval
%    - user picks tolerance
%    - search inside that interval for a root that "works"
%    - if doesn't work, ask for new interval and tolerance
% 2) smarter
%    - user enters interval
%    - program tries increasingly tighter toleranances
%    - good if interval known to contain root
% 3) even better
%    - user enters a known starting point and "direction"
%    - program marches "forward" until sign changes
%    - program remembers interval surrounding sign change
%    - program tries tighter tolerances until finding solution
%------------------------------------------------------------------------------
% I picked #3
% I think I could be handling the boundaries a bit better.
% Also, if the root "hits" zero that's bad! (can't have negative sheet pile)
% Maybe the program could alter the initial interval?
%------------------------------------------------------------------------------
% Refresh the workspace
  clear
  disp(['Starting LHS/RHS...']);
%------------------------------------------------------------------------------
% Given data:
  Gamma = 18;     % kN/m^3, unit weight of soil
  phi   = pi/6;   % 30deg,  measure of soil friction
  L     = 3;      % m,      height of pile above ground cut
  P     = 30;     % kN/m,   force that soil exerts on pile, per unit m
%------------------------------------------------------------------------------
% Guess range for piledepth d (m)
%------------------------------------------------------------------------------
% Get initial guess and direction

% User must enter a positive number
  d = input('Enter your initial guess for pile depth [positive]: ');
  while( ~isnumeric(d) | d < 0 )
     d = input('No! You must enter a positive number: ');
  end

% Get initial direction
  neg=-1;pos=1;
  dir = input('What direction do wish to iterate? [pos,neg] ');
  while( ~(dir==1 | dir==-1) )
     dir = input('No! You must choose either pos or neg: ');
  end

% Should also get boundaries or outer limits.
% i.e., is there some outer interval not to exceed?
%------------------------------------------------------------------------------
% Initialize more variables
  eqn = piledepth(phi,Gamma,P,L,d);  % get initial value of equation
  oldsign = sign(eqn);               % store sign of equation
  target = 0;                        % stopping value
  iterations = 0;                    % count of iterations to reach root
  inc = 0.001;                       % value to increment $d$ in each step
  eps = 0.001;                       % tolerance for stopping iteration
  refine = 10;                       % amount to divide $inc$ if iteration
                                     %   is too "coarse"
%------------------------------------------------------------------------------
% Iterate starting at $d$ and use $dir$ for direction
% Want root of equation, so check $eqn-target$
% Stop iterating when "close" (near eps)
 
  while ( abs(eqn-target) > eps )

     % Need to guess next value
     % increment $d$ by $inc$ (could do this more granular?)
 
     d = d + dir*inc;
     eqn = piledepth(phi,Gamma,P,L,d);
     newsign = sign(eqn);
 
     % If sign changes, root was missed
     if newsign ~= oldsign
	dir = -1*dir;      % change direciton
        inc = inc/refine;  % refine increment
        oldsign = newsign; % store the new sign of eqn
     end

     % Need boundary condition -- what happens if hit zero?
     % (user might have gone in the wrong direction!)
     % Need to "swap" (change) direction of iteration
     if d <= eps
        disp(['Holy crap! We went too far...hitting reverse...']);
        d = eps; % don't want zero, because that's a solution!
	eqn = piledepth(phi,Gamma,P,L,d);
        oldsign = sign(eqn);
        newsign = oldsign;
	dir = pos;
        inc = 0.1; % pick "large" inc to "get back" quicker
     end

     % Below is how I handle going "past" a boundary or taking "too long":
     % If iterations > high value, should break from exhaustion....
     % Really should use an addition condition in the $while$ to avoid $break$
     iterations = iterations + 1;
     if iterations > 1000
        break
     end
  end
%------------------------------------------------------------------------------
% Either repeat analysis or report output
  if iterations > 1000
     disp('This is taking too long...try a different interval!');
     clear
     lhsrhsmethod
  else
     disp('Finished LHS/RHS method!')
     iterations
     d
     eqn=piledepth(phi,Gamma,P,L,d)
 end
%------------------------------------------------------------------------------
% Let the user repeat this method, move on, or quit
  restart('lhsrhsmethod');
%------------------------------------------------------------------------------
