from utils import show
from math import sqrt, pi

from Designer import sst
SST = sst.SST()

Fy = 345.     # ASTM A992
E = 200000.
G = 77000.
phi = 0.9

Dsg = 'W250x101'
Lcol = 5500. # total length of column, in mm
Lx = Lcol       # unbraced length, x-axis
Ly = Lcol       # unbraced length, y-axis
Kx = 1.      # effective length factor
Ky = 1.
Lc = Ly      # unbraced length of compression flange
Cf = 1000.   # factored axial load, kN
Mfx1 = -200.  # moment about strong axis, top and bottom
Mfx2 = 240.  # same sign for double curvature
Mfy1 = 0.    # momemt about weak axis, top and bottom
Mfy2 = 0.
bracedx = True   # braced for bending about strong axis?
bracedy = True   # braced for bending about weak axis

A,D,T,B,W,Ix,Iy,Zx,Zy,Rx,Ry,J,Cw = SST.section(Dsg,'A,D,T,B,W,Ix,Iy,Zx,Zy,Rx,Ry,J,Cw')
show('Dsg,A,D,T,B,W,Rx,Ry,*1E6,Ix,Iy,*1E3,Zx,Zy,J,*1E9,Cw')

fclass = 4
if B/(2.*T) <= 145./sqrt(Fy):
    fclass = 1
elif B/(2.*T) <= 170./sqrt(Fy):
    fclass = 2
elif B/(2.*T) <= 200./sqrt(Fy):
    fclass = 3
    
H = (D - 2.*T)
phiCy = phi*A*Fy * 1E-3
wclass = 4
if (H/W) <= 1100.*(1-0.39*Cf/phiCy)/sqrt(Fy):
    wclass = 1
elif (H/W) <= 1700.*(1-0.61*Cf/phiCy)/sqrt(Fy):
    wclass = 2
elif (H/W) <= 1900.*(1-0.65*Cf/phiCy)/sqrt(Fy):
    wclass = 3
    
sclass = max(fclass,wclass)
show('B/(2.*T),145./sqrt(Fy),170./sqrt(Fy),200./sqrt(Fy),fclass',
     'H,H/W,Cf/phiCy',
     '1100.*(1-0.39*Cf/phiCy)/sqrt(Fy),1700.*(1-0.61*Cf/phiCy)/sqrt(Fy),1900.*(1-0.65*Cf/phiCy)/sqrt(Fy)'
     ,'wclass,sclass,sclass <= 2')
assert sclass <= 2, 'This only works for class 1 or class 2 sections'

Mfx = max(abs(Mfx1),abs(Mfx2))
Mfy = max(abs(Mfy1),abs(Mfy2))
show('Mfx,Mfy')

def kappa(Mf1,Mf2):    
    Mfmax = max(abs(Mf1),abs(Mf2))
    if Mfmax > 0.:
        kap = min(abs(Mf1),abs(Mf2))/Mfmax
        if Mf1*Mf2 > 0:
            return kap
        else:
            return -kap
    return 0.
def U1(L,I,Cf,Mf1,Mf2):
    w1 = max(0.6-0.4*kappa(Mf1,Mf2),0.4)
    Ce = pi*pi*E*I/(L*L) * 1E-3
    return w1,Ce,w1/(1.-Cf/Ce)

Mfx1,Mfx2,kappa(Mfx1,Mfx2),U1(Lx,Ix,Cf,Mfx1,Mfx2)

if bracedx or bracedy:
    Cr = phi*A*Fy * 1E-3
    Mrx = phi*Zx*Fy * 1E-6
    w1,Ce,U1xc = U1(Lx,Ix,Cf,Mfx1,Mfx2)
    U1x = max(1.0,U1xc)    
    lhs_a = Cf/Cr + 0.85*U1x*Mfx/Mrx
    show('Cf,Cr,w1,Ce,U1xc,U1x,Mfx,Mrx,Cf/Cr,0.85*U1x*Mfx/Mrx')
    if Mfy > 0.:
        Mry = phi*Zy*Fy * 1E-6
        Beta = 0.6
        w1y,Cey,U1yc = U1(Ly,Iy,Cf,Mfy1,Mfy2)
        U1y = max(1.0,U1yc)
        lhs_a += Beta*U1y*Mfy/Mry
        show('Beta,w1y,Cey,U1yc,U1y,Mfy,Mry,Beta*U1y*Mfy/Mry')
