Plot

global Q S T t r

Q = input('Enter the Discharge Rate Q (cubic meters per day) = ');


S = input('Enter the Storativity S (dimensionless) = ');


T = input('Enter the Transmissivity T (square meters per day) = ');

r=[.5:2:200];
r1=100;

%At time t = 1 hour
t=1;
draw1=draw(r,Q,S,T,t);
time1=draw(r1,Q,S,T,t)+.5;

%At time t = 1 day
t=24;
draw2=draw(r,Q,S,T,t);
time2=draw(r1,Q,S,T,t)+.5;

%At time t = 30 days
t=720;
draw3=draw(r,Q,S,T,t);
time3=draw(r1,Q,S,T,t)+.5;

%Draw the plot for t = 1 hour; t = 1 day; and t = 30 days
plot(r,-draw1,'y-',r,-draw2,'m--',r,-draw3,'g-.');
title('Cone of Depression');
xlabel('Radial Distance (in meters)');
ylabel('Drawdown');
text(100,-time1,'Time = 1 Hour')
text(100,-time2,'Time = 1 Day')
text(100,-time3,'Time = 30 Days')

function d=draw(r,Q,S,T,t)
%
%  function d = draw(r, Q, S, T, t)
%
%  This is C. V. Theis's drawdown function
%
%  It is a function of the radial distance, r, and the time, t, since the 
%  pumping began. r may be a vector, but t must be a scalar.
%
%  The parameters are:
%     discharge rate, Q, [L^3/T]
%     storativity, S, [dimensionless]
%     transmissivity, T [L^2/T].

u = (S/(4*T*t))*r.*r;
d=(Q/(4*pi*T)).*well(u);


function w=well(u)
%
%  function w= well(u)
%
%  Theis's well function
%

Econst = 0.577215664901532860606512;

n = 10;
term = 1;
for k = n:-1:2  
   term =  1 + ((1-k)/k^2)*u.*term;
end
w = -Econst - log(u) + u.*term;

• Example of a Plot