import math
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate 

T=input("What is the type of beam :  ")
if(T=="Simply supported beam" or T=="simply supported beam" or T=="ssb"):
    b=int(input("Width of wall in mm ="))
    l=int(input("Clear span in mm ="))
    imposedload = int(input("Imposed load in N/mm="))
    Fck=int(input("Grade of concrete in N/mm2="))
    Fy=int(input("Grade of steel in N/mm2 ="))
    
    print("DIMENSIONING OF SECTION:-")
    print("b = width of wall =",b,"mm")
    d=l/20
    print("depth d =",d,"mm")
    provideddepth=int(input("Provided depth="))
    effectivedepth=50
    D=provideddepth+effectivedepth
    print("Total Depth=",D,"mm")
    
    print("EFFECTIVE SPAN:-")
    Leff1=l+d
    Leff2=l+b
    Leff=min(Leff1,Leff2)
    print("Effective span=",Leff,"mm")
        
    print("COMPUTE THE LOADS:-")
    unitweight=25
    Dl=b*D*unitweight/1000
    print("Dead load=",Dl/1000,"KN/m")
    LL=3000
    print("Imposed load=",imposedload/1000,"KN/m")
    W=Dl+LL+imposedload
    print("Total load=",W/1000,"KN/m")
    partialfactor=1.5
    Wu=W*partialfactor
    print("Factored load=",Wu/1000,"KN/m")
    
    print("MOMENT AND SHEAR:-")
    Mu=Wu * Leff *Leff /8
    print("Momemt=",Mu*(10**-6)/1000,"KN-m")
    Vu=Wu * Leff /2
    print("Shear=",Vu*(10**-3)/1000,"KN")
    
    print("CHECK FOR EFFECTIVE DEPTH:-")
    if(Fy == 450):
        req=Mu/(1000*0.138*Fck*b)
    elif(Fy == 500):
        req=Mu/(1000*0.1336*Fck*b)
    Dreq= math.sqrt(req)
    print(" Required depth=",Dreq,"mm")
    if(provideddepth>Dreq):
        print("Provided depth is safe")
    else:
        print("change the depth and do it again or go for doubled reinforced beams")
        
    print("TENSION REINFORCEMENT:-")
    Mulim=(0.1336*Fck*b*provideddepth*provideddepth)/1000000
    print("Mu,lim=",Mulim,"KNm")
    if(Mu<Mulim):
        print("Mu<Mu,lim  , section is under reinforced")
    else:
        print("section is under reinforced")
    A=(0.87*Fy*Fy)/(Fck*b) 
    B=-0.87*Fy*provideddepth
    C=Mu/1000
    Ast1=(-B+(math.sqrt((B*B)-(4*A*C))))/(2*A)
    Ast2=(-B-(math.sqrt((B*B)-(4*A*C))))/(2*A)
    print(Ast1,Ast2)
    Ast=min(Ast1,Ast2)
    print("Ast",Ast,"mm^2")
    d1=20
    ast=(3.14*d1**2/4)
    N0=(Ast//ast)+1
    Astpro=(ast*N0)
    print("provided tension reinforcement=",Astpro,"mm^2")
    print("provide 20mm bars for reinfocement")
    print("DESIGN OF SHEAR:-")
    Tv = Vu/(b*provideddepth*1000)
    print("Nominal shear stress Tv = ",Tv,"N/mm^2")
    Pt=float((100*Astpro)/(b*provideddepth))
    if Fck==20:
        ptp=[0.15,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3]
        tauc=[0.28,0.36,0.48,0.56,0.62,0.67,0.72,0.75,0.79,0.81,0.82,0.82,0.82]
    elif Fck==25:
        ptp=[0.15,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3]
        tauc=[0.29,0.36,0.49,0.57,0.64,0.7,0.74,0.78,0.82,0.85,0.88,0.9,0.92]
    else:
        ptp=[0.15,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3]
        tauc=[0.29,0.37,0.5,0.59,0.66,0.71,0.76,0.8,0.84,0.88,0.91,0.94,0.96]
    print("Provide Pt % is", round(Pt,2))
    Tc= scipy.interpolate.interp1d(ptp,tauc)
    print(Tc(Pt))

    
    Vus=(Vu/1000)-(Tc(Pt) *b*provideddepth)
    print(Vus)
    ds=16
    print("Assume providing 2 legged",ds,"mm diameter bars are used as stirrups")
    Asv=(2*3.14*ds**2)/4
    print("Asv is the total c/s area of stirrups, Asv=",Asv)
    print("For vertical stirrups,   \n Vus = (0.87*Fy*Asv*d)/Sv ")
    Sv=(0.87*Fy*Asv*provideddepth)/Vus
    print("Sv = ",Sv)
    print("provide 2 legged",ds,"mm diameter stirrups at maximum spacing")
    print("As per cl 26.5.1 pg.47")
    spacing=(min((0.75*provideddepth),300))
    print("provide maximum spacing of shear reinforcement =",spacing)
    
    print(" CHECK FOR DEFLECTION")
    print("(l/d) = 20           \n      As per Cl.23.2.1.(C) Pg.37      \n (l/d)max = (l/d) x Modification Factor       \n      Modification factor obtained from fig 4, pg.38 ")
    Fs=0.58*Fy * (Ast/Astpro)
    print("Fs=",Fs,"N/mm^2")
    if Fs<=130:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,2,1.8,1.5,1.37,1.21]
    elif Fs>130 and Fs <=160:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,2,1.58,1.35,1.2,1.1]
    elif Fs>160 and Fs<=220:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,1.65,1.3,1.1,1,0.95]
    elif Fs>220 and Fs<=260:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,1.35,1.03,0.97,0.85,0.8]
    else:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,1.1,0.9,0.8,0.75,0.7]
    kt=scipy.interpolate.interp1d(ptaa,ktt)
    Kt=kt(Pt)
    Kt1= Kt * 20 
    if(Kt1  >  (l/provideddepth)):
        print("section is safe in deflection")
    else:
        print("section is not safe in deflection")
        