else:
    print('Unbraced frame: cross-sectional strength not checked')
    lhs_a = 0.
show('lhs_a,lhs_a <= 1')

assert Mfx > 0. or Mfy > 0., "No end bending moments in either direction"
K = 1.
if Mfx > 0. and Mfy == 0.:
    KLr = K*Lx/Rx
elif Mfx == 0. and Mfy > 0.:
    KLr = K*Ly/Ry
else:
    KLr = K*Ly/Ry
# Cr by 13.3
Fe = pi*pi*E/(KLr*KLr)
lamda = sqrt(Fy/Fe)
n = 1.34
Cr = phi*A*Fy*(1.+lamda**(2*n))**(-1./n) * 1E-3
Mrx = phi*Zx*Fy * 1E-6
U1x = U1(Lx,Ix,Cf,Mfx1,Mfx2)[2] if bracedx else 1.0
lhs_b = Cf/Cr + 0.85*U1x*Mfx/Mrx
show('Cf,K,K*Lx/Rx,K*Ly/Ry,KLr,Fe,lamda,n,Cr,U1x,Mfx,Mrx,Cf/Cr,0.85*U1x*Mfx/Mrx')
if Mfy > 0.:
    Mry = phi*Zy*Fy * 1E-6
    KLry = Ky*Ly/Ry
    Fey = pi*pi*E/(KLry*KLry)
    lamday = sqrt(Fy/Fey)
    Beta = min(0.85,0.6 + 0.4*lamday)
    U1y = U1(Ly,Iy,Cf,Mfy1,Mfy2)[2] if bracedy else 1.0
    lhs_b += Beta*U1y*Mfy/Mry
    show('lamday,Beta,U1y,Mfy,Mry,Beta*U1y*Mfy/Mry')
show('lhs_b,lhs_b <= 1')

KLr = max(Kx*Lx/Rx,Ky*Ly/Ry)
Fe = pi*pi*E/(KLr**2)
lamda = sqrt(Fy/Fe)
n = 1.34
Cr = phi*A*Fy*(1.+lamda**(2*n))**(-1./n) * 1E-3
show('Cf,Kx,Lx,Rx,Ky,Ly,Ry,Kx*Lx/Rx,Ky*Ly/Ry,KLr,Fe,lamda,n,Cr')
L = Lc
P = E*Iy*G*J
Q = Iy*Cw*(pi*E/L)**2
kap = kappa(Mfx1,Mfx2)
omega2 = min(1.75 + 1.05*kap + 0.3*kap**2, 2.5)
Mu = (omega2*pi/L) * sqrt(P+Q) * 1E-6
Mp = Zx*Fy * 1E-6
if Mu <= 0.67*Mp:
    Mrx = phi*Mu
else:
    Mrx = min(1.15*phi*Mp*(1-0.28*Mp/Mu),phi*Mp)
U1x = max(U1(Lx,Ix,Cf,Mfx1,Mfx2)[2],1.0) if bracedx else 1.0
show('L,kap,omega2,Mu,Mp,Mfx,U1x,Mrx,Cf/Cr,0.85*U1x*Mfx/Mrx')
lhs_c = Cf/Cr + 0.85*U1x*Mfx/Mrx
if Mfy > 0.:
    Mry = phi*Zy*Fy * 1E-6
    KLry = Ky*Ly/Ry
    Fey = pi*pi*E/(KLry*KLry)
    lamday = sqrt(Fy/Fey)
    Beta = min(0.85,0.6 + 0.4*lamday)
    U1y = U1(Ly,Iy,Cf,Mfy1,Mfy2)[2] if bracedy else 1.0
    lhs_c += Beta*U1y*Mfy/Mry
    show('Mfy,lamday,Beta,U1y,Mry,Beta*U1y*Mfy/Mry')
show('lhs_c,lhs_c <= 1')

lhs_d = Mfx/Mrx
if Mfy > 0.:
    lhs_d += Mfy/Mry
show('lhs_d, lhs_d < 1')

if lhs_a <= 1 and lhs_b <= 1 and lhs_c < 1 and lhs_d < 1:
    print(Dsg, 'is OK!')
else:
    print(Dsg, 'is NG!')