How to change my Mach number to stop at 2.89 instead of it going higher

so the graph I’m getting here has a line that goes to up pretty much double the amount that I need it to go to. It’s supposed to stop at 2.893 as that’s the highest Mach number my nozzle goes up to but can’t seem to get it to end there. Can someone please have a look at the snippet and tell me where I’m going wrong.

R_star = 0.2476
A_star = R_star^2 * π
gamma = 1.1706

x_values = P_X:0.01:1.5
y_values = R_throat .+ x_values.^2

M = zeros(length(x_values))

for i in 1:length(x_values)
    if x_values[i] == 0
        M0 = 1.0
    elseif x_values[i] < 0
        M0 = 0.5
    else
        M0 = 5
    end

    function f(M)
        return A_star / (y_values[i]^2 * π) * sqrt(((M^2 * γ - M^2 + 2) / (γ + 1))^((γ + 1) / (γ - 1))) - M
    end

    M[i] = find_zero(f, M0)
end

plot(x_values, M, xlabel="x", ylabel="Mach Number", title="Mach Analysis")

1 Like

Is there more code than what you provided here? You defined a variable gamma at the top but then reference a variable γ in your f(M) definition.

Other than that, nothing jumps out at me as an obvious error. It might help if you showed us a formatted/written version of the f(M) equation you’re trying to code to make sure that code actually says what you intended. The for i in 1:length loop with by-index writes/accesses is not very “Julian” code, so there’s probably a more elegant implementation that could be made, but that’s probably not central to the problem you’re currently trying to solve; feel free to ask for follow-up if interested in those tips.

Hello Mike, yes there’s some other coding lines above what I provided that include (\gamma) with the same value which is why i wrote it that way as its got a saved value in there. And can you explain what you meant by wanting to see a formatted/written version of the f(M) equation?

And sure i’d like to see what other more elegant/“Julian” way of writing the for ~~~ for i in 1:length ~~~ loop, as it would make my code better.

If I’m not mis-reading something, I’m interpreting your code as implementing this equation:
f(M) = \frac{A^*}{\pi y^2} \sqrt{ \left( \frac{\gamma M^2 - M^2 + 2}{\gamma+1} \right)^{(\gamma+1)/(\gamma-1)} } - M

Does this look right? It can be easy to make syntax errors in complicated equations.

I have to run right now but I’ll be back later with the other referenced tips.

Yes that’s how the function looks like when typing it on by hand.
And no worries take your time

So, in general, there’s a common code pattern in languages where you:

  • Establish an empty vector
  • Loop for i in 1:length(vector)
    • Do something to v[i]

Julia has really great affordances for array/vector operations and broadcasting and functional programming tools like map that can usually simplify this code and improve performance. This usually means defining a function that operates per element and then broadcasting/map’ing it over the entire vector at once. Example:

x_values = P_X:0.01:1.5
y_values = R_throat .+ x_values.^2

function M(x, y)
    # Transform your loop code into a single function of
    # the corresponding x and y values that determine it.
    return M
end

# This map generates pairs of (x,y) values from these input
# vectors and applies the M function to them, collecting the
# outputs into a new vector
M_values = map(M, x_values, y_values)
2 Likes