% This m-file displays the solution to the two layer aquifer problem 
% for radial symmetry. Data for the animation is stored in the mat-file
% framedata.mat.
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
clear all
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
figure(1):clf;
figure(2):clf;
figure(3):clf;
figure(4);clf;
figure(5);clf;
% working in meters and days here are the parameters
bprime=5;Kprime=1.728;S=5e-5;
K=60.48;b=10;
T=K*b;
B=sqrt(T*bprime/Kprime);
Bsq4=B*B*4;
% looks like the upper limit for t should be 200 days
SoverT=S/T;
r=(B/10)*linspace(1, 40,40);
r=[B/40 r];
roverB=r/B;
framedata=zeros(400,41);
t=linspace(.5,200,400);
load framedata;
figure(1);%polt s vs. r/B for various values of t
h1=axes('Ydir','reverse');
hold on
index=[1 11 41 81 161];
for j=1:5
    plot(roverB(1,1:21),framedata(index(1,j),1:21));
    hold on
end
set(gca,'Box','on');
text(.5,.5,['B = ' num2str(B) ' meters']);
text(.5,.45,'o o o Denotes Infinite Time Solution');
text(.2,.025,'t = 0.5');
text(.27,.05,'5');
text(.4,.07,'20');
text(.6,.07,'40');
xlabel('Dimensionless Distance from Well, r/B');
ylabel('Dimensionless Drawdown, s(r,t)T/Q');
assym=(1/(2*pi))*besselk(0,r/B);
plot(roverB(1,1:21),assym(1,1:21),'o');
text(.2,.65,'Press Enter to Continue');
pause
figure(2);clf%plot s vs. t for various r
h2=axes('Ydir','reverse');
hold on;
index2=[1 2 3 5 10];
for i=1:5
    plot(t,framedata(:,index2(1,i)));
    hold on
end
set(gca,'Box','on');
xlabel('Time, t [days]')
ylabel('Dimensionless Drawdown, s(r,t)T/Q')
text(120,.65,['B = ' num2str(B) ' meters']);
text(20,.6,'r/B = 0.025')
text(30,.4,'r/B = 0.1')
text(35,.28,'r/B = 0.3')
text(40,.17,'r/B = 0.5')
text(50,.07,'r/B = 1')
text(15,.65,'Press Enter to Continue')
pause
t=[0 t];
framedata=[zeros(1,41);framedata];
figure(3);clf;% plot s vs r and t
axis([0 4 0 200 0 0.6]);
[R,T1]=meshgrid(roverB,t);
mesh(R,T1,framedata);
view(45,30)
xlabel('Distance from Well, r/B','Rotation',-20)
ylabel('Time, t [days]', 'Rotation',20)
zlabel('Dimensionless Drawdown, s(r,t)T/Q')
text(1,200,.5,'Press Enter to Continue')
pause
lim=max(max(framedata));

figure(4);clf;%plot animation of well drawdown as time evolves
r1=[-fliplr(r) r]/B;
framedata2=[fliplr(framedata) framedata];
h3=axes('Ydir','reverse');
axis([r1(1,1) r1(1,82) -.6*lim 3*lim]);
hold on
patch([r1(1,1) r1(1,82) r1(1,82) r1(1,1) r1(1,1)],[-.4*lim -.4*lim 0 0 -.4*lim],'m');
hold on
text(1,-.3*lim,['B = ' num2str(B) ' m']);
xpatch=[r1 r1(1,82) r1(1,1) r1(1,1)];
ypatch=[framedata2(1,:) lim lim framedata2(1,1)];
P1=patch(xpatch,ypatch,'c','EraseMode','background');
hold on
patch([r1(1,1) r1(1,82) r1(1,82) r1(1,1) r1(1,1)],[lim lim 3*lim 3*lim lim],'b');
hold on
plot([0 0 .5],[2.9*lim -.5*lim -.5*lim],'k','LineWidth',5);
hold on;
P2=plot([.5 .9 1.3],[-.5*lim -.45*lim -.42*lim],'b','LineWidth',5, 'EraseMode','background');
hold on
text(-3.5,-.1*lim,'Permeable Layer')
text(-3.5,.4*lim,'Clay-Aquitard')
text(-3.5,1.4*lim,'Aquifer')

T1=text(-3.5,-.3*lim,'Press Enter to Continue');
pause
set(T1,'String','               ');
for i=1:401
    ypatch=[framedata2(i,:) 1*lim 1*lim framedata2(i,1)];
    set(P1,'Ydata',ypatch);
    plot([0 0],[2.8*lim -.5*lim],'k','LineWidth',5);
    hold on
    pause(.02)
end
index2=[1 6 20 40];
for i=1:4
    plot(r1,framedata2(index(1,i),:));
    hold on
end
xlabel('Radial Distance from Well, r/B')
ylabel('Dimensionless Drawdown, s(r,t)T/Q')
text(-3.5,-.1*lim,'Permeable Layer')
text(1,.3,['Kprime = ' num2str(Kprime) ' m/day'])
text(1,.2,['bprime = ' num2str(bprime) ' m'])
text(1,2*lim,['S = ' num2str(S)])
text(1,1.6*lim,['b = ' num2str(b),' m'])
text(1,1.8*lim,['K = ' num2str(K) ' m/day'])
set(T1,'String','Press Enter to Continue');
set(gca,'Box','on');

pause
figure(5):clf;% plot 3-d plot of steady-state s(t)
for i=1:40
    for j=1:40
        r=sqrt(i^2+j^2)*.1;
        w(i,j)=(1/(2*pi))*besselk(0,r);
    end
end
y1=linspace(.1,4,40);
ax=(1/(2*pi))*besselk(0,y1);
w1=flipud(w);
w2=fliplr(w);
w3=fliplr(w1);
ax1=fliplr(ax);
ax2=ax';
ax3=flipud(ax2);
W=[w3 ax3 w1;ax1 ax(1,1) ax;w2,ax2 w];
x=linspace(-4,4,81);
y=x;
[X, Y]=meshgrid(x,y);
h5=axes('Zdir','reverse');
hold on
mesh(X,Y,W);
text(2,4,.35,'Distance, x/B','Rotation',33);
text(4,-2,.35,'Distance, y/B', 'Rotation',-11);
view(120,30);
zlabel('Drawdown, sT/Q')
grid on
text(-4,-1,.35,'Press Enter to Continue')
text(3,-4, .3,'Steady-State Cone of Depression')
pause