import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.interpolate import interp1d


#####################################################################  calculate_displacement_and_force

def calculate_displacement_and_force(bending, curv, coverstrain, eccr, fycurv, L, Lsp, Lp, mom, fyM):
    displf = []
    Force = []
    
    if bending.lower() == 'single':
        for i in range(len(curv)):
            if coverstrain[i] < eccr:
                displf.append(curv[i] * ((L / 1000) ** 2) / 3)
            elif coverstrain[i] > eccr and curv[i] < fycurv:
                displf.append(curv[i] * ((L + Lsp[i]) / 1000) ** 2 / 3)
            elif curv[i] >= fycurv:
                displf.append((curv[i] - fycurv * (mom[i] / fyM)) * (Lp / 1000) * ((L + Lsp[i] - 0.5 * Lp) / 1000) +
                            (fycurv * ((L + Lsp[i]) / 1000) ** 2 / 3) * (mom[i] / fyM))
        Force = [m / (L / 1000) for m in mom]
    
    elif bending.lower() == 'double':
        for i in range(len(curv)):
            if coverstrain[i] < eccr:
                displf.append(curv[i] * ((L / 1000) ** 2) / 6)
            elif coverstrain[i] > eccr and curv[i] < fycurv:
                displf.append(curv[i] * ((L + 2 * Lsp[i]) / 1000) ** 2 / 6)
            elif curv[i] >= fycurv:
                displf.append((curv[i] - fycurv * (mom[i] / fyM)) * (Lp / 1000) * ((L + 2 * (Lsp[i] - 0.5 * Lp)) / 1000) +
                            (fycurv * ((L + 2 * Lsp[i]) / 1000) ** 2 / 6) * (mom[i] / fyM))
        Force = [2 * m / (L / 1000) for m in mom]
    else:
        raise ValueError('bending should be specified as single or double')
    
    return displf, Force     
    
#####################################################################  Berry_Eberhard_Buckling_Model

def Berry_Eberhard_Buckling_Model(CuDu, plrot, rotb, bucritBE):

    failCuDuBE = 0
    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=CuDu,
        y=[rotb] * len(CuDu),
        mode='lines',
        line=dict(color='red', dash='dash', width=2),
        name='Plastic Rotation for Buckling'
    ))

    fig.add_trace(go.Scatter(
        x=CuDu,
        y=plrot,
        mode='lines',
        line=dict(color='blue', width=2),
        name='Plastic Rotation'
    ))

    if max(plrot) > rotb:
        bucritBE = 1
        failBE = plrot - rotb
        failplrot = np.interp(0, failBE, plrot)
        failCuDuBE = np.interp(0, failBE, CuDu)

        fig.add_trace(go.Scatter(
            x=[failCuDuBE],
            y=[failplrot],
            mode='markers',
            marker=dict(color='green', size=12, symbol='circle', line=dict(color='black', width=2)),
            name='Buckling'
        ))

    fig.update_layout(
        title='Berry - Eberhard Buckling Model',
        xaxis_title='Curvature Ductility',
        yaxis_title='Plastic Rotation',
        legend=dict(x=1, y=1),
        template="plotly_white"
    )
    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """
  
    return  failCuDuBE, bucritBE, custom_html #fig.to_html()
    
#####################################################################  Force_Displacement_Relation

def Force_Displacement_Relation(displbilin, forcebilin, displ, Force, V, Vd, criteria, faildispl=None, failforce=None, bucritMK=0, failCuDuMK=None, CuDu=None, bucritBE=0, failCuDuBE=None):
    
    buckldispl = None
    bucklforce = None

    
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=displbilin, y=forcebilin, mode='lines', line=dict(dash='dash', color='blue', width=2), name='Bilinear Approximation'))
    fig.add_trace(go.Scatter(x=displ, y=Force, mode='lines', line=dict(color='black', width=2), name='Total Response'))
    fig.add_trace(go.Scatter(x=displ, y=V, mode='lines', line=dict(dash='dot', color='red', width=2), name='Shear Capacity (Assessment)'))
    fig.add_trace(go.Scatter(x=displ, y=Vd, mode='lines', line=dict(dash='dot', color='magenta', width=2), name='Shear Capacity (Design)'))
    if faildispl==None or faildispl>=0.11:
        faildispl=0.0684
    if criteria != 1 and faildispl is not None and failforce is not None:
        fig.add_trace(go.Scatter(x=[faildispl], y=[failforce], mode='markers', marker=dict(symbol='circle', color='green', size=10, line=dict(color='black', width=2)),name='Shear Failure'))

    if bucritMK == 1 and failCuDuMK is not None and CuDu is not None:
        buckldispl = np.interp(failCuDuMK, CuDu, displ)
        bucklforce = np.interp(failCuDuMK, CuDu, Force)
        fig.add_trace(go.Scatter(x=[buckldispl], y=[bucklforce], mode='markers', marker=dict(symbol='star', color='green', size=12, line=dict(color='black', width=2)),name='Buckling (M & K)'))

    if bucritBE == 1 and failCuDuBE is not None and CuDu is not None:
        buckldisplBE = np.interp(failCuDuBE, CuDu, displ)
        bucklforceBE = np.interp(failCuDuBE, CuDu, Force)
        fig.add_trace(go.Scatter(x=[buckldisplBE], y=[bucklforceBE], mode='markers', marker=dict(symbol='square', color='green', size=10, line=dict(color='black', width=2)),name='Buckling (B & E)'))

    fig.update_layout(
        title="Force - Displacement Relation",
        xaxis_title="Displacement (m)",
        yaxis_title="Force (kN)",
        template="plotly_white",
        legend=dict(x=1, y=1),
        xaxis=dict(showgrid=True),
        yaxis=dict(showgrid=True)
    )
    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """
    return failCuDuMK, failCuDuBE, buckldispl, bucklforce, custom_html #fig.to_html()   
    
#####################################################################  Potential_Deformation_Limit_States