else:
    b=int(input("Width of wall in mm ="))
    l=int(input("Clear span in mm ="))
    imposedload = int(input("Imposed load in N/mm="))
    Fck=int(input("Grade of concrete in N/mm2="))
    Fy=int(input("Grade of steel in N/mm2 ="))
    
    print("DIMENSIONING OF SECTION:-")
    print("b = width of wall =",b,"mm")
    d=l/7
    print("depth d =",d,"mm")
    provideddepth=int(input("provided depth d in mm="))
    effectivedepth=50
    D=provideddepth+effectivedepth
    print("Total Depth=",D,"mm")
    
    print("EFFECTIVE SPAN:-")
    Leff1=l+(0.5*provideddepth)
    Leff2=l+b
    Leff=min(Leff1,Leff2)
    print("Effective span=",Leff,"mm")
        
    print("COMPUTE THE LOADS:-")
    unitweight=25
    Dl=b*D*unitweight/1000
    print("Dead load=",Dl/1000,"KN/m")
    LL=3000
    print("Imposed load=",imposedload/1000,"KN/m")
    W=Dl+LL+imposedload
    print("Total load=",W/1000,"KN/m")
    partialfactor=1.5
    Wu=W*partialfactor
    print("Factored load=",Wu/1000,"KN/m")
    
    print("MOMENT AND SHEAR:-")
    Mu=Wu * Leff *Leff /2
    print("Momemt=",Mu*(10**-6)/1000,"KN-m")
    Vu=Wu * Leff 
    print("Shear=",Vu*(10**-3)/1000,"KN")
    
    print("CHECK FOR EFFECTIVE DEPTH:-")
    if(Fy==450):
         req=Mu/(1000*0.1336*Fck*b)
    elif(Fy==500):  
        req=Mu/(1000*0.1336*Fck*b)
    Dreq= math.sqrt(req)
    print(" Required depth=",Dreq,"mm")
    if(provideddepth>Dreq):
        print("Provided depth is safe")
    else:
        print("change the depth and do it again")
        
    print("TENSION REINFORCEMENT:-")
    Mulim=(0.1336*Fck*b*provideddepth*provideddepth)/1000
    print("Mu,lim=",Mulim*(10**-3),"KN-m")
    if(Mu<Mulim):
        print("Mu<Mu,lim  , section is under reinforced")
    else:
        print("section is over reinforced")
    A=(0.87*Fy*Fy)/(Fck*b) 
    B=-0.87*Fy*provideddepth
    C=Mu/1000
    Ast1=(-B+(math.sqrt((B*B)-(4*A*C))))/(2*A)
    Ast2=(-B-(math.sqrt((B*B)-(4*A*C))))/(2*A)
    print(Ast1,Ast2)
    Ast=min(Ast1,Ast2)
    print("Ast",Ast,"mm^2")
    d1=20
    ast=(3.14*d1**2/4)
    N0=(Ast/ast)
    Astpro=(ast*N0)
    print("provided tension reinforcement=",Astpro,"mm^2")
    
    print("DESIGN OF SHEAR:-")
    Tv = Vu/(b*provideddepth*1000)
    print("Nominal shear stress Tv = ",Tv,"N/mm^2")
    Pt=float((100*Astpro)/(b*provideddepth))
    print(Pt)

    if Fck==20:
        ptp=[0.15,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3]
        tauc=[0.28,0.36,0.48,0.56,0.62,0.67,0.72,0.75,0.79,0.81,0.82,0.82,0.82]
    elif Fck==25:
        ptp=[0.15,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3]
        tauc=[0.29,0.36,0.49,0.57,0.64,0.7,0.74,0.78,0.82,0.85,0.88,0.9,0.92]
    else:
        ptp=[0.15,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3]
        tauc=[0.29,0.37,0.5,0.59,0.66,0.71,0.76,0.8,0.84,0.88,0.91,0.94,0.96]
    print("Provide Pt % is", round(Pt,2))
    Tc= scipy.interpolate.interp1d(ptp,tauc)
    print(Tc(Pt))

    Vus=(Vu/1000)-(Tc(Pt) *b*provideddepth)
    print(Vus)
    ds=16
    print("Assume providing 2 legged",(ds/2),"mm diameter bars are used as stirrups")
    Asv=(2*3.14*ds**2)/4
    print("Asv is the total c/s area of stirrups, Asv=",Asv)
    print("For vertical stirrups,   \n Vus = (0.87*Fy*Asv*d)/Sv ")
    Sv=(0.87*Fy*Asv*provideddepth)/Vus
    print("Sv = ",Sv)
    print("provide 2 legged",ds/2,"mm diameter stirrups at maximum spacing")
    print("As per cl 26.5.1 pg.47")
    spacing=(min((0.75*provideddepth),300))
    print("provide maximum spacing of shear reinforcement =",spacing)
        
    print(" CHECK FOR DEFLECTION")
    Fs=0.58*Fy * (Ast/Astpro)
    print("Fs=",Fs,"N/mm^2")
    if Fs<=130:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,2,1.8,1.5,1.37,1.21]
    elif Fs>130 and Fs <=160:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,2,1.58,1.35,1.2,1.1]
    elif Fs>160 and Fs<=220:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,1.65,1.3,1.1,1,0.95]
    elif Fs>220 and Fs<=260:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,1.35,1.03,0.97,0.85,0.8]
    else:
        ptaa=[0,0.4,0.8,1.2,1.6,2.0]
        ktt=[2,1.1,0.9,0.8,0.75,0.7]
    kt=scipy.interpolate.interp1d(ptaa,ktt)
    Kt=kt(Pt)
    Kt1= Kt * 20 
    if(Kt1  >  (l/provideddepth)):
        print("section is safe in deflection")
    else:
        print("section is not safe in deflection")



