% Vertical Stress Caused By A Strip Load (Finite Width And Infinite Length)
% Hello and welcome to my Mining Engineering ES Project/Code.
% The basic function of the code relies on the idea of determining the
% pressure under a sidewalk.  The sidewalk has constant width and weight
% but infinite length.  This code determines different pressures as we plot
% out deeper and further away.  From this project, all we want is a plot.
% Lastly, we recommend the user input the 3 inputs in the following order:
% 10, 1, 1.  :)

% q = load/unit area 
% B = 1 (width of the road, not length)
% Delta_z would be the depth under the ground surface
% i = counter for rows and shows depth or z
% z = depth measured from ground surface
% Delta_x would be the distance away from the sidewalk
% j = equal to x, is a counter for columns
% x = distance from sidewalk
% delta_p = shows the pressure at depth "z" and distance "x"
% from the center of the load

clear all
close all
q=input('Enter the load: '); % This would be the weight of your sidewalk.
B=1; % This is the width of your sidealk; it's a constant.
delta_z=input('Enter a positive integer for delta_z: '); % Delta_z would be the depth under the ground surface
i=11:-1:1; % There are 11 rows.
z=(i-1)*delta_z; % This is the depth under the ground surface ranging from row 1-11.
delta_x=input('Enter a positive integer for delta_x: '); % Delta_x would be the distance away from the sidewalk
j=1:1:11; %There are 11 columns.
x=(j-1)*delta_x; % This is the distance away from the sidewalk ranging from column 1-11.



for i=1:11 % I is from rows 1-11.
    for j=1:11 % J is from columns 1-11.
        delta_pq(i,j)=(q/pi)*(atan((x(j)+(B/2))/z(i))-atan((x(j)-(B/2))/z(i))-((B*z(i)*(x(j)^2-z(i)^2-(B^2/4)))/((x(j)^2+z(i)^2-(B^2/4))^2+B^2*z(i)^2))); % The master code (formula)!!
    end
    aa=size(delta_pq(i,:)); % An addition tp the master code to help graph.
    y=10:-1:0; % An additional command so that it won't graph below 0.
    qq=size(y);
end
    contour(delta_pq); % The required plot.

    
% The following 3 lines of code were supposed to be a mirror image of the 
% above code; just reversed on the other side of the X-axis.
% However, both the head TA and one of the other TAs (along w/ myself) just
% couldn't figure it out.
m2=fliplr(delta_pq);
m=m2+delta_pq;
contour(m2)