b=110;Ef=144000;efu=0.0181;ffu=Ef*efu;b1=0.8;ecu=0.0033;Af=142.5;d=195;
Efv=107000;efvu=0.00668;ffvu=Efv*efvu;cp=50;
t=0.000001;%the shear stress at the calculation point
fyx=efu*Ef;fyy=efvu*Efv;fc=-28;smx=180-2.7116*(600-cp)^2/10000;smy=150;fcr=0.33*(-fc)^0.5;ec=-0.0035;Ec=22*(abs(fc)/10)^0.33*1000;ecr=fcr/Ec;agg=10;

psy=0.0025; %%%%%%%%%%%%%%%%%%%SELF DEFINED
e1=0.00000001;
c1=2;
n=0;c2=1;%concrete crush parameters
c3=10;
tana=2*2.7116*(600-cp)/10000;%the angle of inclined flexural reinforcement
at=atan(tana);

du=195-2.7116*(600-400)^2/10000;

syms x ef Mu
[Mu,ef,x]=solve(b1*(-fc)*b*x==Af*Ef*ef, b1*(-fc)*b*x*(du-x/2)==Mu, x*ef==(du-x)*ecu);
ef=ef(ef>0);
x=x(x>0);
Mu=Mu(Mu>0);
x=eval(x);
Mu=eval(Mu);
ef=eval(ef);


P=Mu/400*3/2;%shear force according to max bending moment
Vm=P*2/3;%assign the shear force as Vm
V=Vm;

while abs(c3)>5,
    
Mc=V*cp;%The calculation section's moment 150 is the distance away from support
dc=195-2.7116*(600-cp)^2/10000;%The effective depth of the calculation section which is 150mm from the support;
%calculate the concrete strain/the tensile force/the compression force 
psx=Af/b/dc;
syms xc efc ecc

[ecc,efc,xc]=solve(b1*Ec/2*(-ecc)*b*xc==Af*Ef*efc, b1*Ec/2*(-ecc)*b*xc*(dc-xc/2*b1)==Mc, xc*efc==(dc-xc)*(-ecc));
efc=efc(efc>0);
xc=xc(xc>0);
ecc=ecc(ecc<0);
xc=eval(xc);
ecc=eval(ecc);
efc=eval(efc);


Cc=b1*Ec/2*(-ecc)*b*xc;%the compression force of the compression area
efc=Cc/Ef/Af;
Vf=1;




Tl=efc*Af*Ef;
Tv=Tl*tana;%vertical component of tensile force of flexural reinforcement
Ftd=1.3*(V-0.5*Vf-Tv);
ex=efc+Ftd/Ef/Af;

while abs(c1)>0.1,
if e1>ecr,
    fc1=fcr/(1+(200*e1)^0.5);
else
    fc1=Ec*e1;
end
A=-Efv*psy*(e1-ex);
B=-t;
C=Efv*e1*psy+fc1;
a=((-B-(B^2-4*A*C)^0.5)/2/A);
theta=atan(a);
rxy=2*(e1-ex)*tan(theta);
ey=e1-(e1-ex)*tan(theta)^2;
e2=ex-(e1-ex)*tan(theta)^2;
if -e2+ec>0,
    c2=0;
    break
end
fcx=fc1-t/tan(theta);
fcy=fc1-t*tan(theta);
fc2=fc1-t*(tan(theta)+1/tan(theta));

if 1/(0.8-0.34*e1/ec)<1,
    fc2max=fc*(1/(0.8-0.34*e1/ec));
else
    fc2max=fc;
end

fc2a=fc2max*(2*(e2/ec)-(e2/ec)^2);

c1=fc2-fc2a;
if abs(c1)>0.1,
    e1=e1+0.000001;
end
n=n+1;
end

%judge if the reinforcement yield or what


if ey<efvu,
    fsy=ey*Efv;
else
    c2=0;
   break
    
end
if ex<efu,
    fsx=ex*Ef;
else
    c2=0;
   break
    
end

