%crank_problem.m
thetaindegrees=15*[0:23];
theta=(pi/12)*[0,  1,  2,  3,      4,      5,  6,     7,   8,     9,  10,   11,  12, 13, 14, 15, 16,     17,    18,    19, 20,    21,  22,     23]';
Poftheta=     [5,  5,  5,  5,  7.333,  9.666, 12,  12.5,  13,  13.5,  14, 14.5,  15, 15, 15, 15, 15,  13.25,  11.5,  9.75,  8,  7.25,  6.5,  5.75]';
#%figure(1)
#%plot(thetaindegrees,Poftheta,"-o")
a=0.090;
b=0.200;
acostheta=a*cos(theta);
asintheta=a*sin(theta);
x = acostheta + (b.^2 - asintheta.^2).^(1/2);
bcosphi=x-acostheta;
bsinphi=asintheta;
xposofpointB=acostheta;
yposofpointB=asintheta;
xposofpointC=x;

figure(2)
plot(thetaindegrees,xposofpointB,"-x",thetaindegrees,yposofpointB,"-o",thetaindegrees,xposofpointC,"-+")
xlabel("theta (degrees)")
title("Some plots: the x-component of point B, the y-component of point B, and the x-component of point C,\n all in meters, measured from an origin at point A, versus theta")
axis([0,345,-0.1,0.3])
line([90,90],[-0.1,0.3],"color","k")
line([180,180],[-0.1,0.3],"color","k")
line([270,270],[-0.1,0.3],"color","k")
legend("xposofpointB (m)","yposofpointB (m)", "xposofpointC (m)", "location", "east")
legend("boxon")

%torque: -M = tau = r x F = rx*Fy - ry*Fx
%where rx=a*cos(theta), ry=a*sin(theta),
%      Fx=-P(theta), Fy=P(theta)*sin(phi)/cos(phi)
Fx=-1.*Poftheta;
Fy=-1*Fx.*bsinphi./bcosphi;
tauoftheta = acostheta.*Fy - asintheta.*Fx;

figure(3)
plot(thetaindegrees,0.1*Poftheta,"-x",thetaindegrees,tauoftheta,"-o")
axis([0,345,-1.6,1.6])
line([0,345],[0,0],"color","k")
line([90,90],[-1.6,1.6],"color","k")
line([180,180],[-1.6,1.6],"color","k")
line([270,270],[-1.6,1.6],"color","k")
xlabel("theta (degrees)")
legend("0.1*P (N)","M (N*m)","location","northeast")
legend("boxon")
title("Crank moment M and piston force P, both ploted versus crank angle theta.\n  M is in newton*meters.  P is in 10s of newtons.")