function SS_time = hillslope_hw4 (D,I)
% written by: S. Mudd, 2/6/2006
% an evolving hillslope model
% this takes a channel incision rate, I (in m/yr)
% and a hillslope diffusivity, D (in m^2/yr)
% and returns the time it takes for the hillslope to reach 
% steady state, SS_time

dx = 1;             % spacing of the x locations in meters
n_nodes = 41;       % number of nodes in the model

dt = 10;           % time spacing
thresh = I/100;

time_per_frame = 1000;               % time elapsed between frames (yrs)
framespacing = time_per_frame/dt;    %Conditional for taking frames


%create the vector containing x the x locations
for i = 1:n_nodes
    x(i) = dx*(i-1);
end

% create Z and Znew vectors
% these are the intial conditions
for i = 1:n_nodes
    Z(i) = 30;
    Znew(i) = 30;
end

% initialize E(1), so that we don't close the while loop on the first step
E(1) = 0;
j = 0;

% entering the while loop
while thresh < abs(E(1)-I)
    % in the while loop, we arenet looping through an index j, so we'll put that index 
    % in here
    j = 1+j;            % add 1 to whatever j was in the last timestep   
    t_ime(j) = j*dt;       % update the time elapsed
    
    % get the elevation of the channel
    Z(n_nodes) = Z(n_nodes) - I*dt;
    Znew(n_nodes) = Z(n_nodes);
    
    % get the erosion rate at the channel
    E(n_nodes) = I;
    
    % get the elevations of the intermediate nodes
    for i = 2:n_nodes-1
        Znew(i) = dt*D/(dx*dx)*(Z(i+1) - 2*Z(i) +Z(i-1)) + Z(i);
        E(i) = -(Znew(i)-Z(i))/dt;
    end
            
    % get the last node
    Znew(1) = 2*dt*D/(dx*dx)*(Z(2) - Z(1)) + Z(1);
    E(1) = -(Znew(1)-Z(1))/dt;
    
    % reset Z
    Z = Znew;
    
    % now assign Z_chan, which is the elevation of the hillslope surface above the 
    % channel. Recall that the channel is at i = n_nodes 
    Z_bl = Z-Z(n_nodes);
    
    % plot only at times of spacing framespacing
    temp2 = mod(j,framespacing);
    if temp2 == 0
        subplot (2,1,1)
        bar(x,Z_bl);
        title(['time is ' num2str(t_ime(j))])
        axis([0 dx*(n_nodes-1) 0 0.6*x(n_nodes)*x(n_nodes)*I/D]);
        xlabel('distance from divide, x (m)')
        ylabel('elevation of hillslope surface (m)')
        
        subplot (2,1,2)
        plot(x,E,'Linewidth',2);
        title(['time is ' num2str(t_ime(j))])
        axis([0 dx*(n_nodes-1) 0 I*1.2]);
        xlabel('distance from divide, x (m)')
        ylabel('erosion rate (m)')
        
        pause(0.1);
  end
end

% calculate slope and curvature at steady state
Z_bl(n_nodes) = Z_bl(n_nodes) - I*dt;
for i = 2:n_nodes-1
    x_sc(i-1) = x(i);                           % slope and curvature are misssing the 
                                                % boudnary nodes, so we have a new x 
                                                % vector that is of different size
    slope(i-1) = (Z_bl(i+1)-Z_bl(i-1))/dx;  % calcualte slope
    curv(i-1) = (Z_bl(i+1)-2*Z_bl(i)+Z_bl(i-1))/(dx*dx);  % calcualte curvature
end

%plot the topography, slope and curvature of the hillslope
figure
subplot(3,1,1)
plot(x,Z_bl,'LineWidth',2);
title(['D = ' num2str(D) ' and I = ' num2str(I)])
axis([0 dx*(n_nodes-1) 0 max(Z_bl)]);
%xlabel('distance from divide, x (m)')
ylabel('elevation of hillslope surface (m)')

subplot(3,1,2)
plot(x_sc,slope,'g','LineWidth',2);
a = axis;
axis([0 dx*(n_nodes-1) a(3) 0]);
%xlabel('distance from divide, x (m)')
ylabel('slope (dimensionless)')

subplot(3,1,3)
plot(x_sc,curv,'r','LineWidth',2);
a = axis;
axis([0 dx*(n_nodes-1) -I*1.2/D -I*0.8/D]);
xlabel('distance from divide, x (m)')
ylabel('curvature (1/m)')

% return the time when the model reached steady state
SS_time = t_ime(j);