dfc1=fc1-psy*(fyy-fsy);
s=1/(sin(theta)/smx+cos(theta)/smy);
w=e1*s;
vcimax=(-fc)^0.5/(0.31+24*w/(agg+16));
if dfc1<=0,
    vci=0;
    fci=0;
else
    F=dfc1/tan(theta)-0.18*vcimax;
    if F<=0,
        fci=0;
        vci=dfc1/tan(theta);
    else
        D=0.82/vcimax;
        E=1/tan(theta)-1.64;
        fci=(-E-(E^2-4*D*F)^0.5)/2/D;
        vci=(fci+dfc1)/tan(theta);
    end
end

fsycr=fsy+(fc1+fci-vci*tan(theta))/psy;
fsxcr=fsx+(fc1+fci+vci/tan(theta))/psx;

%check failure mode
kc=1.64-1/tan(theta);
if kc<0,
    kc=0;
end


fm1=vcimax*(0.18+0.3*kc^2)*tan(theta)+psy*(fyy-fsy)-fc1;
fm2=(-fc2max+fc2);
fm3=(fyx-fsxcr);
fm4=(fyy-fsycr);



while ffvu>fsycr,

t=t+0.01;


e1=0.000001;

c1=2;n=0;c2=1;

while abs(c1)>0.1,
if e1>ecr,
    fc1=fcr/(1+(200*e1)^0.5);
else
    fc1=Ec*e1;
end
A=-Efv*psy*(e1-ex);
B=-t;
C=Efv*e1*psy+fc1;
a=((-B-(B^2-4*A*C)^0.5)/2/A);
theta=atan(a);
rxy=2*(e1-ex)*tan(theta);
ey=e1-(e1-ex)*tan(theta)^2;
e2=ex-(e1-ex)*tan(theta)^2;
if -e2+ec>0,
    c2=0;
    break
end
fcx=fc1-t/tan(theta);
fcy=fc1-t*tan(theta);
fc2=fc1-t*(tan(theta)+1/tan(theta));

if 1/(0.8-0.34*e1/ec)<1,
    fc2max=fc*(1/(0.8-0.34*e1/ec));
else
    fc2max=fc;
end

fc2a=fc2max*(2*(e2/ec)-(e2/ec)^2);

c1=fc2-fc2a;
if abs(c1)>0.1,
    e1=e1+0.000001;
end
n=n+1;
end

%judge if the reinforcement yield or what




if ey<efvu,
    fsy=ey*Efv;
else
   break;
end
if ex<efu,
    fsx=ex*Ef;
else
    c2=0;
   break
    
end

dfc1=fc1-psy*(fyy-fsy);
s=1/(sin(theta)/smx+cos(theta)/smy);
w=e1*s;
vcimax=(-fc)^0.5/(0.31+24*w/(agg+16));
if dfc1<=0,
    vci=0;
    fci=0;
else
    F=dfc1/tan(theta)-0.18*vcimax;
    if F<=0,
        fci=0;
        vci=dfc1/tan(theta);
    else
        D=0.82/vcimax;
        E=1/tan(theta)-1.64;
        fci=(-E-(E^2-4*D*F)^0.5)/2/D;
        vci=(fci+dfc1)/tan(theta);
    end
end

fsycr=fsy+(fc1+fci-vci*tan(theta))/psy;
fsxcr=fsx+(fc1+fci+vci/tan(theta))/psx;

%check failure mode
kc=1.64-1/tan(theta);
if kc<0,
    kc=0;
end

fm1=vcimax*(0.18+0.3*kc^2)*tan(theta)+psy*(fyy-fsy)-fc1;
fm2=(-fc2max+fc2);
fm3=(fyx-fsxcr);
fm4=(fyy-fsycr);

end

%calculate the shear reinforcement capacity
Vf=t*b*xc/2+(dc-xc)*b*t;
Vv=Tv+1.3*(V-0.5*Vf-Tv)*sin(at);

c3=Vf+Vv-V;

if Vf+Vv<V,
    V=0.9*V;
else
    V=1.1*V;
end

end
 if c3<0,
    V=V/0.9;
else
    V=V/1.1;
end  
      
        
P=V*3/2;