gad.m

这是几何代数的matlab工具包

     input('');
     b/(b*b)    %w
     GAprompt; %/
     disp('>> % 	Inverse of unit vector:');
     % 	Inverse of unit vector:
     fprintf(1,'>> 1/e1     ');
     input('');
     1/e1    %w
     GAprompt; %/
     disp('>> %	The geometric product is invertible; ');
     %	The geometric product is invertible; 
     disp('>> %	From (a*b) and b, retrieve a = e1+e3 :');
     %	From (a*b) and b, retrieve a = e1+e3 :
     disp('>> a*b');
     a*b
     disp('>> b');
     b
     fprintf(1,'>> (a*b)/b     ');
     input('');
     (a*b)/b    %w
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 8 ) 
     GAps = '>>>> ';
     disp('>> % 	PROJECTION, REJECTION');
     % 	PROJECTION, REJECTION
     GAfigure; clc; %/
     disp('>> % 	PROJECTION, REJECTION');
     % 	PROJECTION, REJECTION
     disp('>> %');
     %
     global x a A p r; %/
     clf; %/
     disp('>> % 	Take two vectors x and a:');
     % 	Take two vectors x and a:
     fprintf(1,'\n');     disp('>> x = e1 + e2 + e3;');
     x = e1 + e2 + e3;
     disp('>> a = e2/4 + e3/2;');
     a = e2/4 + e3/2;
     draw(x,'r'); %/
     random=unit(pi*e1+pi/exp(1)*e2+exp(1)*e3); %/
     GAtext(0.7*x+0.1*unit(grade((random^x)/x,1)),'x'); %/
     draw(a,'g'); %/
     GAtext(0.7*a-0.1*unit(grade((random^a)/a,1)),'a'); %/
     va = [-0.05 1.1 0 1.25 -0.25 1.5]; %/
     axis(va); %/
     axis('vis3d'); %/
     axis off; %/
     GAprompt; %/
     GAorbiter(60,3); %/
     axis(va); %/
     GAprompt; %/
     disp('>> % 	projection of x onto a: ');
     % 	projection of x onto a: 
     fprintf(1,'\n');     disp('>> p = inner(x,a)/a ');
     p = inner(x,a)/a 
     draw(p,'m'); %/
     draw(a,'g'); %/
     DrawPolyline({x,p},'k'); %/
     GAtext(0.7*p-0.1*unit(grade((random^p)/p,1)),'(x \bullet a) / a'); %/ 
     axis(va); %/
     GAprompt; %/
     disp('>> % 	rejection of x by a:');
     % 	rejection of x by a:
     fprintf(1,'\n');     disp('>> r = (x^a)/a  ');
     r = (x^a)/a  
     draw(r,'m'); %/
     DrawPolyline({x,r},'k'); %/
     GAtext(0.5*r-0.1*unit(grade((random^r)/r,1)),'r = x \wedge a) / a'); %/ 
     axis(va); %/
     GAprompt; %/
     axis(va); %/
     GAorbiter(360,10); %/ 
     axis(va); %/
     GAprompt; %/
     disp('>> % 	Explanation: ');
     % 	Explanation: 
     axis(va); %/
     GAprompt; %/
     disp('>> %	span x^a,');
     %	span x^a,
     DrawBivector(x,a,'y'); %/
     GAtext((x+a)/2 + 0.1*(x^a)/I3,'x \wedge a'); %/
     axis(va); %/
     disp('>> % 	reshape orthogonally,');
     % 	reshape orthogonally,
     GAprompt; %/
     DrawBivector(r,a,'y'); %/
     GAtext((r+a)/2 + 0.1*(x^a)/I3,'x \wedge a'); %/
     axis(va); %/
     GAprompt; %/
     title('x \wedge a = r \wedge a = r \bullet a + r \wedge a = r a','FontSize',12);%/
     disp('>> % 	division by a then gives the rejection.');
     % 	division by a then gives the rejection.
     axis(va); %/
     GAprompt; %/
     fprintf(1,'\n');     disp('>> %	Projection and rejection formulas also work for bivectors:	');
     %	Projection and rejection formulas also work for bivectors:	
     GAprompt; %/
     clf; %/
     disp('>> x = e1/3 - e2/2 + e3;');
     x = e1/3 - e2/2 + e3;
     disp('>> A = e1^(e2/3 + e3/2);');
     A = e1^(e2/3 + e3/2);
     draw(x,'r'); %/
     axis off; %/
     axis([-0.5 0.5 -0.8 0.2 -0.4 1]); %/
     random=unit(pi*e1+pi/exp(1)*e2+exp(1)*e3); %/
     GAtext(0.7*x-0.1*unit(grade((random^x)/x,1)),'x'); %/
     draw(A,'g'); %/
     GAtext(-0.4*unit(grade(inner(random,A)/A,1))+0.1*unit(A/I3),'A'); %/ 
     va = [-0.5 0.5 -0.81 0.31 -0.4 1]; %/
     axis(va); %/
     axis('vis3d'); %/
     GAprompt; %/
     GAorbiter(360,4); %/
     axis off; %/
     disp('>> % 	projection of x onto A: ');
     % 	projection of x onto A: 
     fprintf(1,'\n');     disp('>> p = inner(x,A)/A  ');
     p = inner(x,A)/A  
     draw(p,'m'); %/
     DrawPolyline({x,p},'k'); %/
     GAtext(0.9*p+0.1*unit(grade((random^p)/p,1)),'(x \bullet A) / A'); %/ 
     axis(va); %/
     GAprompt; %/
     disp('>> % 	rejection of x by a:');
     % 	rejection of x by a:
     fprintf(1,'\n');     disp('>> r = (x^A)/A  ');
     r = (x^A)/A  
     draw(r,'m'); %/
     DrawPolyline({x,r},'k'); %/
     GAtext(0.7*r-0.1*unit(grade((random^r)/r,1)),'(x \wedge A) / A'); %/ 
     axis(va); %/
     GAprompt; %/
     GAorbiter(340,10); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 9 ) 
     GAps = '>>>> ';
     disp('>> %	ROTATIONS IN A PLANE');
     %	ROTATIONS IN A PLANE
     GAfigure; clc; %/
     disp('>> %	ROTATIONS IN A PLANE');
     %	ROTATIONS IN A PLANE
     disp('>> %');
     %
     global a b c R x; %/
     disp('>> clf;');
     clf;
     va = [-1 1 -0.2 1.3 -0.5 0.5]; %/
     disp('>> a = unit(e1+e2/3);');
     a = unit(e1+e2/3);
     disp('>> b = unit(e1+e2/2);');
     b = unit(e1+e2/2);
     draw(a,'b'); %/
     GAtext(1.1*a,'a'); %/
     draw(b,'b'); %/
     bt = GAtext(1.1*b,'b'); %/
     GAview([0 90]); %/
     axis off; %/
     axis(va); %/
     GAprompt; %/
     disp('>> % 	Two vectors determine a rotation/dilation');
     % 	Two vectors determine a rotation/dilation
     fprintf(1,'>> R = b/a     ');
     input('');
     R = b/a    %w
     disp('>> %	This is an operator. applicable to vectors in the plane:');
     %	This is an operator. applicable to vectors in the plane:
     disp('>> c = 1.1*e2-e1/2;');
     c = 1.1*e2-e1/2;
     draw(c,'m');  %/
     GAtext(1.05*c,'c');  %/
     axis(va); %/
     fprintf(1,'>> x = R*c	    ');
     input('');
     x = R*c	   %w
     draw(x,'r'); %/
     GAtext(x-0.2*e1,'x = (b/a)c');%/
     axis(va); %/
     disp('>> %	x is to c as b is to a');
     %	x is to c as b is to a
     GAprompt;  %/
     disp('>> %	Repeated application to a:');
     %	Repeated application to a:
     title('R = b/a ','Color','r'); %/
     GAprompt; %/
     draw(R*a,'r');  		%/
     delete(bt); 		%/
     GAtext(1.1*R*a,'R a'); 		%/
     axis(va); %/
     GAprompt; %/
     draw(R*R*a,'r'); 		%/
     GAtext(1.1*R*R*a,'R^2 a'); 		%/
     axis(va); %/
     GAprompt; %/
     draw(R*R*R*a,'r'); 		%/
     GAtext(1.1*R*R*R*a,'R^3 a'); 		%/
     axis(va); %/
     GAprompt; %/
     draw(R*R*R*R*a,'r'); 		%/
     GAtext(1.1*R*R*R*R*a,'R^4 a'); 		%/
     axis(va); %/
     GAprompt; %/
     draw(R*R*R*R*R*a,'r'); 		%/
     GAtext(1.1*R*R*R*R*R*a,'R^5 a'); 		%/
     axis(va); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 10 ) 
     GAps = '>>>> ';
     disp('>> %	ROTATION BY ROTORS');
     %	ROTATION BY ROTORS
     GAfigure; clc; %/
     disp('>> %	ROTATION BY ROTORS');
     %	ROTATION BY ROTORS
     global plane angle i phi Raxis; %/
     clf; %/
     disp('>> %');
     %
     disp('>> % 	Making a rotor:');
     % 	Making a rotor:
     disp('>> %');
     %
     disp('>> plane = e1^e2;');
     plane = e1^e2;
     disp('>> angle = pi/4;');
     angle = pi/4;
     i = plane; %/
     phi = angle; %/
     fprintf(1,'>> R = gexp(-plane*angle/2)   ');
     input('');
     R = gexp(-plane*angle/2) %w %%
     GAprompt;
     disp('>> % 	Applying the rotor:');
     % 	Applying the rotor:
     disp('>> x = e1 + e3;');
     x = e1 + e3;
     draw(x,'r'); %/
     axis([-1 1 -1 1 -1 1]); %/
     axis off; %/
     GAtext(0.7*x+0.1*unit(grade(inner(x,plane)/plane,1)),'x'); %/
     draw(plane*angle); %/
     GAview([15 30]); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAtext(0.1*plane/I3-0.3*unit(grade(inner(x,plane)/plane,1)),'i \phi'); %/
     fprintf(1,'>> r = R*x/R  ');
     input('');
     r = R*x/R %w
     draw(r,'m'); %/
     axis([-1 1 -1 1 -1 1]); %/
     title(['rotor R = e^{-i \phi /2}'],'Color','b'); %/
     GAtext(0.9*r+0.1*unit(grade(inner(r,plane)/plane,1)),'R x R^{-1}'); %/
     GAprompt; %/
     title(''); %/
     GAorbiter(-360,5); %/
     disp('>> % In 3-space, you could characterize the rotation by an axis:');
     % In 3-space, you could characterize the rotation by an axis:
     fprintf(1,'>> Raxis = plane/I3  ');
     input('');
     Raxis = plane/I3 %w
     draw(Raxis,'k'); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAtext(1.1*Raxis,'i \phi / I_3'); %/
     GAprompt; %/
     disp('>> % Now rotate a bivector.');
     % Now rotate a bivector.
     disp('>> B = 0.7*x^(e1+e2);');
     B = 0.7*x^(e1+e2);
     draw(B,'w'); %/
     axis([-1 1 -1 1 -1 1]); %/
     Blabel = -0.75*unit(meet(B,i)); %/
     GAtext(Blabel,'B'); %/
     fprintf(1,'>> RB = R*B/R  ');
     input('');
     RB = R*B/R %w
     draw(RB,'y'); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAtext(R*Blabel/R,'R B R^{-1}'); %/
     GAprompt; %/
     GAorbiter(-400,10); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 11 ) 
     GAps = '>>>> ';
     disp('>> %	QUATERNIONS IN GEOMETRIC ALGEBRA');
     %	QUATERNIONS IN GEOMETRIC ALGEBRA
     GAfigure; clc; %/
     disp('>> %	QUATERNIONS IN GEOMETRIC ALGEBRA');
     %	QUATERNIONS IN GEOMETRIC ALGEBRA
     global i j k u q bivector Rangle Raxis; %/
     clf; %/
     disp('>> % The basic ''vectors'' in quaternions are unit bivectors.');
     % The basic 'vectors' in quaternions are unit bivectors.
     fprintf(1,'>> i = e1*I3     ');
     input('');
     i = e1*I3    %w
     fprintf(1,'>> j = -e2*I3     ');
     input('');
     j = -e2*I3    %w
     fprintf(1,'>> k = e3*I3     ');
     input('');
     k = e3*I3    %w
     disp('>> % The quaternion product is the geometric prodcut:');
     % The quaternion product is the geometric prodcut:
     fprintf(1,'>> i*i 		 ');
     input('');
     i*i 		%w
     fprintf(1,'>> j*j 		 ');
     input('');
     j*j 		%w
     fprintf(1,'>> k*k 		 ');
     input('');
     k*k 		%w
     fprintf(1,'>> i*j 		 ');
     input('');
     i*j 		%w
     fprintf(1,'>> i*j*k 		 ');
     input('');
     i*j*k 		%w
     disp('>> % A (unit) quaternion is a rotor:');
     % A (unit) quaternion is a rotor:
     fprintf(1,'>> q = 1 + i +j +k  ');
     input('');
     q = 1 + i +j +k %w
     GAprompt; %/
     fprintf(1,'>> u = q/norm(q) 	 ');
     input('');
     u = q/norm(q) 	%w
     bivector = sLog(u); %/
     Rangle = norm(bivector); %/
     Raxis = bivector/I3; %/
     disp('>> % A quaternion can be applied to a vector, bivector etc.,');
     % A quaternion can be applied to a vector, bivector etc.,
     disp('>> % without converting it to a matrix first');
     % without converting it to a matrix first
     disp('>> % (and without normalization).');
     % (and without normalization).
     clf; %/
     x = e1; %/
     draw(x,'b'); %/
     axis off; %/
     GAtext(1.1*x,'x','b'); %/
     draw(bivector,'r'); %/
     draw(Raxis,'k'); %/
     label = -0.5*unit(grade(inner(x + q*x/q,bivector)/bivector,1))+ 0.1*unit(Raxis); %/
     GAtext(label, 'log(q)'); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAview([45 30]); %/
     GAprompt; %/
     fprintf(1,'>> Rx = q*x/q   ');
     input('');
     Rx = q*x/q  %w
     draw(Rx,'g'); %/
     GAtext(1.1*Rx,'q x q^{-1}','k'); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAprompt; %/
     GAorbiter(360,10); %/
     GAprompt; %/
     disp('>> % And it can be applied directly to bivectors:');
     % And it can be applied directly to bivectors:
     disp('>> B = x^(e2+e3); ');
     B = x^(e2+e3); 
     draw(B,'b'); %/
     axis([-1 1 -1 1 -1 1]); %/
     label = grade(meet(B,bivector),1); %/
     GAtext(0.75*label,'B','b'); %/
     fprintf(1,'>> RB = q*B/q   ');
     input('');
     RB = q*B/q  %w
     draw(RB,'g'); %/
     GAtext(0.75*q*label/q,'q B q^{-1}','k'); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAprompt; %/
     GAorbiter(360,10); %/
     disp('>> ');
     
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 12 ) 
     GAps = '>>>> ';
     disp('>> % 	INTERPOLATION OF ORIENTATIONS');
     % 	INTERPOLATION OF ORIENTATIONS
     GAfigure; clc; %/