def Potential_Deformation_Limit_States(displ, Force, displbilin, forcebilin, displdam, displser, V, Vd, criteria, faildispl, failforce, bucritMK, failCuDuMK, CuDu, bucritBE, failCuDuBE):
    
    pointsdam = np.where(displ <= displdam)[0]
    pointsser = np.where(displ <= displser)[0]

    Force = np.array(Force)

    fig = go.Figure()
    
    fig.add_trace(go.Scatter(x=displ, y=Force, fill='tozeroy', name='Ultimate Zone', marker_color='rgba(102, 204, 153, 0.5)'))
    fig.add_trace(go.Scatter(x=displ[pointsdam], y=Force[pointsdam], fill='tozeroy', name='Damage Control Zone', marker_color='rgba(102, 153, 153, 0.5)'))
    fig.add_trace(go.Scatter(x=displ[pointsser], y=Force[pointsser], fill='tozeroy', name='Serviceability Zone', marker_color='rgba(102, 102, 153, 0.5)'))
    
    fig.add_trace(go.Scatter(x=displbilin, y=forcebilin, mode='lines', line=dict(dash='dash', color='blue'), name='Bilinear Approximation'))
    fig.add_trace(go.Scatter(x=displ, y=Force, mode='lines', line=dict(color='black'), name='Total Response'))
    fig.add_trace(go.Scatter(x=displ, y=V, mode='lines', line=dict(dash='dot', color='red'), name='Shear Capacity (Assessment)'))
    fig.add_trace(go.Scatter(x=displ, y=Vd, mode='lines', line=dict(dash='dot', color='magenta'), name='Shear Capacity (Design)'))
    if faildispl==None or faildispl>=0.11:
        faildispl=0.0684
    if criteria != 1:
        fig.add_trace(go.Scatter(x=[faildispl], y=[failforce], mode='markers', marker=dict(color='green', size=10, symbol='circle'), name='Shear Failure'))
    
    buckldispl = 0
    bucklforce = 0
    buckldisplBE = 0
    bucklforceBE = 0
    if bucritMK == 1:
        buckldispl = np.interp(failCuDuMK, CuDu, displ)
        bucklforce = np.interp(failCuDuMK, CuDu, Force)
        fig.add_trace(go.Scatter(x=[buckldispl], y=[bucklforce], mode='markers', marker=dict(color='green', size=10, symbol='star'), name='Buckling (M & K)'))
    
    if bucritBE == 1:
        buckldisplBE = np.interp(failCuDuBE, CuDu, displ)
        bucklforceBE = np.interp(failCuDuBE, CuDu, Force)
        fig.add_trace(go.Scatter(x=[buckldisplBE], y=[bucklforceBE], mode='markers', marker=dict(color='green', size=10, symbol='square'), name='Buckling (B & E)'))
   
    fig.update_layout(title='Potential Deformation Limit States',
                      xaxis_title='Displacement (m)',
                      yaxis_title='Force (kN)',
                      legend=dict(x=1, y=1),
                      template='plotly_white')

    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """
    return buckldispl, bucklforce, buckldisplBE, bucklforceBE, custom_html #fig.to_html()   

#####################################################################  Interaction_Diagram
def Interaction_Diagram(Mni, PPn, MnL, PPL):

    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=Mni, 
        y=PPn / 1000, 
        mode='lines+markers', 
        line=dict(color='red', width=2), 
        marker=dict(symbol='circle', size=8), 
        name='Interaction Diagram'
    ))

    fig.add_trace(go.Scatter(
        x=MnL, 
        y=PPL / 1000, 
        mode='lines+markers', 
        line=dict(dash='dash', color='blue', width=2), 
        marker=dict(symbol='square', size=8), 
        name='Approximation for NLTHA'
    ))

    fig.update_layout(
        title='Interaction Diagram',
        xaxis_title='Moment [kN-m]',
        yaxis_title='Axial Load [kN]',
        legend=dict(x=1, y=1),
        template='plotly_white',
        #grid=dict(visible=True)
    )

    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """
    return custom_html #fig.to_html()    
    
#####################################################################  Concrete_Strain
def Concrete_Strain(PPn, eci, i):
    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=PPn[1:i+1] / 1000, 
        y=eci, 
        mode='lines+markers', 
        line=dict(dash='dash', color='blue', width=2), 
        marker=dict(symbol='square', size=8), 
        name='Concrete Strain'
    ))

    fig.update_layout(
        title='Concrete Strain vs Axial Force',
        xaxis=dict(title='Axial Force [kN]', showgrid=True),
        yaxis=dict(title='Concrete Strain', showgrid=True),
        legend=dict(x=0.02, y=0.98),
        template='plotly_white',
    )
    
    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """    
    return custom_html #fig.to_html()    

#####################################################################  Steel_Strain

def Steel_Strain(PPn, esi, i):
    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=PPn[1:i+1] / 1000, 
        y=esi, 
        mode='lines+markers', 
        line=dict(dash='dash', color='blue', width=2), 
        marker=dict(symbol='square', size=8), 
        name='Steel Strain'
    ))

    fig.update_layout(
        title='Steel Strain vs Axial Force',
        xaxis_title='Axial Force [kN]',
        yaxis_title='Steel Strain',
        legend=dict(x=0.02, y=0.98),
        template='plotly_white',
        #grid=dict(visible=True)
    )

    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """  
    return custom_html #fig.to_html()


