import numpy as np

g = np.array([1, 4/3, 4/3, 4/3, 4/3, 4/3, 4/3, 1])
b = np.array([2, 3, 1, 2, 3, 1, 1, 1])
p = np.array([1, 2, 2, 2, 2, 2, 2, 1])
h = np.empty_like(g)
sizing = np.array([3, 2, 2, 2, 2, 2, 2, 3])
min_sizing = np.array([3,4,4,4,4,4,4,3])

C_in = np.array([3,0,0,0,0,0,0,0,45], dtype=np.float32)

N = len(C_in) - 1
H = C_in[-1] / C_in[0]
G = np.prod(g)
B = np.prod(b)

P = np.sum(p)

for i in range(4):
    F = G*B*H
    f_opt = F**(1/N)
    D = N*f_opt + P
    
    for j in range(N)[::-1]:
        C_in[j] = g[j]*b[j]*C_in[j+1]/f_opt

    for j in range(N)[::-1]:
        C_in[j] = np.rint(C_in[j]/sizing[j])*sizing[j]
        if C_in[j] < min_sizing[j]:
            C_in[j] = min_sizing[j]
        h[j] = C_in[j+1] / C_in[j]
    H = np.prod(h)