#CONICAL BAR LOADED BY A MOMENT ABOUT THE x-AXIS
#
#Stresses vary with R^{-3} so stress functions must vary with
#R^{-1} (A and E) and R^{-2} (B). 
#
#The required function for B is R^{-2}P_1^1(cos(beta))sin(theta), 
#but this form would give m>n for A,E, so we need to use the log(R+z)
#form of equations (22.37).
#
#To get the appropriate non-axisymmetric function, take log(R+z)
#and differentiate with respect to y, obtaining 
#
phi:=unapply(A*sin(beta)*sin(theta)/(R*(1+cos(beta))),R,theta,beta);
omega:=unapply(B*sin(beta)*sin(theta)/R^2,R,theta,beta);
psi:=unapply(E*sin(beta)*cos(theta)/(R*(1+cos(beta))),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);
stb1:=unapply(stb,R,beta);
sRt1:=unapply(sRt,R,beta);
#
c1:=simplify(sbR1(R,beta0)/sin(theta));
c2:=simplify(sbb1(R,beta0)/sin(theta));
c3:=simplify(stb1(R,beta0)/cos(theta));
#
t:=unapply(simplify(sbR1(R,beta)*sin(theta)
+sRt1(R,beta)*cos(theta)*cos(beta)),theta);
c4:=unapply(int(t(theta),theta=0..2*Pi),beta);
c5:=simplify(int(c4(beta)*R^3*sin(beta),beta=0..beta0));
#
#
solution:=simplify(solve({c1=0,c2=0,c3=0,c5=M},{A,B,E}));
#
sRR2:=simplify(subs(solution, sRR1(R,beta)));
stt2:=simplify(subs(solution, stt1(R,beta)));
sbb2:=simplify(subs(solution, sbb1(R,beta)));
stb2:=simplify(subs(solution, stb1(R,beta)));
sRt2:=simplify(subs(solution, sRt1(R,beta)));