% %draw Simple Cylinder and shows how the trajectory intersects with the
% cylinder 
r=1;    % inner radius of the beaker
r1=2;   % outer radius of the beaker
i=1;
while i<=100
%getting random points
thetas=2*180*rand();
z1=rand();
R1=r+(r1-r)*(rand())^.5;  %creates random uniform distribution
x1=R1*cosd(thetas);
y1=R1*sind(thetas);

%getting random trajectory
R2=1.5;           % R2 is the length of the trajectory
theta= acosd(1-2*rand());      %uniform theta distribution in degrees. REMEBER this is in FLOATING POINT
Phi= 2*180*rand();             %uniform phi distribution in degrees. REMEBER this is in FLOATING POINT 
x2p= R2*sind(theta)*cosd(Phi);     % measured in the source coordinate
y2p= R2*sind(theta)*sind(Phi);     %measured in the source coordinate
z2p= R2*cosd(theta);             %measured in the source coordinate
x2=x2p+x1;                       % measured in the cylinder coordinate
y2=y2p+y1;                       % measured in the cylinder coordinate
z2=z2p+z1;                      % measured in the cylinder coordinate

%plotting
plot3([x1,x2],[y1,y2],[z1,z2]);
hold on
i=i+1;
end 

%plots a cylinder in 3d which can be rotated in 3d
cylinder(1);     % a cylinder with radius 1
colormap([1,1,1]);
hold on
daspect([1 1 1])    %streches the sphere
axis ([-1.6 1.6 -1.6 1.6])
xlabel('x-axis')
ylabel('y-axis')
zlabel('z-axis')
grid on
rotate3d on