function pts=asimp(init,final,alphastart)
global m0 m1 m2 r0 r1 r2;
%
% Using Simpson's rule this function integrates the constraint
% arising from conservation of angular momentum ti calculate
% the chane in pahse corresponding to moving along a line-segment
% in shape space.
% The input are the initial and terminal shapes, and the initial
% value of the phase (attitude alpha). The output is an array of
% all the points in the total space that were computed.
%
%
% the number of integration steps is controlled by the change
% in shape space:
%
dist=norm(final-init);
nn=2*(1+floor(dist*10));
hstep=dist/nn;
%
% Subdivide the line segment in shape space
% Temporarily use the third row to create the subdivision.
%
pts(1:3,:)=[diag(init)*ones(2,nn+1);(0:nn)/nn];
pts(1:2,:)=pts(1:2,:)+(final-init)*pts(3,:);
%
% Use Simpson's rule to integrate the change of phase. Note that both the
% numerator and the first term in the denominator are independent of the
% current shape (theta1,theta2). The first commands, in parallel, calculate
% the derivative of alpha at the partition points. The for-loop then sums
% the increments, saving the intermediate alpha values.
%
pts(3,1)=alphastart;
num=4*[m1*r1^2,m2*r2^2]*[final-init];
denom1=(m0+3*m1+3*m2)*r0^2+3*m1*r1^2+3*m2*r2^2;
denom2=6*r0*[m1*r1,m2*r2]*cos([pts(1:2,:)]);
increm=-num./(denom1+denom2)*hstep;
for i=1:nn/2
% Note that the first column contains the zeroth-point. The third column
% column contains the firt new point. The second column contains the inter-
% mediate values used in each step of Simpson's rule.
pts(3,2*i+1)=pts(3,2*i-1)+(increm(2*i-1)+4*increm(2*i)+increm(2*i+1))/6;
end;
% Finally discard half the points that were only needed for
% Simpson's rule's oversampling.
pts=pts(:,1:2:nn);