% Clear the Command Window and Workspace and close any open figures
clear;
clc;
close all;

% Define viewing size and step size
x1 = -5:0.01:-1.01;
x2 = -0.99:0.01:0.99;
x3 = 1.01:0.01:5;
y = -5:0.01:5;

% Other variables
V = 100;
R = 1;
Q = 3/2 * R * V;

% Create meshgrids
[X1, Y1] = meshgrid(x1, y);
[X2, Y2] = meshgrid(x2, y);
[X3, Y3] = meshgrid(x3, y);

% Calculate (nondimensionalized) thetas
theta1_1 = atan((Y1/R) ./ (X1/R + R));
theta2_1 = -atan((Y1/R) ./ (X1/R - R));
theta1_2 = atan((Y2/R) ./ (X2/R + R));
theta2_2 = -atan((Y2/R) ./ (X2/R - R));
theta1_3 = atan((Y3/R) ./ (X3/R + R));
theta2_3 = -atan((Y3/R) ./ (X3/R - R));

% Calculate stream function
psi1 = V*Y1/R + Q * (theta1_1 + theta2_1);
psi2 = V*Y2/R + Q * (theta1_2 + theta2_2);
psi3 = V*Y3/R + Q * (theta1_3 + theta2_3);

% Plot contours
figure;
hold on;
[~, g] = contour(X1, Y1, psi1, 25);
g.LineColor = 'black';
[~, g] = contour(X2, Y2, psi2, 25);
g.LineColor = 'black';
[~, g] = contour(X3, Y3, psi3, 25);
g.LineColor = 'black';