#                        CONICAL BAR IN TENSION




phi:=unapply(A*ln(R*(1+cos(beta))),R,theta,beta);
psi:=0;
omega:=unapply(B/R,R,theta,beta);


read ABErtb;

sRR1:=unapply(sRR,R,beta);
stt1:=unapply(stt,R,beta);
sbb1:=unapply(sbb,R,beta);
sbR1:=unapply(sbR,R,beta);

c1:=sbR1(R,beta0);
c2:=sbb1(R,beta0);
f:=unapply(simplify(sRR1(Rb,beta)*cos(beta)-sbR1(Rb,beta)*sin(beta)),beta);
c3:=simplify(2*Pi*int(f(beta)*Rb^2*sin(beta),beta=0..beta0));

solution:=solve({c1=0,c2=0,c3=F},{A,B});

sRR2:=simplify(subs(solution,sRR1(R,beta)));
stt2:=simplify(subs(solution,stt1(R,beta)));
sbb2:=simplify(subs(solution,sbb1(R,beta)));
sbR2:=simplify(subs(solution,sbR1(R,beta)));

R:=sqrt(r^2+z^2);

phi:=unapply(A*ln(R+z),r,theta,z);
psi:=0;
omega:=unapply(B/R,r,theta,z);

read ABErtz;

r:=x*tan(beta0);

szz1:=subs(solution,szz);

szz2:=unapply(simplify(szz1*Pi*tan(beta0)^2/F),x,z,beta0,nu);

#This plot normalizes the axial stress distribution szz to the uniform
#value F/A.

plot({szz2(x,1,Pi/180,0.3),szz2(x,1,Pi/3,0.3),
szz2(x,1,Pi/12,0.3),szz2(x,1,Pi/6,0.3)},x=0..1,axes=boxed);