% Terzaghi 1D Consolidation
% Maggie Parks 11/28/12

clear; close all

% Geometry of soil column
N = 42; % N sections
zmax = 21 ; % Bottom of soil column (m)
imax = 2000; % Number of iterations

zwater = 3.5 ; % Water table  (m)
zsoil1 = 3.5 ; % Bottom of Layer 1 (m)
zsoil2 = 7 ; % Bottom of Layer 2 (m)
zsoil3 = 21 ; % Bottom of Layer 3 (m)

% Soil constants
cv = 9.3 ; %Coefficient of consolidation (m^2/min)
gammawater = 9.81 ; % Density of water (kN/m^2)
gammasoil1 = 18.1 ; % Soil layer 1 (kN/m^2)
gammasoil2 = 19.6 ; % Soil layer 2 (kN/m^2)
gammasoil3 = 14.7 ; % Soil layer 3 (kN/m^2)

% Initial and Boundary Conditions
z = linspace( 0 , zmax , N+1); % depth
stress = linspace( 0 , zmax , N+1); % Stress with depth vector n+1 long




% Calc stress with depth
for i = 1 : N+1 ;
    if z(i) < zsoil1 ;
        stress(i) = z(i).* gammasoil1 ;
    elseif z(i) < zsoil2 ;
        stress(i) = (zsoil1 * gammasoil1) + ((z(i) - zsoil1) * gammasoil2)...
            - ((z(i) - zwater) * gammawater);
    else
        stress(i) = (zsoil1 * gammasoil1) + ((zsoil2 - zsoil1) * gammasoil2)...
            + ((z(i) - zsoil2) * gammasoil3) - ((z(i) - zwater) * gammawater);
    end
end


    
% Calc excess pore pressure
dsigma = 100; % Applied stress increment (kPa)
finalstress = stress+dsigma;
uexcess = linspace( 0 , zmax , N+1); % Excess pore pressure
uexcess(:) = 0 ;
uexcess( N*zsoil2/zmax + 2 : N) =dsigma ;

% Now, calculate the dissipation of pore pressure
unew = uexcess;
plotmaxx = 400;
f=0;
stressnew = finalstress - uexcess;
for i = 1:imax ;
    for j = (N*zsoil2/zmax + 2):N ;
        unew(j)=0.5*(unew(j-1)+unew(j+1));
    end
       
    
    if mod(i+9,10) == 0;
        f=f+1;
        
        plot(stress,z)
        
        set(gca,'YDir','reverse','XAxisLocation','top');
        axis([0 plotmaxx 0 zmax])
        xlabel('Stress (kPa)')
        ylabel('Depth (m)')
        hold all
        
        plot(0:plotmaxx,zsoil1,'-')
        plot(0:plotmaxx,zsoil2,'-')
        plot(stressnew,z,'--r')
        
        plot(finalstress,z)
        
        %saveas(gcf,['Layers',sprintf('%03d',f),'.png']);
        pause(0.01)
    end
    
    uexcess = unew ;
    stressnew = finalstress - uexcess;
    
hold off

end