x=(0,0,b,b,0)
y=(0,D,D,0,0)
x1=(effectivedepth/2,effectivedepth/2,b-(effectivedepth/2),b-(effectivedepth/2),effectivedepth/2)
y1=(effectivedepth/2,D-(effectivedepth/2),D-(effectivedepth/2),effectivedepth/2,effectivedepth/2)
plt.plot(x,y)
plt.plot(x1,y1)
plt.scatter(effectivedepth,(effectivedepth/2 )+5,s=ds)
plt.scatter(effectivedepth,D-(effectivedepth/2 )-5,s=ds)
plt.scatter(b-(effectivedepth),D-(effectivedepth/2 )-5,s=ds)
plt.scatter(b-(effectivedepth),(effectivedepth/2)+5,s=ds)
x2=(b+100,b+100,l+b+100,l+b+100,b+100)
x3=(b+100,l+b+100)
y3=(effectivedepth/2,effectivedepth/2)
y2=(0,D,D,0,0)
x4=(b+100,l+b+100)
y4=(D-(effectivedepth/2),D-(effectivedepth/2))
plt.plot(x2,y2)
plt.plot(x3,y3)
plt.plot(x4,y4)
x5=(1600,1600,2800,2800,4000,4000,5200,5200,6400,6400)
y5=(25,525,525,25,25,525,525,25,25,525)
x6=(1600,2800,4000,5000)
y6=(300,300,300,600)
plt.plot(x5,y5)
plt.plot(x6,y6,marker='x')
plt.title("DETAILING OF BEAMS")
plt.show()
