gad.m

这是几何代数的matlab工具包

     disp('>> % 	INTERPOLATION OF ORIENTATIONS');
     % 	INTERPOLATION OF ORIENTATIONS
     global RA RB u v; %/
          clf; %/
     disp('>> %	Problem: interpolate two orientations.');
     %	Problem: interpolate two orientations.
     disp('>> %');
     %
     disp('>> % 	An orientation can be characterized ');
     % 	An orientation can be characterized 
     disp('>> %	by a rotation from a standard pose.');
     %	by a rotation from a standard pose.
     disp('>> %	Let the orientations be RA and RB.');
     %	Let the orientations be RA and RB.
     disp('>> %');
     %
     u = e1+e2-e3; %/
     v = e1+e3; %/
     view = [-1 2 -1 2 -2 1]; %/
     disp('>> % 	Initial orientation RA (applied to a bivector u^v):');
     % 	Initial orientation RA (applied to a bivector u^v):
     fprintf(1,'\n');     disp('>> RA = gexp(-I3*e1*pi/2/2);');
     RA = gexp(-I3*e1*pi/2/2);
     DrawBivector(RA*u/RA,RA*v/RA,'b');  %/
     GAtext(1.1* RA*(u+v)/RA,'" R_A"','b'); %/
     axis(view); axis off; %/
     GAview([30 30]); %/
     GAprompt; %/
     disp('>> % 	Final orientation (applied to u^v):');
     % 	Final orientation (applied to u^v):
     fprintf(1,'\n');     disp('>> RB = gexp(-I3*e2*pi/2/2);');
     RB = gexp(-I3*e2*pi/2/2);
     DrawBivector(RB*u/RB,RB*v/RB,'g');  axis(view); %/
     GAtext(1.1* RB*(u+v)/RB,'" R_B"','k'); %/
     GAprompt; %/
     disp('>> % 	Interpolation through division of total rotor:');
     % 	Interpolation through division of total rotor:
     fprintf(1,'\n');     fprintf(1,'>> Rtot =  RB/RA  ');
     input('');
     Rtot =  RB/RA %w
     disp('>> % which is done through incremental rotor R:');
     % which is done through incremental rotor R:
     fprintf(1,'\n');     n = 8;                          %/
     disp('>> R = gexp(sLog(Rtot)/n)');
     R = gexp(sLog(Rtot)/n)
     axisR = unit(GAZ(-sLog(R)/I3));   %/ 
     draw(axisR,'r'); %/
     title('R = exp( log( R_B/R_A ) / n)','Color','r'); %/
     draw(-sLog(Rtot)/n,'r'); %/
     axis(view); %/
     GAprompt; %/
     disp('>> % 	Execute the interpolations from RA[u^v] to RB[u^v]');
     % 	Execute the interpolations from RA[u^v] to RB[u^v]
          Ri = RA;  %/
          for i=1:n-1  %/
     	 disp(['i = ', num2str(i) ':   R' num2str(i) ' = R*R' num2str(i-1)]); %/
              Ri = R*Ri; %/
              DrawBivector(Ri*u/Ri,Ri*v/Ri); %/
              axis(view); %/
              drawnow; %/
          end %/
     GAprompt; %/
     title(''); %/
     GAorbiter(125,10); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 13 ) 
     GAps = '>>>> ';
     disp('>> %	HOMOGENEOUS REPRESENTATION');
     %	HOMOGENEOUS REPRESENTATION
     GAfigure; clc; %/
     disp('>> %	HOMOGENEOUS REPRESENTATION');
     %	HOMOGENEOUS REPRESENTATION
     disp('>> %');
     %
     disp('>> % 	For 2-space representation of (e1,e2) plane,');
     % 	For 2-space representation of (e1,e2) plane,
     disp('>> %  	add 1 more dimension e.');
     %  	add 1 more dimension e.
     global e p P q Q plane u d M; %/
     clf; %/
     e=e3; %/
     size = 1; %/
     va = [-size size -size size 0 2]; %/
     plane = {e+size*(e1+e2),e+size*(e1-e2),e+size*(-e1-e2),e+size*(-e1+e2)}; %/
     disp('>> e = e3; ');
     e = e3; 
     I2 = e1^e2; %/
     draw(e,'k') %/
     GAview([120 5]); %/
     GAtext(e/2-0.1*e2,'e'); %/
     GAtext(1.1*e,'0','k'); %/
     axis off; %/
     DrawPolygon(plane,'w'); %/
     axis(va); %/
     GAprompt; %/
     disp('>> % A (weighted) point may be represented by a vector.');
     % A (weighted) point may be represented by a vector.
     disp('>> p = e1/3+e2/2;');
     p = e1/3+e2/2;
     draw(p,'b'); %/
     GAtext(0.5*p+0.1*e,'p','b'); %/
     axis(va); %/
     GAprompt; %/
     disp('>> P = e+p;             ');
     P = e+p;             
     DrawHomogeneous(e,P,'b','r'); %/
     GAtext(0.5*P + 0.1*unit(p),'P','b'); %/
     axis(va); %/
     GAprompt; %/
     disp('>> % 	We might as well denote label at the point in 2-space.');
     % 	We might as well denote label at the point in 2-space.
     GAtext(1.1*P,'P','r'); %/
     GAprompt; %/
     disp('>> % 	Tilting for a better view.');
     % 	Tilting for a better view.
     GAtilt(10,5); %/
     GAprompt; %/
     clf; %/
     e=e3; %/
     size = 1; %/
     plane = {e+size*(e1+e2),e+size*(e1-e2),e+size*(-e1-e2),e+size*(-e1+e2)}; %/
     draw(e,'k'); %/
     GAtext(1.1*e,'0'); %/
     axis off; %/
     GAview([120 15]); %/
     disp('>> p = e1/4+e2/2;');
     p = e1/4+e2/2;
     disp('>> q = 2*e2/3-e1/2;');
     q = 2*e2/3-e1/2;
     disp('>> P = e+p;');
     P = e+p;
     disp('>> Q = e+q;');
     Q = e+q;
     DrawPolygon(plane,'w'); %/
     DrawHomogeneous(e,P,'b','r'); %/
     GAtext(1.1*P,'P','r'); %/
     DrawHomogeneous(e,Q,'g','r'); %/
     GAtext(1.1*Q,'Q','r'); %/
     axis(va); %/
     GAprompt; %/
     disp('>> % 	The bivector formed by P and Q can be used');
     % 	The bivector formed by P and Q can be used
     disp('>> % 	to represent the line element from P to Q.');
     % 	to represent the line element from P to Q.
     fprintf(1,'>> P^Q  ');
     input('');
     P^Q %w
     DrawBivector(P,Q,'y'); %/
     GAtext(0.25*(P+Q)-0.1*unit((P^Q)/I3),'P \wedge Q'); %/
     DrawSimplex({P,Q},'n','r'); %/
     axis(va); %/
     GAprompt; %/
     disp('>> %	The bivector can be reshaped to P^(Q-P):');
     %	The bivector can be reshaped to P^(Q-P):
     DrawBivector(P,Q-P,'g'); %/
     draw(Q-P,'r'); %/
     GAtext((Q-P)/2-0.1*e,'Q-P','r'); %/
     axis(va); %/
     title('P \wedge Q = P \wedge (Q-P)','Color','r'); %/
     disp('>> %	A line element: characterized by 2 points, ');
     %	A line element: characterized by 2 points, 
     disp('>> %	or by point and direction.');
     %	or by point and direction.
     GAprompt; %/
     GAtilt(-20,5); %/
     GAtilt(20,7.5); %/
     disp('>> %	The projective split of the bivector');
     %	The projective split of the bivector
     disp('>> %	retrieves the line parameters.');
     %	retrieves the line parameters.
     disp('>> %');
     %
     disp('>> %	The tangent vector:');
     %	The tangent vector:
     GAprompt; %/
     disp('>> PQ = P^Q');
     PQ = P^Q
     fprintf(1,'>> u = inner(e,PQ)     ');
     input('');
     u = inner(e,PQ)    %w
     DrawBivector(e,u,'m'); %/
     perp = grade((e^u)/norm(e^u)/I3,1); %/
     GAtext( (e+u)/2 + 0.1*perp,'e \wedge u'); %/
     DrawSimplex({e,e+u},'n','m'); %/
     GAtext( u/2+1.05*e+0.05*perp,'u'); %/
     axis(va); %/
     GAprompt; %/
     disp('>> %	The moment:');
     %	The moment:
     fprintf(1,'>> M = inner(e,e^PQ)     ');
     input('');
     M = inner(e,e^PQ)    %w
     d = M/u;  %/
     DrawBivector(d,u,'m'); %/
     GAtext( (d+u)/2 + 0.05*grade((d^u)/norm(d^u)/I3,1),'M'); %/
     axis(va); %/
     GAprompt; %/
     disp('>> %	The perpendicular support vector:');
     %	The perpendicular support vector:
     fprintf(1,'>> d = M/u     ');
     input('');
     d = M/u    %w
     DrawSimplex({e,e+d},'n','m'); %/
     GAtext( d/2+1.05*e+0.05*unit(u),'d'); %/
     axis(va); %/
     GAprompt; %/
     GAtilt(70,5); %/
     disp('>> ');
     
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 14 ) 
     GAps = '>>>> ';
     disp('>> % 	INTERSECTION OF SUBSPACES');
     % 	INTERSECTION OF SUBSPACES
     GAfigure; clc; %/
     disp('>> % 	INTERSECTION OF SUBSPACES');
     % 	INTERSECTION OF SUBSPACES
     disp('>> %');
     %
     global e P Q R S PQ RS M MM; %/
     e = e3;         %/
     size = 1; 	%/
     IP = {e+size*(-e1-e2),e+size*(-e1+e2),e+size*(e1+e2),e+size*(e1-e2)}; %/
     disp('>>      % LINE INTERSECTION AS MEET OF BIVECTORS');
          % LINE INTERSECTION AS MEET OF BIVECTORS
          clf; %/
     disp('>> % extra dimension of the homogeneous embedding');
     % extra dimension of the homogeneous embedding
     disp('>>      e = e3;            ');
          e = e3;            
          draw(e,'k'); %/
          GAview([30 15]); %/
     	axis off; %/
          DrawPolygon(IP,'w'); %/
     disp('>> % points P,Q,R,S');
     % points P,Q,R,S
     disp('>>      P = e- e1/3+0.9*e2;    ');
          P = e- e1/3+0.9*e2;    
          DrawHomogeneous(e,P,'n','b'); GAtext(1.07*P,'P'); %/
     disp('>>      Q = e+ 0.9*e1+e2/2;    ');
          Q = e+ 0.9*e1+e2/2;    
          DrawHomogeneous(e,Q,'n','b'); GAtext(1.07*Q,'Q'); %/
     disp('>>      R = e- e1/2-e2/4;  ');
          R = e- e1/2-e2/4;  
          DrawHomogeneous(e,R,'n','b'); GAtext(1.07*R,'R'); %/
     disp('>>      S = e+ 0.9*(e1+e2);');
          S = e+ 0.9*(e1+e2);
          DrawHomogeneous(e,S,'n','b'); GAtext(1.07*S,'S'); %/
          axis([-size size -size size (-2*size+1.2) 1.2]); %/
     disp('>> % lines');
     % lines
     disp('>>      PQ = join(P,Q);  ');
          PQ = join(P,Q);  
          %% DrawHomogeneous(e,PQ,'n','c'); %/
     GAprompt;
          DrawSimplex({P,Q},'n','c'); %/
          draw(PQ,'c'); %/
          axis([-size size -size size (-2*size+1.2) 1.2]); %/
     disp('>>      RS = join(R,S);  ');
          RS = join(R,S);  
          %% DrawHomogeneous(e,RS,'n','g'); %/
     GAprompt;
          DrawSimplex({R,S},'n','g'); %/
          draw(RS,'g'); %/
          axis([-size size -size size (-2*size+1.2) 1.2]); %/
          GAprompt; %/
     disp('>> % intersection of those lines');
     % intersection of those lines
     disp('>>      MM = meet(PQ,RS)       ');
          MM = meet(PQ,RS)       
          DrawHomogeneous(e,MM,'n','m'); %/
     	draw(MM,'m'); %/
     disp('>>      M = MM/inner(e,MM); ');
          M = MM/inner(e,MM); 
          DrawHomogeneous(e,M,'n','m'); GAtext(1.07*M,'M'); %/
          axis([-size size -size size (-2*size+1.2) 1.2]); %/
          GAprompt; %/
          GAtilt(-20,5); %/
          GAtilt(40,10); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 15 ) 
     GAps = '>>>> ';
     disp('>> %	PAPPUS'' THEOREM');
     %	PAPPUS' THEOREM
     GAfigure; clc; %/
     disp('>> %	PAPPUS'' THEOREM');
     %	PAPPUS' THEOREM
     disp('>> %');
     %
     disp('>> % 	We use Pappus merely as an example doing geometry.');
     % 	We use Pappus merely as an example doing geometry.
     disp('>> %');
     %
     disp('>> %	Pappus'' theorem (actually from projective geometry):');
     %	Pappus' theorem (actually from projective geometry):
     disp('>> %	-- Specify two pairs of collinear points');
     %	-- Specify two pairs of collinear points
     disp('>> % 	-- Cross-connect the points and intersect corresponding lines.');
     % 	-- Cross-connect the points and intersect corresponding lines.
     disp('>> % 	-> The result is three collinear points.');
     % 	-> The result is three collinear points.
     disp('>> %');
     %
     disp('>> % Two points and a factor to determine collinear point:');
     % Two points and a factor to determine collinear point:
     global P1 P2 P3 Q1 Q2 Q3 H1 H2 H3 A1 A2 A3; %/
     clf;	%/
     e = e3;	%/
     P1 = e - e1; %/
     P2 = e  + e2/4; %/
     mu = 1.5; %/
     P3 = (1-mu)*P1 + mu*P2; %/
     Q1 = e - 3*e1/4 -e2; %/
     Q2 = e + e1/6 - 1.1*e2; %/
     nu = 1.5; %/
     Q3 = (1-nu)*Q1 + nu*Q2; %/
     maxx = max([inner(e1,P1),inner(e1,P2),inner(e1,P3), inner(e1,Q1),inner(e1,Q2),inner(e1,Q3)]); %/
     minx = min([inner(e1,P1),inner(e1,P2),inner(e1,P3), inner(e1,Q1),inner(e1,Q2),inner(e1,Q3)]); %/
     maxy = max([inner(e2,P1),inner(e2,P2),inner(e2,P3), inner(e2,Q1),inner(e2,Q2),inner(e2,Q3)]); %/
     miny = min([inner(e2,P1),inner(e2,P2),inner(e2,P3), inner(e2,Q1),inner(e2,Q2),inner(e2,Q3)]); %/
     disp('>> ');
     
     GAmouse = 1; %/
     disp('>> % 	First line');
     % 	First line
     DrawHomogeneous(e,P1,'n','r'); GAtext(1.1*P1,'P_1','r'); %/
     DrawHomogeneous(e,P2,'n','r'); GAtext(1.1*P2,'P_2','r'); %/
     DrawHomogeneous(e,P3,'n','r'); GAtext(1.1*P3,'P_3','r'); %/
     DrawPolyline({P1,P2,P3},'r'); %/
     axis off; %/
     GAview([0 90]); %/
     axis([minx maxx miny maxy 0 1]);  %/
     disp('>> % 	Next: second line');
     % 	Next: second line
     GAprompt; %/
     DrawHomogeneous(e,Q1,'n','m'); GAtext(1.1*Q1,'Q_1','m'); %/
     DrawHomogeneous(e,Q2,'n','m'); GAtext(1.1*Q2,'Q_2','m'); %/
     DrawHomogeneous(e,Q3,'n','m'); GAtext(1.1*Q3,'Q_3','m'); %/
     DrawPolyline({Q1,Q2,Q3},'m'); %/
     axis([minx maxx miny maxy 0 1]);  %/
     disp('>> % 	Next: intersect the connection lines');
     % 	Next: intersect the connection lines
     GAprompt; %/
     DrawPolyline({P1,Q2},'k'); %/
     DrawPolyline({P2,Q1},'k'); %/
     disp('>> H1 = meet(join(P1,Q2),join(P2,Q1));');
     H1 = meet(join(P1,Q2),join(P2,Q1));
     disp('>> A1 = H1/inner(H1,e3);');
     A1 = H1/inner(H1,e3);
     DrawHomogeneous(e3,A1,'n','g'); GAtext(1.1*A1,'A_1','b'); %/
     axis([minx maxx miny maxy 0 1]);  %/
     GAprompt; %/
     DrawPolyline({P3,Q2},'k'); %/
     DrawPolyline({P2,Q3},'k'); %/
     H2 = meet(join(P2,Q3),join(P3,Q2)); %/
     A2 = H2/inner(H2,e3); %/
     DrawHomogeneous(e3,A2,'n','g'); GAtext(1.1*A2,'A_2','b'); %/
     axis([minx maxx miny maxy 0 1]);  %/
     GAprompt; %/
     DrawPolyline({P1,Q3},'k'); %/
     DrawPolyline({P3,Q1},'k'); %/
     H3 = meet(join(P1,Q3),join(P3,Q1)); %/
     A3 = H3/inner(H3,e3); %/
     DrawHomogeneous(e3,A3,'n','g'); GAtext(1.1*A3,'A_3','b'); %/
     axis([minx maxx miny maxy 0 1]);  %/
     GAprompt; %/
     DrawPolyline({A1,A2},'b') %/
     axis([minx maxx miny maxy 0 1]);  %/
     disp('>> % 	Next: check whether A1, A2, A3 are on a line:');
     % 	Next: check whether A1, A2, A3 are on a line:
     GAprompt; %/
     error = grade( A1^A2^A3/I3, 0); %/
     title(['A_1 \wedge A_2 \wedge A_3 = ' num2str(error) ' * I_3'],'Color','b'); %/
     GAmouse =0; %/
     disp(' ');
     disp('End of GAD sequence.  Returning to Matlab.');
end
catch ; end