% Solves the Theis model for a single layer confined aquifer with
% radial symmetry.
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
% Copyright Raymond G. Jacquot amd Thomas V. Edgar
clear all

% working in meters and days here are the parameters
K=8.7;b=5.30;
T=K*b;S=2e-5;
r=linspace(1,100,41);
framedata=zeros(400,41);
t=linspace(.005,2,400);
Sover4T=S/(4*T);
u=zeros(800,41);
a=[-.57721566 .99999193 -.24991055 .05519968 -.00976004 .00107857];
for j=1:400% time loop
    for i=1:41 % space loop
        u(j,i)=(r(1,i)*r(1,i)*Sover4T)/(t(1,j));
        u1=u(j,i);
        E1=a(1,1)+a(1,2)*u1+a(1,3)*(u1^2)+a(1,4)*(u1^3)+a(1,5)*(u1^4)+a(1,6)*(u1^5)-log(u1);
        I=(1/(4*pi))*(E1);
        framedata(j,i)=I;
    end
end
figure(1);clf;%polt s vs. r for various values of t
h1=axes('Ydir','reverse');
hold on
index=[1 121 241 400 ];
for j=1:40:361    
    plot(r,framedata(j,:));
    hold on;
end
box on
xlabel('Radial Distance from Well, r [m]')
ylabel('Dimensionless Drawdown, s(r,t)T/Q')
text(60,1.2,['b = ' num2str(b) ' m']);
text(60,1.1,['K = ' num2str(K) ' m/day']);
text(60,1.,['S = ' num2str(S)]);
t1=5;
text(10,1.3,'Press Enter to Continue')
pause
hold off

figure(2);clf;%plot s vs. t for various r
h2=axes('Ydir','reverse');
axis([0 2 0 1.4]);
hold on
index2=[1 3 5 10 15];
for i=1:4:41
    plot(t,framedata(:,i))
    hold on
end
box on
text(1.4,.3,['b = ' num2str(b) ' m']);
text(1.4,.2,['K = ' num2str(K) ' m/day']);
text(1.4,.1,['S = ' num2str(S)]);
xlabel('Time, t [days]')
ylabel('Dimensionless Drawdown, s(r,t)T/Q')
text(.5,1.3,'Press Enter to Continue')
pause
t=[0 t];
framedata=[zeros(1,41);framedata];
figure(3);clf;% plot s vs r and t
axis([0 100 0 4 0 0.6]);
[R,T1]=meshgrid(r,t);
mesh(R,T1,framedata);
view(45,30)
xlabel('Distance from Well, r [m]','Rotation',-20)
ylabel('Time, t [days]', 'Rotation',20)
zlabel('Dimensionless Drawdown, s(r,t)T/Q')
text(30,4,1,'Press Enter to Continue')
pause
lim=max(max(framedata));

figure(4);clf;%plot animation of well drawdown as time evolves
h3=axes('Ydir','reverse');
r1=[-fliplr(r) r];
framedata2=[fliplr(framedata) framedata];
axis([r1(1,1) r1(1,82) -1.1*lim 2*lim])
hold on
patch([r1(1,1) r1(1,82) r1(1,82) r1(1,1) r1(1,1)],[-.8*lim -.8*lim 0 0 -.8*lim],'m');
hold on
xpatch=[r1 r1(1,82) r1(1,1) r1(1,1)];
ypatch=[framedata2(1,:) 2*lim 2*lim framedata2(1,1)];
P1=patch(xpatch,ypatch,'b','EraseMode','Background');
hold on
plot([0 0],[1.9*lim -lim],'k','LineWidth',5);%well pipe
plot([-2 20],[-lim -lim],'k','LineWidth',5)%horiz. well pipe
plot([20 25 30],[-lim -.95*lim -.85*lim],'b','Linewidth',5)% water spout
hold on
text(-80,-.4*lim,'Impermeable layer')
text(-80,1.5*lim,'Aquifer')
T1=text(-80,-.7*lim,'Press Enter to Animate');
pause
set(T1,'String',' ')
for i=1:401
    ypatch=[framedata2(i,:) 2*lim 2*lim framedata2(i,1)];
    set(P1,'Ydata',ypatch);
    plot([0 0],[1.9*lim -lim],'k','LineWidth',5);
    hold on
    pause(.02)
end
text(-80,.3*lim,['t = ' num2str(t(1,401)) 'days']); 
xlabel('Radial Distance from Well, r [m]')
ylabel('Dimensionless Drawdown, s(r,t)T/Q')
text(40,1.5*lim,['b = ' num2str(b) ' m']);
text(40,1.7*lim,['K = ' num2str(K) ' m/day']);
text(40,1.9*lim,['S = ' num2str(S)]);
set(T1,'String','Press Enter to Continue');
set(gca,'Box','on')
pause

figure(5);clf;
for j=1:40:401 
    semilogx(r,framedata(j,:));
    hold on;
end
xlabel('Radial Distance from Well, r [m]')
ylabel('Dimensionless Drawdown, s(r,t)T/Q')
text(2,.2,'Press Enter to Continue')
pause