I am trying to model a streamline with a pumping well. I have an equation describing the position of the streamline, however, the analytical solution fails at the well location. What I want is to locate the cells that misbehave in matrix U, which equal (xi+or-2) and (yi+or-2) and make them equal to xi and yi. I tried if (..) function but I havent been successful. Any aid is appreciated. Thanks, Eric
theta = 135; %degrees
Q = 10000; %ft^3/day
xi = 10; %Well X Position
yi = -10;%Well Y position
q = 1.5 B= 300Psi = @(m,d) q*sin(theta)*m - q*cos(theta)*d - Q./(2*pi*B)*atan(((m-xi)./(d-yi))); %
D = []; %matrix to contain results from Psi(x,y)
m = 9;%starting x position of particle
d = -9;%starting y position of particle
E= Psi(m,d); %evaluation of Psi for initial position
D = [D;E,m,d];%store Psi in D
%Psi2 equation used in the fzero equation
Psi2 = @(x,y) q*sin(theta)*x - q*cos(theta)*y - Q./(2*pi*B)*atan(((x-xi)./(y-yi))) - D(1,1);% where D(1,1) is the initial Psi value obtained from above.
U = [];%matrix created to contain results of fzero(Psi2).
for y = (d-25):(abs(d)+25); %Y values of particles
x0 = [-5000 5000]; %initial guess for fzero
[x1]= fzero(@(x) Psi2(x,y),x0); %fzero calculation
U=[U;x1, y]; end% if U(:,1) == (xi-2) && U(:,2) == (yi-2)
% U(:,1) = xi
% elseif U(:,1)>x1;
% x1=x1;
% else U(:,1)<x1;
% x1=x1;% end
plot(U(:,1), U(:,2), 'k')xlabel 'x (ft)'ylabel 'y (ft)'hold on plot (xi, yi,'r X')plot (m,d,'g O')
Best Answer