I am a structural engineer and trying to plot a moment-curvature interaction diagram via python. There, in the code, are so much iteration required. I am aware that Matlab is one of the best platforms in terms of scientific calculations and iterations. However, I want to learn the essentials of python while using my profession. I will briefly explain the code.

There are couple forces acting on the section; compression and tension. Compression forces are due to the concrete layer and the upper steel layer while the tensions are caused by the lower steel layers. The c is the neutral axis where these forces equal the sum of external forces which is also the axial force acting on the member. For the first case, it assumed that the total axial load (N) acting on the member is equal to zero. This will result that the sum of all compression and the tension forces must be zero. This could be achieved by manipulating the neutral axis depth, c. As the sum of forces closes to the intended value N=0, the related c value will determine all the values. Is there any problem that can be fixed in this code? Is it possible to make the code run faster? I will attach the code below, thank you.

```
def moment_curvature():
import numpy as np
import matplotlib.pyplot as plt
c_list = []
M_list = []
K_list = []
strain = np.linspace(0, 0.0038, 10)
for i in range(len(strain)):
ec_i = strain[i]
c_list.append(axis_depth(ec_i)[0])
M_list.append(axis_depth(ec_i)[1])
K_list.append(axis_depth(ec_i)[2])
plt.plot(K_list, M_list)
plt.xlabel('K (1/mm)')
plt.ylabel('M (kN.m)')
plt.show()
#######
def axis_depth(ec_i):
import numpy as np
###
h = 1000
bw = 500
fc = 40
#fcd = fc / 1.5
fs = 420
fyd = fs / 1.15
Es = 200000
ecu = 0.0038
e0 = 0.002
######
neutral_axis = np.linspace(1, h, h*100)
for k in range(len(neutral_axis)):
c = neutral_axis[k]
### concrete layers
sections = 100
delta = c / sections
area = delta * bw
### first layer of concrete: y_1, ec_1, fc_1, Fc_1
y_1 = delta / 2
ec_1 = ec_i * (1 - (y_1 / c))
if ec_1 <= e0:
fc_1 = fc * ((2 * ec_1 / e0)-(ec_1 / e0) ** 2)
elif e0 < ec_1 <= ecu:
fc_1 = fc * (1 - 0.15 * ((ec_1 - e0) / (ecu - e0 )))
Fc_1 = fc_1 * area * 1e-3
moment_arm = (h / 2) - y_1
moment_1 = Fc_1 * moment_arm * 1e-3
###
concrete_force = Fc_1
concrete_moment = moment_1
y = delta / 2
for i in range(sections - 1):
y += delta
ec_layer = ec_i * (1 - (y / c))
fc_layer = concrete_strain(ec_layer, e0, ecu, fc)
Fc_layer = fc_layer * area * 1e-3
concrete_force += Fc_layer
moment_arm = (h / 2) - y
moment_layer = Fc_layer * moment_arm * 1e-3
concrete_moment += moment_layer
### steel reinforcement properties: area, h_s, d_s, y_s
steel_area = [1256.64, 628.32, 1256.64]
#h_s = [50, 70, 950]
d_s = [950, 930, 50]
y_s = [-450, -430, 450]
es_1 = ec_i * (1 - (d_s[0] / c))
es_2 = ec_i * (1 - (d_s[1] / c))
es_3 = ec_i * (1 - (d_s[2] / c))
steel_stress = [es_1 * Es, es_2 * Es, es_3 * Es]
steel_force = 0
steel_moment = 0
for j in range(len(steel_stress)):
if steel_stress[j] >= 0:
if abs(steel_stress[j]) >= fyd:
steel_stress[j] = fyd
else:
if abs(steel_stress[j]) >= fyd:
steel_stress[j] = -fyd
steel_f = (steel_stress[j] * steel_area[j] * 1e-3)
steel_force += steel_f
steel_moment += (steel_f * y_s[j] * 1e-3)
###
N = 0
total_force = concrete_force + steel_force
total_moment = steel_moment + concrete_moment
curvature = ec_i / c
if (N - 0.1) < total_force < (N + 0.1):
break
return [c, total_moment, curvature]
######
def concrete_strain(ec_layer, e0, ecu, fc):
if ec_layer <= e0:
fc_layer = fc * ((2 * ec_layer / e0) - (ec_layer / e0) ** 2)
elif e0 < ec_layer <= ecu:
fc_layer = fc * (1 - 0.15 * ((ec_layer - e0) / (ecu - e0 )))
return fc_layer
######
moment_curvature()
```