using Pkg

function isinstalled(pkgname)
    deps = Pkg.dependencies()
    installed = false
    for (uuid, dep) in deps
        if dep.name == pkgname
            installed = true
            break
        end
    end
    return installed
end

if !isinstalled("PyPlot")
    println("PyPlot package is required but not found. Attempting to install ...")
    Pkg.add("PyPlot")
end
using PyPlot

"""
    potential(x, y[, U, R])

Compute the velocity potential for the incompressible inviscous 2D
flow around a cylinder. `x` and `y` are the cartesian coordinates of
the point where the pptential is calculated. (The origin of the
reference frame is at the center (axis) of the cylinder.) `R` is the
radius of the cylinder, `U` is the speed of the flow at infinity.

When the function is used as `potential(x, y)`, the default values `U
= 1.0` and `R = 1.0` are used. That is equivalent to using the
dimensionless coordinates `x/R` and `y/R` and measuring the potential
function in the units of `U*R`.

"""
function potential(x, y, U = 1.0, R = 1.0)
    r2 = x^2 + y^2
    R2 = R^2
    if (r2 < R2)
        # mask the inside of the cylinder
        psi = NaN
    else
        psi = U * (1 + R2 / r2) * x
    end
    return psi
end

"""
    stream(x, y[, U, R])

Compute the stream function for the incompressible inviscous 2D flow
around a cylinder. `x` and `y` are the cartesian coordinates of the
point where the stream function is calculated. (The origin of the
reference frame is at the center (axis) of the cylinder.) `R` is the
radius of the cylinder, `U` is the speed of the flow at infinity.

When the function is used as `stream(x, y)`, the default values `U =
1.0` and `R = 1.0` are used. That is equivalent to using the
dimensionless coordinates `x/R` and `y/R` and measuring the stream
function in the units of `U*R`.

"""
function stream(x, y, U = 1.0, R = 1.0)
    r2 = x^2 + y^2
    R2 = R^2
    if (r2 < R2)
        # mask the inside of the cylinder
        psi = NaN
    else
        psi = U * (1 - R2 / r2) * y
    end
    return psi
end

"""
    pressure(x, y[, U, R, rho])

Compute the pressure for the incompressible inviscous 2D flow around a
cylinder. `x` and `y` are the cartesian coordinates of the point where
the pressure is calculated. (The origin of the reference frame is at
the center (axis) of the cylinder.) `R` is the radius of the cylinder,
`rho` is the fluid density, `U` is the speed of the flow at infinity.

When the function is used as `pressure(x, y)`, the default values `U =
1.0`, `R = 1.0`, and `rho = 1.0` are used. That is equivalent to using
the dimensionless coordinates `x/R` and `y/R` and measuring the
pressure in the units of `rho*U^2`.

"""
function pressure(x, y, U = 1.0, R = 1.0, rho = 1.0)
    x2 = x^2
    y2 = y^2
    r2 = x2 + y2
    R2 = R^2
    if (r2 < R2)
        # mask the inside of the cylinder
        p = NaN
    else
        p = 1.0 / 2.0 * rho * U^2 * R2 / r2^2 * (2 * x2 - 2 * y2 - R2)
    end
    return p
end


# parameters

xmax = 3.0
ymax = 3.0
np = 200

# arrays

y = range(-ymax, ymax, length = np)
x = range(-xmax, xmax, length = np)
u = zeros(np, np)
v = zeros(np, np)
p = zeros(np, np)

# Calculate
for i = 1:np
    for j = 1:np
        u[i, j] = stream(x[j], y[i])
        v[i, j] = potential(x[j], y[i])
        p[i, j] = pressure(x[j], y[i])
    end
end

# Contour levels of stream function are the streamlines of the flow

# Plot the streamlines
# automatic selection of contour levels
ncont = 42
lws = 0.5
contour(x, y, u, ncont, colors = "black", linewidths = lws, linestyles="solid")
axis("scaled")
xlabel(L"$x/R$")
ylabel(L"$y/R$")

# Plot equipotential lines
ncont = 5
contour(x, y, v, ncont, colors = "black", linewidths = lws, linestyles="dotted")

# Plot the pressure
im = imshow(p, extent = [-ymax, ymax, -xmax, xmax], cmap = "bwr")
cbar = colorbar(im)
cbar.set_label(label = L"$\qquad\frac{P}{\rho \, U^2}$", rotation = 0, y = 0.54, fontsize = 14)
clim(-1.4, 1.4)

# Sketch the cylinder
R = 1.
phi = range(0.0, 2 * pi, length = np)
plot(R .* cos.(phi), R .* sin.(phi), linewidth = 2, color = "black")

if !isinteractive()
    extension = ".pdf"
    out = splitext(basename(@__FILE__))
    saveout = string(out[1], extension)
    savefig(saveout)
end