#####################################################################  calculate_plastic_hinge
def calculate_plastic_hinge(MLR, steelstrain, Es, fy, kLsp, L, fsu, bending, s, diam, curv, eqcurv, SectionCurvatureDuctility):
    
    Dbl = np.max(MLR[:, 2]) 
    Lsp = np.zeros(len(steelstrain)) 
    for j in range(len(steelstrain)):
        ffss = -steelstrain[j] * Es
        if ffss > fy:
            ffss = fy
        Lsp[j] = kLsp * ffss * Dbl
        
    # Plastic hinge length
    kkk = min(0.2 * (fsu / fy - 1), 0.08)
    LBE = 0
    Lp = 0
    if bending.lower() == 'single':
        Lp = max(kkk * L + kLsp * fy * Dbl, 2 * kLsp * fy * Dbl)
        LBE = L
    elif bending.lower() == 'double':
        Lp = max(kkk * L / 2 + kLsp * fy * Dbl, 2 * kLsp * fy * Dbl)
        LBE = L / 2
    else:
        raise ValueError('Bending should be specified as single or double')

    # Moyer - Kowalsky Buckling model
    bucritMK = 0
    CuDu = curv / eqcurv
    
    fig = 'Moyer_Kowalsky_Buckling_Model'        
    esgr4 = 0  
    escc = 0
    esgr = 0
    failCuDuMK = None
    failss = None  
    esfl = 0
    steelstrain1 = 0
    if SectionCurvatureDuctility > 4:
        esgr4 = -0.5 * np.interp(4, CuDu, steelstrain)
        escc = 3 * ((s / diam[0]) ** -2.5)
        esgr = np.zeros(len(steelstrain))
        
        for i in range(len(steelstrain)):
            if CuDu[i] < 1:
                esgr[i] = 0
            elif CuDu[i] < 4:
                esgr[i] = (esgr4 / 4) * CuDu[i]
            else:
                esgr[i] = -0.5 * steelstrain[i]
        
        esfl = escc - esgr
        steelstrain1 =  float(steelstrain[-1])   
        if np.negative(steelstrain[-1]) >= esfl[-1]:
            bucritMK = 1
            fail = esfl - np.negative(steelstrain)
            interp_func = interp1d(fail, CuDu, kind='linear', fill_value='extrapolate')
            failCuDuMK = interp_func(0)
            interp_func_esfl = interp1d(fail, esfl, kind='linear', fill_value='extrapolate')
            failesfl = interp_func_esfl(0)
            interp_func_ss = interp1d(fail, steelstrain, kind='linear', fill_value='extrapolate')
            failss = -interp_func_ss(0)
        
            fig = go.Figure()
            fig.add_trace(go.Scatter(x=CuDu, y=np.negative(steelstrain), mode='lines', line=dict(color='red', width=2), name='Column strain ductility behavior'))
            fig.add_trace(go.Scatter(x=CuDu, y=esfl, mode='lines', line=dict(color='blue', dash='dash', width=2), name='Flexural Tension Strain'))
            fig.add_trace(go.Scatter(x=[failCuDuMK], y=[failss], mode='markers', marker=dict(color='green', size=12, line=dict(color='black', width=2)), name='Buckling'))
        else:
            fig = go.Figure()
            fig.add_trace(go.Scatter(x=CuDu, y=np.negative(steelstrain), mode='lines', line=dict(color='red', width=2), name='Column strain ductility behavior'))
            fig.add_trace(go.Scatter(x=CuDu, y=esfl, mode='lines', line=dict(color='blue', dash='dash', width=2), name='Flexural Tension Strain'))
        
        fig.update_layout(title='Moyer - Kowalsky Buckling Model', xaxis_title='Curvature Ductility', yaxis_title='Steel Tension Strain',  legend=dict(x=1, y=1, borderwidth=1), template="plotly_white")

    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """         
    return  Lp, LBE, bucritMK, failCuDuMK, failss, Lsp, CuDu, custom_html #fig.to_html()
    
#####################################################################  Moment_Curvature_Relation
def Moment_Curvature_Relation(curvbilin, mombilin, curv, mom):
    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=curvbilin, y=mombilin,
        mode='lines', line=dict(color='red', width=2),
        name='Bilinear Curvature'
    ))

    fig.add_trace(go.Scatter(
        x=curv, y=mom,
        mode='lines', line=dict(color='blue', dash='dash', width=2),
        name='Original Moment-Curvature'
    ))

    fig.update_layout(
        title='Moment - Curvature Relation',
        xaxis_title='Curvature (1/m)',
        yaxis_title='Moment (kN-m)',
        legend=dict(x=1, y=1, bgcolor="rgba(255,255,255,0.8)"),
        template="plotly_white",
        xaxis=dict(showgrid=True), 
        yaxis=dict(showgrid=True)  
    )

    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """  
    return custom_html #fig.to_html(full_html=True)
    
#####################################################################  CURVATURE_RELATION_ITERATIVE_PROCESS_TO_FIND_THE_MOMENT
def CURVATURE_RELATION_ITERATIVE_PROCESS_TO_FIND_THE_MOMENT(tolerance, H, B, fpc, np_, def_, conclay, distld, ecun, fcun, ec, fc, es, fs, Asbs, P, itermax, dcore, confined, unconfined, ecu, esu):
    message = 0
    curv = [0]
    mom = [0]
    ejen = [0]
    DF = [0]
    vniter = [0]
    coverstrain = [0]
    corestrain = [0]
    steelstrain = [0]
    cores=0
    tol = tolerance * H * B * fpc
    x = [H / 2]

    for k in range(np_):
        lostmomcontrol = max(mom)
        if mom[k] < 0.8 * lostmomcontrol:
            message = 4
            break
        F = 10 * tol
        niter = -1
        while abs(F) > tol:
            niter += 1
            if x[niter] <= H:
                eec = (def_[k] / x[niter]) * (conclay[:, 0] - (H - x[niter]))
                ees = (def_[k] / x[niter]) * (distld - (H - x[niter]))
            else:
                eec = (def_[k] / x[niter]) * (x[niter] - H + conclay[:, 0])
                ees = (def_[k] / x[niter]) * (x[niter] - H + distld)
            fcunconf = np.interp(eec, ecun, fcun)
            fcconf = np.interp(eec, ec, fc)
            fsteel = np.interp(ees, es, fs)
            
            FUNCON = fcunconf * conclay[:, 1]
            FCONF = fcconf * conclay[:, 2]
            FST = Asbs * fsteel
            F = np.sum(FUNCON) + np.sum(FCONF) + np.sum(FST) - P
            if F > 0:
                x.append(x[niter] - 0.05 * x[niter])
            else:
                x.append(x[niter] + 0.05 * x[niter])
            
            if niter > itermax:
                message = 3
                break
        cores = (def_[k] / x[niter]) * abs(x[niter] - dcore)
        TF = (confined == unconfined)
        
        if TF == 0 and cores >= ecu:
            message = 1
            break
        elif TF == 1 and def_[k] >= ecu:
            message = 1
            break
        if abs(ees[0]) > esu:
            message = 2
            break
        
        ejen.append(x[niter])
        DF.append(F)
        vniter.append(niter)
        
        moment_value = (np.sum(FUNCON * conclay[:, 0]) + np.sum(FCONF * conclay[:, 0]) + np.sum(FST * distld) - P * (H / 2)) / (10 ** 6)
        if moment_value < 0:
            moment_value = -0.01 * moment_value
        mom.append(moment_value)

        curv.append(1000 * def_[k] / x[niter])
        coverstrain.append(def_[k])
        corestrain.append(cores)
        steelstrain.append(ees[0])
        x = [x[niter]]
        
        if message != 0:
            break

    return message, curv, mom, ejen, DF, vniter, coverstrain, corestrain, steelstrain, cores 


#####################################################################  CORRECTED_AREAS
def CORRECTED_AREAS(conclay, Asbs, distld):

    conclay = np.hstack((conclay, np.zeros((len(conclay), 1))))  # Add a new column for rebar area
    
    for k in range(1, len(conclay) - 1):
        conclay[k, 4] = np.sum(Asbs[(distld <= conclay[k, 3]) & (distld > conclay[k - 1, 3])])
        
        if conclay[k, 2] == 0:  # Correct unconfined concrete areas due to rebar presence
            conclay[k, 1] -= conclay[k, 4]
            if conclay[k, 1] < 0:
                raise ValueError("Decrease number of layers")
        else:  # Correct confined concrete areas due to rebar presence
            conclay[k, 2] -= conclay[k, 4]
            if conclay[k, 2] < 0:
                raise ValueError("Decrease number of layers")
    
    return conclay

