%width and depth
b=110;dm=195;
%longitudinal bars
Ef=144000;efu=0.0181;ffu=Ef*efu;Af=142.5;
%concrete strength and factors
fc=28;
Ec=22*(fc/10)^0.33*1000;
b1=0.85-0.05*(fc-28)/7;%concrete factor
if b1<0.65,
    b1=0.65;
end
ecu=0.003;
%bending capacity calculation point
a=400;
du=195-2.7116*(600-a)^2/10000; %effective depth at the loading point

%calculate the nominal flexural strength (Mn) of critical section 
pfb=0.85*b1*fc/ffu*(Ef*ecu)/(Ef*ecu+ffu);
pf=Af/b/du;
if pf>pfb,
    ff=((Ef*ecu)^2/4+0.85*b1*fc/pf*Ef*ecu)^0.5-0.5*Ef*ecu;
    Mn=pf*ff*(1-0.59*pf*ff/fc)*b*du^2;
else
    Mn=Af*ffu*(du-b1/2*(ecu/(ecu+efu))*d);
end

%applied load calculation according to nominal flexural strength
Pu=Mn/a*1.5;

% material property of shear links 
Efv=107000;efvu=0.004;ffvu=Efv*efvu;
%shear reinforcement ratio%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pfv=0.00;

%critical section
av=50;
dv=195-2.7116*(600-av)^2/10000;
pf2=Af/b/dv;%reinforcement raito of critical cross section

%shear contribution of concrete
nf=Ef/Ec;
k=(2*pf2*nf+(pf2*nf)^2)^0.5-pf2*nf;
c=k*dv;
Vc=0.4*(fc)^0.5*b*c;
%shear contribution of shear reinforcement;
Vf=pfv*ffvu*b*dv;

%shear resistance of concrete and shear reinforcement
Vv=Vc+Vf;


%shear contribution of longitudinal bars
tana=2*2.7116*(600-av)/10000;
at=atan(tana);

Thf=1;%assume the flexural tensile force at the critical cross section
Vtf=tana*Thf;%calculate the vertical component
Mv=Thf*(dv-Thf/(0.85*fc*b)/2);%calculate the bending moment
Va=Mv/(av);%calculate the applied shear force when Th is assumed;
Ftd=1.3*(Va-0.5*Vf-Vtf);
Vt=tana*(Thf+Ftd*cos(at));%calculate the vertical component


Vu=(Vt+Vv);%calculate the total shear resistance;
C0=abs(Va-Vu);%the difference between applied load and shear resistance
%iteratation to find the answer.
while C0>10,
    Thf=Thf+1;
    Vtf=tana*Thf;%calculate the vertical component
Mv=Thf*(dv-Thf/(0.85*fc*b)/2);%calculate the bending moment
Va=Mv/(av);%calculate the applied shear force when Th is assumed;
Ftd=1.3*(Va-0.5*Vf-Vtf);
Vt=tana*(Thf+Ftd*cos(at));%calculate the vertical component


Vu=(Vt+Vv);%calculate the total shear resistance;
C0=abs(Va-Vu);%the difference between applied load and shear resistance
%iteratation to find the answer.
end

Pvu=Vu*1.5;