построение нескольких массивов из цикла for

Question

Я новичок здесь (newb ie), поэтому вопрос, возможно, очевиден. Я закодировал форму мембраны. Теперь я изменяю начальный угол от 0 до 180 градусов с интервалом (1 градус) и строю мембраны на том же графике с вложенным l oop. Однако, если я нанесу мембраны под различным начальным углом с первым l oop, мембрана будет иметь другую форму, чем с вложенным l oop. Первый l oop показывает правильную форму, а вложенный l oop делает то, что мне не нужно.

Я не могу понять в коде, как это происходит.

import numpy as np
import matplotlib.pyplot as plt
import math

rho = 1000
g = 9.81
H_up = 5.15  # H - 2.7 # H = 7.85 #Water depth
H_down = 2.15  # D - 2.7 # D = 4.85
h_d = 5.3

p = h_d * rho * g
T = h_d/2 * p
alfa = h_d / (H_up-H_down)
T = 1/2 * alfa * rho * g * h_d**2
phi_0 = math.acos(rho * g / 2 * (H_up ** 2 - H_down ** 2) / T - 1)

def geometry(p, phi):
    rho = 1000
    g = 9.81
    H_up = 5.15  # H - 2.7
    H_down = 2.15  # D - 2.7
    h_d = 5.3  # [m] gate height

    alfa = h_d / (H_up-H_down)#1.5 # 1*10**10
    T = 1/2 * alfa * rho * g * h_d**2 #50000  # [N/m] 1/2 * rho * g * (H - D)**2 / (1 + np.cos(phi1)) #  h_d/2 * h_d * rho * g
    t = 1.6 * 10 ** -2  # [m] thickness sheet
    sigma = T / t #[N/m^2]
    E =  3800 * 10**3 /t #3200 / t * 10 ** 3  # [N/m**2] 900*10**6
    epsilon = sigma / E
    f = (1 + epsilon)

    dphi = - f * p / T
    dx = f * np.cos(phi)
    dy = f * np.sin(phi)
    return (dphi, dx, dy)

steps = 100
dS = 30./steps

phi = np.zeros(steps)
x = np.zeros(steps)
y = np.zeros(steps)
p = np.zeros(steps)
case = 1

for i in range (len(y)-1):
    phi[0] = phi_0  #geometry(p[i-1], phi[i])[3]
    phi[i] = phi[i-1] + geometry(p[i-1], 0)[0] * dS
    x[i] = x[i - 1] + geometry(p[i-1], phi[i])[1] * dS
    x[i+1] = x[i]
    y[i] = y[i - 1] + geometry(p[i-1], phi[i])[2] * dS
    y[i+1] = y[i]

    if y[i] > H_up and case == 1:
        case = 2
    if y[i] < H_down and case == 2:
        case = 3
    if y[i] < 0 and case == 3:
        break

    if case == 1:
        p[i] = (H_up - (H_up - y[i])) * rho * g
    if case == 2:
        p[i] = H_up * rho * g
    if case == 3:
        p[i] = (H_down - y[i]) * rho * g

plt.plot(x,y)
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()

phi1 = np.linspace(0, np.pi, 180)

phi = np.zeros(steps)
x = np.zeros(steps)
y = np.zeros(steps)
p = np.zeros(steps)
case = 1

for j in range (len(phi1)):
    phi2 = phi1[j]
    for i in range (len(y)-1):
        phi[0] = phi2 #geometry(p[i-1], phi[i])[3]
        phi[i] = phi[i-1] + geometry(p[i-1], phi[i])[0] * dS
        x[i] = x[i - 1] + geometry(p[i-1], phi[i])[1] * dS
        y[i] = y[i - 1] + geometry(p[i-1], phi[i])[2] * dS

        if y[i - 1] > H_up and case == 1:
            case = 2
        if y[i - 1] < H_down and case == 2:
            case = 3
        if y[i - 1] < 0 and case == 3 and max(y) > 10:
            plt.plot(x, y)
            phi = np.zeros(steps)
            x = np.zeros(steps)
            y = np.zeros(steps)
            p = np.zeros(steps)
            break

        if case == 1:
            p[i] = (H_up - (H_up - y[i])) * rho * g
        if case == 2:
            p[i] = H_up * rho * g  # upstream waterlevel - bottom # 40000  # [N/m^2]
        if case == 3:
            p[i] = (H_down - y[i]) * rho * g
plt.show()