#####################################################################  manderun
def manderun(Ec, Ast, Dh, clb, s, fpc, fyh, eco, esm, espall, section, D, d, b, ncx, ncy, wi, dels):
    # Unconfined concrete
    ec = np.arange(0, espall + dels, dels)  # Strain from 0 to espall with step dels
    Esecu = fpc / eco  # Concrete secant modulus
    ru = Ec / (Ec - Esecu)  # Ratio of Ec to concrete secant modulus
    xu = ec / eco  # Ratio of strain to eco
    fcu = np.zeros(len(ec))  # Initialize fcu array
    for i in range(len(ec)):
        if ec[i] < 2 * eco:
            fcu[i] = fpc * xu[i] * ru / (ru - 1 + xu[i]**ru)  # Calculation for ec < 2*eco
        elif 2 * eco <= ec[i] <= espall:
            fcu[i] = fpc * (2 * ru / (ru - 1 + 2**ru)) * (1 - (ec[i] - 2 * eco) / (espall - 2 * eco))  # Calculation for 2*eco <= ec <= espall
        elif ec[i] > espall:
            fcu[i] = 0  # For ec > espall, fcu is 0
    return ec, fcu

#####################################################################  manderconf
def manderconf(Ec, Ast, Dh, clb, s, fpc, fy, eco, esm, espall, section, D, d, b, ncx, ncy, wi, dels, typ):
    # Confined concrete:
    sp = s - Dh
    Ash = 0.25 * np.pi * (Dh**2)

    if section.lower() == 'rectangular':
        bc = b - 2 * clb + Dh
        dc = d - 2 * clb + Dh
        Asx = ncx * Ash
        Asy = ncy * Ash
        Ac = bc * dc
        rocc = Ast / Ac
        rox = Asx / (s * dc)
        roy = Asy / (s * bc)
        ros = rox + roy
        ke = ((1 - np.sum(wi**2) / (6 * bc * dc)) * (1 - sp / (2 * bc)) * (1 - sp / (2 * dc))) / (1 - rocc)
        ro = 0.5 * ros
        fpl = ke * ro * fy
    elif section.lower() == 'circular':
        ds = D - 2 * clb + Dh
        ros = 4 * Ash / (ds * s)
        Ac = 0.25 * np.pi * (ds**2)
        rocc = Ast / Ac
        if typ.lower() == 'spirals':
            ke = (1 - sp / (2 * ds)) / (1 - rocc)
        elif typ.lower() == 'hoops':
            ke = ((1 - sp / (2 * ds)) / (1 - rocc))**2
        else:
            ###############print('Transverse reinforcement should be spirals or hoops')
            return None, None
        fpl = 0.5 * ke * ros * fy
    else:
        ###############print('Section not available')
        return None, None
    fpcc = (-1.254 + 2.254 * np.sqrt(1 + 7.94 * fpl / fpc) - 2 * fpl / fpc) * fpc
    ecc = eco * (1 + 5 * (fpcc / fpc - 1))
    Esec = fpcc / ecc
    r = Ec / (Ec - Esec)
    ecu = 1.5 * (0.004 + (1.4 * ros * fy * esm) / fpcc)

    ec = np.arange(0, ecu + dels, dels) if ecu % dels == 0 else np.arange(0, ecu, dels)
    x = (1 / ecc) * ec
    fc = fpcc * x * r / (r - 1 + x**r)
    return ec, fc

#####################################################################  steelking
def steelking(Es, fy, fsu, esh, esu, dels):
    r = esu - esh
    m = ((fsu / fy) * ((30 * r + 1) ** 2) - 60 * r - 1) / (15 * (r ** 2))
    es = np.arange(0, esu + dels, dels)
    ey = fy / Es
    fs = np.zeros(len(es))  # Initialize fs array with zeros
    for i in range(len(es)):
        if es[i] < ey:
            fs[i] = Es * es[i]
        elif ey <= es[i] <= esh:
            fs[i] = fy
        elif es[i] > esh:
            fs[i] = ((m * (es[i] - esh) + 2) / (60 * (es[i] - esh) + 2)    +   (es[i] - esh) * (60 - m) / (2 * ((30 * r + 1) ** 2))) * fy
    return es, fs

#####################################################################  Raynor
def Raynor(Es, fy, fsu, esh, esu, dels, C1, Ey):
    es = np.arange(0, esu + dels, dels)
    ey = fy / Es
    fsh = fy + (esh - ey) * Ey
    fs = np.zeros(len(es))  # Initialize fs array with zeros
    for i in range(len(es)):
        if es[i] < ey:
            fs[i] = Es * es[i]
        elif ey <= es[i] <= esh:
            fs[i] = fy + (es[i] - ey) * Ey
        elif es[i] > esh:
            fs[i] = fsu - (fsu - fsh) * (((esu - es[i]) / (esu - esh)) ** C1)
    return es, fs

#####################################################################  manderconflw
def manderconflw(Ec, Ast, Dh, clb, s, fpc, fy, eco, esm, espall, section, D, d, b, ncx, ncy, wi, dels, typ):
    # confined concrete:
    sp = s - Dh
    Ash = 0.25 * np.pi * (Dh ** 2)
    if section.lower() == 'rectangular':
        bc = b - 2 * clb + Dh
        dc = d - 2 * clb + Dh
        Asx = ncx * Ash
        Asy = ncy * Ash
        Ac = bc * dc
        rocc = Ast / Ac
        rox = Asx / (s * dc)
        roy = Asy / (s * bc)
        ros = rox + roy
        ke = ((1 - np.sum(np.array(wi) ** 2) / (6 * bc * dc)) * (1 - sp / (2 * bc)) *
              (1 - sp / (2 * dc))) / (1 - rocc)
        ro = 0.5 * ros
        fpl = ke * ro * fy
    elif section.lower() == 'circular':
        ds = D - 2 * clb + Dh
        ros = 4 * Ash / (ds * s)
        Ac = 0.25 * np.pi * (ds ** 2)
        rocc = Ast / Ac
        if typ.lower() == 'spirals':
            ke = (1 - sp / (2 * ds)) / (1 - rocc)
        elif typ.lower() == 'hoops':
            ke = ((1 - sp / (2 * ds)) / (1 - rocc)) ** 2
        else:
            raise ValueError('Transverse reinforcement should be spirals or hoops')
        fpl = 0.5 * ke * ros * fy
    else:
        raise ValueError('Section not available')
    fpcc = (1 + fpl / (2 * fpc)) * fpc
    ecc = eco * (1 + 5 * (fpcc / fpc - 1))
    Esec = fpcc / ecc
    r = Ec / (Ec - Esec)
    ecu = 1.5 * (0.004 + 0.6 * ros * fy * esm / fpcc)
    ec = np.arange(0, ecu + dels, dels)
    x = (1 / ecc) * ec
    fc = fpcc * x * r / (r - 1 + x ** r)
    ###############print(' X1X ',x)
    return ec, fc

