%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;
a1=0.85-0.0015*fc;
if a1<0.67,
    a1=0.67;
end
b1=0.97-0.0025*fc;%concrete factor
if b1<0.67,
    b1=0.67;
end
ecu=0.0035;
%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 
pf=Af/b/du;
ff=((Ef*ecu)^2/4+a1*b1*fc/pf*Ef*ecu)^0.5-0.5*Ef*ecu;
c=Af*ff/(a1*fc*b1*b);
Mn=Af*ff*(du-b1*c/2);


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

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

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

%shear contribution of concrete
k3=(dv/0.9/a)^0.5;
if k3>1,
    k3=1;
end

k4=1+(Ef*pf2)^(1/3);
Vc=0.05*k3*k4*(fc)^(1/3)*b*dv;

if Vc>0.22*(fc)^0.5*b*dv,
    Vc=0.22*(fc)^0.5*b*dv;
else
    if Vc<0.11*(fc)^0.5*b*dv,
        Vc=0.11*(fc)^0.5*b*dv;
    end
end
%%arch effect
ka=2.5*(dv/0.9/a);
if ka<1,
    ka=1;
else
    if ka>2.5,
        ka=2.5;
    end  
end
Vc=Vc*ka;

%member size
ks=750/(450+dv/0.9);
if ks<1,
    Vc=Vc*ks;
end

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

Thf=1;%assume the flexural tensile force at the critical cross section
Vtf=tana*Thf;%calculate the vertical component
Mv=Thf*dv;%calculate the bending moment

Va=Mv/(av);%calculate the applied shear force when Th is assumed;

Vfa=1;
Ftd=1.3*(Va-0.5*Vfa-Vtf);% additional tensile force
Vt=tana*(Thf+Ftd*cos(at));%calculate the vertical component
el=(Mv/dv+(Va-Vt))/2/(Ef*Af);
theta=30+7000*el;%angle of compressive stress in concrete
ko=cot(theta/180*pi);
Vf=pfv*ffvu*b*dv*ko;
while abs(Vf-Vfa)>1,
    Vfa=Vfa+1;
    Ftd=1.3*(Va-0.5*Vfa-Vtf);% additional tensile force
    Vt=tana*(Thf+Ftd*cos(at));%calculate the vertical component
    el=(Mv/dv+(Vf-Vt))/2/(Ef*Af);
    theta=30+7000*el;%angle of compressive stress in concrete
    ko=cot(theta/180*pi);
    Vf=pfv*ffvu*b*dv*ko;
end

Vu=(Vt+Vf+Vc);%calculate the total shear resistance;
C0=abs(Va-Vu);%the difference between applied load and shear resistance
%iteratation to find the answer.
while C0>100,
Thf=Thf+100;
Vtf=tana*Thf;%calculate the vertical component
Mv=Thf*(dv-Thf/(a1*fc*b)/2);%calculate the bending moment

Va=Mv/(av);%calculate the applied shear force when Th is assumed;

Vfa=1;
Ftd=1.3*(Va-0.5*Vfa-Vtf);% additional tensile force
Vt=tana*(Thf+Ftd*cos(at));%calculate the vertical component
el=(Mv/dv+(Vf-Vt))/2/(Ef*Af);
theta=30+7000*el;%angle of compressive stress in concrete
ko=cot(theta/180*pi);
Vf=pfv*ffvu*b*dv*ko;
while abs(Vf-Vfa)>1,
    Vfa=Vfa+1;
    Ftd=1.3*(Va-0.5*Vfa-Vtf);% additional tensile force
    Vt=tana*(Thf+Ftd*cos(at));%calculate the vertical component
    el=(Mv/dv+(Vf-Vt))/2/(Ef*Af);
    theta=30+7000*el;%angle of compressive stress in concrete
    ko=cot(theta/180*pi);
    Vf=pfv*ffvu*b*dv*ko;
end

Vu=(Vt+Vf+Vc);%calculate the total shear resistance;
C0=abs(Va-Vu);%the difference between applied load and shear resistance
end

Pvu=Vu*1.5;