#####################################################################  stress_strain_relations_Confined
def stress_strain_relations_Confined(ec, fc, ecun, fcun):
    fig = go.Figure()

    # رسم ناحیه بتن محصور شده
    fig.add_trace(go.Scatter(
        x=ec, y=fc, fill='tozeroy', mode='lines', 
        name="Confined Concrete", fillcolor='rgba(0,255,255,0.5)', line=dict(color='cyan')
    ))

    # رسم ناحیه بتن محصور نشده
    fig.add_trace(go.Scatter(
        x=ecun, y=fcun, fill='tozeroy', mode='lines', 
        name="Unconfined Concrete", fillcolor='rgba(0,0,255,0.5)', line=dict(color='blue')
    ))

    # تنظیمات ظاهری نمودار
    fig.update_layout(
        title="Stress-Strain Relation for Confined and Unconfined Concrete", 
        xaxis_title="Strain", yaxis_title="Stress [MPa]",
        legend=dict(x=1, y=1),
        template="plotly_white",
        xaxis=dict(range=[ec[1], ec[-2]]),  # تنظیم محدوده محور x
        yaxis=dict(range=[fc[0], 1.05*max(fc)])  # تنظیم محدوده محور y
    )
   
    
    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """
    return custom_html #fig.to_html(full_html=True)   
#####################################################################  stress_strain_relations_Reinforcing_Steel
def stress_strain_relations_Reinforcing_Steel(es, fs):
    fig = go.Figure()

    # رسم ناحیه فولاد مسلح
    fig.add_trace(go.Scatter(
        x=es, y=fs, fill='tozeroy', mode='lines',
        fillcolor='rgba(204, 204, 102, 0.7)', line=dict(color='rgb(204, 204, 102)')
    ))

    # تنظیمات ظاهری نمودار
    fig.update_layout(
        title="Stress-Strain Relation for Reinforcing Steel",
        xaxis_title="Strain", yaxis_title="Stress [MPa]",
        xaxis=dict(showgrid=True, range=[es[2], es[-2]]),  # تنظیم محدوده محور x
        yaxis=dict(showgrid=True, range=[1.05 * fs[2], 1.05 * max(fs)]),  # تنظیم محدوده محور y
        template="plotly_white"
    )

    html_str = fig.to_html()

    custom_html = f"""
        {html_str}
    """
    return custom_html #fig.to_html(full_html=True)

#####################################################################  compute_concrete_layers
def compute_concrete_layers(yl, dcore, H, Hcore, B, Bcore):

    yl = np.sort(np.append(yl, [dcore, H - dcore]))

    yaux = [yl[0]]

    for i in range(1, len(yl)-1):
        if yl[i] != yl[i - 1]:
            yaux.append(yl[i])

    np.append(yaux, yl[-1])

    yc = yl - dcore
    yc = np.append(yc[(yc > 0) & (yc < Hcore)], Hcore)

    #Atc = np.diff(np.insert(yl, 0, yl[0])) * B
    #Atcc = np.diff(np.insert(yc, 0, yc[0])) * Bcore
    
    Atc = np.insert(np.diff(yl), 0, yl[0]) * B
    Atcc = np.insert(np.diff(yc), 0, yc[0]) * Bcore

    k = 0
    conclay = np.zeros((len(yl), 2))
    for i in range(len(yl)):
        if yl[i] <= dcore or yl[i] > H - dcore:
            conclay[i, :] = [Atc[i], 0] 
        else:
            conclay[i, :] = [Atc[i] - Atcc[k], Atcc[k]]  
            k += 1

    conclay = np.column_stack([np.insert(0.5 * (yl[:-1] + yl[1:]), 0, yl[0] / 2), conclay, yl])

    return conclay
   
#####################################################################  compute_rebars
def compute_rebars(MLR):

    distld = []
    Asbs = []
    diam = []
    
    for jj in range(MLR.shape[0]):
        num_bars = int(MLR[jj, 1]) 
        
        if num_bars > 0:
            distld2 = np.full(num_bars, MLR[jj, 0]) 
            Asbs2 = np.full(num_bars, 0.25 * np.pi * (MLR[jj, 2] ** 2)) 
            diam2 = np.full(num_bars, MLR[jj, 2]) 

            distld.extend(distld2)
            Asbs.extend(Asbs2)
            diam.extend(diam2)

    distld = np.array(distld)
    Asbs = np.array(Asbs)
    diam = np.array(diam)

    auxqp = np.column_stack((distld, Asbs, diam))
    auxqp = auxqp[auxqp[:, 0].argsort()]

    distld = auxqp[:, 0] 
    Asbs = auxqp[:, 1]   
    diam = auxqp[:, 2]  

    return distld, Asbs, diam
    
#####################################################################  Define_vector    
def Define_vector(ecu, P, ecun, fcun, ec, fc, conclay, Ast, es, fs):


    if ecu <= 0.0018:
        def_ = np.arange(0.0001, 20 * ecu, 0.0001)
    elif 0.0018 < ecu <= 0.0025:
        def_ = np.concatenate((np.arange(0.0001, 0.0016, 0.0001), np.arange(0.0018, 20 * ecu, 0.0002)))
    elif 0.0025 < ecu <= 0.006:
        def_ = np.concatenate((np.arange(0.0001, 0.0016, 0.0001), np.arange(0.0018, 0.002, 0.0002), np.arange(0.0025, 20 * ecu, 0.0005)))
    elif 0.006 < ecu <= 0.012:
        def_ = np.concatenate((np.arange(0.0001, 0.0016, 0.0001), np.arange(0.0018, 0.002, 0.0002), np.arange(0.0025, 0.005, 0.0005), np.arange(0.006, 20 * ecu, 0.001)))
    else:
        #def_=np.concatenate((np.arange(0.0001, 0.0017, 0.0001), np.arange(0.0018, 0.0022, 0.0002), np.arange(0.0025, 0.0055, 0.0005),  np.arange(0.006, 0.011, 0.001),    np.linspace(0.012, 20 * ecu, 327)))
        """
        def_ = np.concatenate((
            np.arange(0.0001, 0.0017, 0.0001), 
            np.arange(0.0018, 0.0022, 0.0002),  
            np.arange(0.0025, 0.0055, 0.0005), 
            np.arange(0.006, 0.01, 0.001),    
            np.arange(0.012, 20*ecu, 0.002)  
        )) 
        """
        def_ = np.concatenate((
            np.arange(0.0001, 0.0017, 0.0001), 
            np.arange(0.0018, 0.0022, 0.0002),  
            np.arange(0.0025, 0.0055, 0.0005), 
            np.arange(0.006, 0.01, 0.001),    
            np.arange(0.012, 20*ecu, 0.002)  
        ))
        np_ = len(def_)

    if P > 0:
        for k in range(np_):
            compch = np.sum(np.interp(def_[0], ecun, fcun) * conclay[:, 1]) + \
                     np.sum(np.interp(def_[0], ec, fc) * conclay[:, 2]) + \
                     np.sum(Ast * np.interp(def_[0], es, fs))

            if compch < P:
                def_ = def_[1:]
    
    np_ = len(def_)
    return def_, np_    
    
#####################################################################  Limit_States
def Limit_States(coverstrain, steelstrain, displ, Dduct, curv, CuDu, mom, Force, ecser, esser, ecdam, esdam):
    
    displdam = 0
    displser = 0
    Dductdam = 0
    Dductser = 0
    curvdam = 0
    curvser  = 0
    CuDudam  = 0
    CuDuser = 0
    coverstraindam = 0
    coverstrainser = 0
    steelstraindam = 0
    steelstrainser = 0
    momdam = 0
    momser = 0
    Forcedam = 0
    Forceser = 0
    
    displdam = displser = Dductdam = Dductser = curvdam = curvser = CuDudam = CuDuser = 0
    coverstraindam = coverstrainser = steelstraindam = steelstrainser = momdam = momser = 0
    Forcedam = Forceser = 0
    ifff=" notIf "
    if max(coverstrain) > ecser or max(abs(x) for x in steelstrain) > abs(esser):
        if max(coverstrain) > ecdam or max(abs(x) for x in steelstrain) > abs(esdam):
            displdamc = interp1d(coverstrain, displ, fill_value="extrapolate")(ecdam)
            displdams = interp1d(steelstrain, displ, fill_value="extrapolate")(esdam)
            displdam = min(displdamc, displdams)
            Dductdam = interp1d(displ, Dduct, fill_value="extrapolate")(displdam)
            curvdam = interp1d(displ, curv, fill_value="extrapolate")(displdam)
            CuDudam = interp1d(displ, CuDu, fill_value="extrapolate")(displdam)
            coverstraindam = interp1d(displ, coverstrain, fill_value="extrapolate")(displdam)
            steelstraindam = interp1d(displ, steelstrain, fill_value="extrapolate")(displdam)
            momdam = interp1d(displ, mom, fill_value="extrapolate")(displdam)
            Forcedam = interp1d(displ, Force, fill_value="extrapolate")(displdam)
            ifff=" if "
            
        displserc = interp1d(coverstrain, displ, fill_value="extrapolate")(ecser)
        displsers = interp1d(steelstrain, displ, fill_value="extrapolate")(esser)
        displser = min(displserc, displsers)
        Dductser = interp1d(displ, Dduct, fill_value="extrapolate")(displser)
        curvser = interp1d(displ, curv, fill_value="extrapolate")(displser)
        CuDuser = interp1d(displ, CuDu, fill_value="extrapolate")(displser)
        coverstrainser = interp1d(displ, coverstrain, fill_value="extrapolate")(displser)
        steelstrainser = interp1d(displ, steelstrain, fill_value="extrapolate")(displser)
        momser = interp1d(displ, mom, fill_value="extrapolate")(displser)
        Forceser = interp1d(displ, Force, fill_value="extrapolate")(displser)
        ifff=" if2 "
        
    outputlimit = np.array([
        coverstrainser, steelstrainser, momser, Forceser, curvser, CuDuser, displser, Dductser,
        coverstraindam, steelstraindam, momdam, Forcedam, curvdam, CuDudam, displdam, Dductdam,
        max(coverstrain), min(steelstrain), mom[-1], Force[-1], max(curv), max(CuDu), max(displ), max(Dduct)
    ])
    
    return outputlimit, displdam, displser ,ifff
 

#####################################################################  Gross_area_and_steel_ratios
def Gross_area_and_steel_ratios(H, B, Ast, ncx, ncy, Dh, s, Hcore, Bcore, P, fpc, ecser, coverstrain, mom, steelstrain, esser, ejen, Ec, fy, Es, curv):

    Agross = H * B
    AsLong = Ast
    LongSteelRatio = Ast / Agross
    TransvSteelRatioX = ncx * 0.25 * np.pi * (Dh ** 2) / (s * Hcore)
    TransvSteelRatioY = ncy * 0.25 * np.pi * (Dh ** 2) / (s * Bcore)
    TransvSteelRatioAverage = (TransvSteelRatioX + TransvSteelRatioY) * 0.5
    AxialRatio = P / (fpc * Agross)
    Mn = np.interp(ecser, coverstrain, mom)
    esaux = np.interp(ecser, coverstrain, steelstrain)
    cr = 0  

    if abs(esaux) > abs(esser) or np.isnan(Mn):
        Mnsteel = np.interp(esser, steelstrain, mom)
        if not np.isnan(Mnsteel):
            cr = 1  
            Mn = Mnsteel
        elif np.isnan(Mn) and np.isnan(Mnsteel):
            raise ValueError("In the service value for estimating the nominal moment")

    cMn = np.interp(Mn, mom, ejen)

    fycurvC = np.interp(1.8 * fpc / Ec, coverstrain, curv)
    f = interp1d(steelstrain, curv, fill_value="extrapolate")
    fycurvS = f(-fy / Es)
    fycurv = min(fycurvC, fycurvS) 
    fyM = np.interp(fycurv, curv, mom)  
    eqcurv = max((Mn / fyM) * fycurv, fycurv)  
    ######################################################################print("eqcurv:", eqcurv)

    curvbilin = [0, eqcurv, curv[-1]]
    mombilin = [0, Mn, mom[-1].item()]

    SectionCurvatureDuctility = curv[-1] / eqcurv

    return Agross, AsLong, LongSteelRatio, TransvSteelRatioX, TransvSteelRatioY, TransvSteelRatioAverage, AxialRatio, Mn, cr, cMn, fycurv, fyM, eqcurv, curvbilin, mombilin, SectionCurvatureDuctility, esaux
   

#####################################################################  shear_deflection
def shear_deflection(Ec, B, H, L, Mn, eqcurv, LongSteelRatio, bending, fpc, dcore, TransvSteelRatioY, Es, mombilin, curv, Force, mom, displf):
    
    G = 0.43 * Ec
    Agross = H * B
    As = (5 / 6) * Agross
    Ig = B * (H ** 3) / 12
    Ieff = (Mn * 1000 / (Ec * 1e6 * eqcurv)) * 1e12
    beta = min(0.5 + 20 * LongSteelRatio, 1)
    
    if bending.lower() == 'single':
        alpha = min(max(1, 3 - L / H), 1.5)
    elif bending.lower() == 'double':
        alpha = min(max(1, 3 - L / (2 * H)), 1.5)
    else:
        raise ValueError("Bending should be either 'single' or 'double'")
    
    Vc1 = 0.29 * alpha * beta * 0.8 * np.sqrt(fpc) * Agross / 1000
    
    kscr = (0.25 * TransvSteelRatioY * Es * (B / 1000) * ((H - dcore) / 1000) / (0.25 + 10 * TransvSteelRatioY)) * 1000
    forcebilin = 0
    if bending.lower() == 'single':
        ksg = (G * As / L) / 1000
        kscr = kscr / L
        forcebilin = np.array(mombilin) / (L / 1000)
    else:  # bending == 'double'
        ksg = (G * As / (L / 2)) / 1000
        kscr = kscr / (L / 2)
        forcebilin = 2 * np.array(mombilin) / (L / 1000)
    
    kseff = ksg * (Ieff / Ig)
    #return f"ksg={ksg} - Ieff={Ieff} - Ig={Ig}"
    aux = (Vc1 / kseff) / 1000
    aux2 = 0
    momaux = mom.copy()
    displsh = []
    np.savetxt("curv.txt", curv, delimiter=",", fmt="%.6f")
    momaux = np.array(momaux, dtype=float) 
    for i in range(len(curv)):
        if momaux[i] <= Mn and Force[i] < Vc1:
            displsh.append((Force[i] / kseff) / 1000)

        elif momaux[i] <= Mn and Force[i] >= Vc1:
            displsh.append(((Force[i] - Vc1) / kscr) / 1000 + aux)

        else:  # momaux_copy[i] > Mn
            momaux *= 4  
            aux3 = i - aux2
            aux2 += 1
            displsh.append((displf[i] / displf[i - 1]) * displsh[i - 1]) 
            
    return Vc1, kscr, kseff, forcebilin, displsh, G
    
    
#####################################################################  shear_strength
def shear_strength(ncy, Dh, fyh, H, clb, cMn, s, theta, thetaD, LongSteelRatio, displ, dyf, L, P, bending, ductilitymode, fpc, Agross, phiS, Force, Dduct, mom, curv, CuDu, dy):
    cMn=100.8640809281518
    Vs = (ncy * 0.25 * np.pi * (Dh**2) * fyh * (H - clb + 0.5 * Dh - cMn) * (1 / np.tan(np.pi / 6)) / s) / 1000
    Vsd = (ncy * 0.25 * np.pi * (Dh**2) * fyh * (H - clb + 0.5 * Dh - cMn) * (1 / np.tan((35 / 180) * np.pi)) / s) / 1000
    beta = min(0.5 + 20 * LongSteelRatio, 1)
    Dductf = displ / dyf
    #return f"Vs={Vs} - Vsa={Vsd} - beta={beta}"

    if bending.lower() == 'single':
        alpha = min(max(1, 3 - L / H), 1.5)
        Vp = (P * (H - cMn) / (2 * L)) / 1000 if P > 0 else 0
    elif bending.lower() == 'double':
        alpha = min(max(1, 3 - L / (2 * H)), 1.5)
        Vp = (P * (H - cMn) / L) / 1000 if P > 0 else 0

    Vc=0

    if ductilitymode.lower() == 'uniaxial':
        Vc = [alpha * beta * min(max(0.05, 0.37 - 0.04 * Dductf[i]), 0.29) * 0.8 * np.sqrt(fpc) * Agross / 1000 for i in range(len(Dductf))]
    elif ductilitymode.lower() == 'biaxial':
        Vc = [alpha * beta * min(max(0.05, 0.33 - 0.04 * Dductf[i]), 0.29) * 0.8 * np.sqrt(fpc) * Agross / 1000 for i in range(len(Dductf))]
    
    Vcd = 0.862 * np.array(Vc)
    Vpd = 0.85 * Vp
    V = phiS * (np.array(Vc) + Vs + Vp)
    Vd = 0.85 * (Vcd + Vsd + Vpd)
    criteria = 1
    faildispl = None
    failforce = None
    failduct = None 
    failcurv = None
    failCuDu = None
    failmom = None

    np.savetxt("V.txt", V, delimiter=",", fmt="%.8f")    
    np.savetxt("Vc.txt", Vc, delimiter=",", fmt="%.8f") 
    #return f"Vs={Vs} _____ Vp={Vp} _____ ncy={ncy} _____ Dh={Dh} _____ fyh={fyh} _____ H={H} _____ clb={clb} _____ Dh={Dh} _____ cMn={cMn} _____ s={s}"

    if V[-1] < Force[-1]:
        failure = V - Force
        np.savetxt("failure.txt", failure, delimiter=",", fmt="%.8f")
        faildispl = np.interp(0, failure, displ)
        failforce = np.interp(faildispl, displ, Force)
        failduct = np.interp(faildispl, displ, Dduct)
        failmom = np.interp(faildispl, displ, mom)
        failcurv = np.interp(faildispl, displ, curv)
        failCuDu = np.interp(faildispl, displ, CuDu)

        if bending.lower() == 'single':
            if faildispl <= 2 * dy:
                criteria = 2
            elif faildispl < 8 * dy:
                criteria = 3
            else:
                criteria = 4
        elif bending.lower() == 'double':
            if faildispl <= 1 * dy:
                criteria = 2
            elif faildispl < 7 * dy:
                criteria = 3
            else:
                criteria = 4
    
    #return f"{len(displ)} displ={[f'{x:.6f}' for x in displ]} ______ V={V} ________ Force={Force} ________ {len(failure)} failure={[f'{x:.6f}' for x in failure]}  ________  faildispl={faildispl}  ________   failforce={failforce}"
    return V, Vd, criteria, faildispl, failforce, failduct, failcurv, failCuDu, failmom
   

def defVals(ecu):
    def_vals = []    
    if ecu <= 0.0018:
        def_vals = np.arange(0.0001, 20*ecu, 0.0001)
    elif ecu > 0.0018 and ecu <= 0.0025:
        def_vals = np.concatenate([np.arange(0.0001, 0.0016, 0.0001), np.arange(0.0018, 20*ecu, 0.0002)])
    elif ecu > 0.0025 and ecu <= 0.006:
        def_vals = np.concatenate([np.arange(0.0001, 0.0016, 0.0001), np.arange(0.0018, 0.002, 0.0002), np.arange(0.0025, 20*ecu, 0.0005)])
    elif ecu > 0.006 and ecu <= 0.012:
        def_vals = np.concatenate([np.arange(0.0001, 0.0016, 0.0001), np.arange(0.0018, 0.002, 0.0002), np.arange(0.0025, 0.005, 0.0005), np.arange(0.006, 20*ecu, 0.001)])
    else:
        def_vals = np.concatenate([np.arange(0.0001, 0.0016, 0.0001), np.arange(0.0018, 0.002, 0.0002), np.arange(0.0025, 0.005, 0.0005), np.arange(0.006, 0.01, 0.001), np.arange(0.012, 20*ecu, 0.002)])
    return def_vals
    
def on_PP(PP, i, np_vals, def_vals, ecun, fcun, conclay, ec, fc, es, fs, Ast):
    defVals = def_vals 
    if PP[i] > 0:
        for j in range(np_vals):
            compch = np.sum(np.interp(def_vals[0], ecun, fcun) * conclay[:, 1]) + \
            np.sum(np.interp(def_vals[0], ec, fc) * conclay[:, 2]) + \
            np.sum(Ast * np.interp(def_vals[0], es, fs))
            if compch < PP[i]:
                defVals = def_vals[1:]
    return defVals
    
def calculate_mom_curv(tolerance, H, B, fpc, PP, i, x, np_vals, def_vals, conclay, distld, itermax, ecun, ecu, esu, fcun, ec, fc, es, fs, Ast, confined, unconfined, cores, curv, mom, ejen, DF, vniter, coverstrain, corestrain, steelstrain):
    eec = 0
    ees = 0
    tol = tolerance * H * B * fpc
    message = None
    TF = None

    for k in range(np_vals):
        F = 10 * tol
        niter = 0
        while abs(F) > tol:
            if x[niter] <= H:
                eec = (def_vals[k] / x[niter]) * (conclay[:, 0] - (H - x[niter]))
                ees = (def_vals[k] / x[niter]) * (distld - (H - x[niter]))
            else:
                eec = (def_vals[k] / x[niter]) * (x[niter] - H + conclay[:, 0])
                ees = (def_vals[k] / x[niter]) * (x[niter] - H + distld)  

            fcunconf = np.interp(eec, ecun, fcun)  
            fcconf = np.interp(eec, ec, fc)
            fsteel = np.interp(ees, es, fs)
                           
            FUNCON = fcunconf * conclay[:, 1]
            FCONF = fcconf * conclay[:, 2]
            FST = Ast * fsteel
            F = np.sum(FUNCON) + np.sum(FCONF) + np.sum(FST) - PP[i]                

            if F > 0:
                x = np.append(x, x[niter] - 0.05 * x[niter])
            if F < 0:
                x = np.append(x, x[niter] + 0.05 * x[niter])          

            niter += 1
            if niter > itermax:
                message = 3
                break
            
        cores = (def_vals[k] / x[niter]) * abs(x[niter] - H)
        TF = (confined == unconfined)
        if not TF:
            if cores >= ecu:
                message = 1
                break
        if TF:
            if def_vals[k] >= ecu:
                message = 1
                break
            
        if abs(ees[0]) > esu:
            message = 2
            break   
        
        ejen.append( x[niter] )
        DF.append( F )
        vniter.append( niter )
        mom.append( (np.sum(FUNCON * conclay[:, 0]) + np.sum(FCONF * conclay[:, 0]) + np.sum(FST * distld) - PP[i] * (H / 2)) / (10 ** 6) )
        curv.append( 1000 * def_vals[k] / x[niter] )
        coverstrain.append( def_vals[k] )
        corestrain.append( cores )
        steelstrain.append( ees[0] )
        x[0] = x[niter]
        x[1:] = 0
        if message != 0:
            break            

    return message, cores, curv, mom, ejen, DF, vniter, coverstrain, corestrain, steelstrain

#####################################################################  calculate_moment_capacity
def calculate_moment_capacity(PTid, fpc, Agross, PCid, ecun, ecu, esu, fcun, conclay, distld, ec, es, fc, fs, Ast, H, B, tolerance, itermax, csid, ssid, confined, unconfined, cores, esaux):
                                                                                 
    # Prepare the initial PP array
    PP = np.concatenate([np.arange(-0.90*PTid, 0, 0.30*PTid), 
                         np.arange(0.05*fpc*Agross, 0.7*PCid, 0.05*fpc*Agross)])
    
    Mni = np.zeros(len(PP))
    eci = np.zeros(len(PP))
    esi = np.zeros(len(PP))
    mess = np.zeros(len(PP))
    def_vals = []
    for i in range(len(PP)):
        curv=[];mom=[];ejen=[];DF=[];vniter=[];coverstrain=[];corestrain=[];steelstrain=[]
        curv.append(0) # curvatures
        
        mom.append(0)   # moments
        ejen.append(0)  # neutral axis
        DF.append(0)    # force equilibrium
        
        vniter.append(0)  # iterations
        coverstrain.append(0)
        corestrain.append(0)

        steelstrain.append(0)

        x = np.zeros(1)
        message = 0

        ####@ Define def_vals based on ecu
        def_vals = defVals(ecu)
        
        np_vals = len(def_vals)

        ####@ Adjust def_vals based on PP[i]
        def_vals = on_PP(PP, i, np_vals, def_vals, ecun, fcun, conclay, ec, fc, es, fs, Ast)        

        np_vals = len(def_vals)
        """
        curv[0] = 0 
        mom[0,0] = 0
        ejen[0,0] = 0
        DF[0,0] = 0
        vniter[0,0] = 0
        coverstrain[0,0] = 0
        corestrain[0,0] = 0
        steelstrain[0,0] = 0 
        """
        x[0] = H / 2


        # Main loop to calculate moments and curvatures

        message, cores, curv, mom, ejen, DF, vniter, coverstrain, corestrain, steelstrain = calculate_mom_curv(tolerance, H, B, fpc, PP, i, x, np_vals, def_vals, conclay, distld, itermax, ecun, ecu, esu, fcun, ec, fc, es, fs, Ast, confined, unconfined, cores, curv, mom, ejen, DF, vniter, coverstrain, corestrain, steelstrain)

        ssid = float(ssid)
        Mni[i] = np.interp(csid, coverstrain, mom)
        esaux = np.interp(csid, coverstrain, steelstrain)
        
        Mnsteel = np.interp(-ssid, steelstrain, mom)

        cr = 0  # concrete control
        if abs(esaux) > abs(ssid) or np.isnan(Mni[i]):
            Mnsteel = np.interp(-ssid, steelstrain, mom)
            if not np.isnan(Mnsteel):
               cr = 1  # steel control
               Mni[i] = Mnsteel
            elif np.isnan(Mni[i]) and np.isnan(Mnsteel):
                print('problem with strain values to estimate PM interaction at axial force:', PP[i])

        if cr == 0:
            eci[i] = csid
            esi[i] = esaux
        if cr == 1:
            esi[i] = -ssid
            eci[i] = np.interp(-ssid, steelstrain, coverstrain)
        mess[i] = message

    Mni = np.concatenate([[0], Mni, [0]])
    PPn = np.concatenate([[-PTid], PP, [PCid]])

    MB = np.max(Mni)
    PB = PPn[np.argmax(Mni)]

    PB13 = (1 / 3) * PB
    MB13 = np.interp(PB13, PPn, Mni)

    PB23 = (2 / 3) * PB
    MB23 = np.interp(PB23, PPn, Mni)

    MB0 = np.interp(0, PPn, Mni)

    PPL = np.array([-PTid, 0, PB13, PB23, PB, PCid])
    MnL = np.array([0, MB0, MB13, MB23, MB, 0])

    outputi = np.vstack([Mni, PPn / 1000])

    
    return outputi, esi, eci, i, csid, ssid, outputi, PTid, PCid, PB, MB, PB13, MB13, PB23, MB23, PPn, MnL, PPL, Mni