gablock.m

这是几何代数的matlab工具包

function GAblock(GAn)
%GAblock: run sample code in tutorial.

try
if ( GAn == 1 ) 
  GAps = 'GAblock >> ';
     disp('>>      % ORTHOGONALIZATION');
     % ORTHOGONALIZATION
     disp('>>      clf;');
     clf;
     disp('>>      u = e1+e2;');
     u = e1+e2;
     disp('>>      v = 0.3*e1 + 0.6*e2 - 0.8*e3;');
     v = 0.3*e1 + 0.6*e2 - 0.8*e3;
     disp('>>      w = e1 -0.2*e2 + 0.5*e3;');
     w = e1 -0.2*e2 + 0.5*e3;
     disp('>>      up = u;');
     up = u;
     disp('>>      vp = (v^up)/up;');
     vp = (v^up)/up;
     disp('>>      wp = (w^up^vp)/(up^vp);');
     wp = (w^up^vp)/(up^vp);
     disp('>>      draw(u); draw(v); draw(w);                 %% The original vectors ...');
     draw(u); draw(v); draw(w);                 %% The original vectors ...
  GAprompt;
     disp('>>      draw(up,''r''); draw(vp,''r''); draw(wp,''r'');  %% ... and orthognalized');
     draw(up,'r'); draw(vp,'r'); draw(wp,'r');  %% ... and orthognalized
  GAprompt;
     disp('>>      GAorbiter');
     GAorbiter
  disp(' ');    disp('End of GAblock sequence.  Returning to Matlab.');
elseif ( GAn == 2 ) 
  GAps = 'GAblock >> ';
     disp('>>      % ROTATION EXERCISE');
     % ROTATION EXERCISE
     disp('>>      clf;');
     clf;
     disp('>>      R1 = gexp(-I3*e1*pi/2/2);');
     R1 = gexp(-I3*e1*pi/2/2);
     disp('>>      R2 = gexp(-I3*e2*pi/2/2);');
     R2 = gexp(-I3*e2*pi/2/2);
     disp('>>      R =  R2*R1;');
     R =  R2*R1;
     disp('>>      a = e1+e2;');
     a = e1+e2;
     disp('>>      Ra = R*a/R;');
     Ra = R*a/R;
     disp('>>      RRa = R*Ra/R;');
     RRa = R*Ra/R;
     disp('>>      draw(a); draw(Ra,''m''); draw(RRa,''r''); %% Draw the objects');
     draw(a); draw(Ra,'m'); draw(RRa,'r'); %% Draw the objects
  GAprompt;
     disp('>>      axisR = unit( -GAZ(sLog(R))/I3 );');
     axisR = unit( -GAZ(sLog(R))/I3 );
     disp('>>      draw(axisR,''g''); % Draw the axis of rotation');
     draw(axisR,'g'); % Draw the axis of rotation
  GAps = 'Refer to the tutorial before continuing >> ';
  disp(' ');  GAprompt;
  GAps = 'GAblock >> ';
     disp('>>      draw(dual(axisR),''g'');');
     draw(dual(axisR),'g');
  disp(' ');    disp('End of GAblock sequence.  Returning to Matlab.');
elseif ( GAn == 3 ) 
  GAps = 'GAblock >> ';
     disp('>>      % INTERPOLATION OF ORIENTATIONS');
     % INTERPOLATION OF ORIENTATIONS
     disp('>>      clf;');
     clf;
     disp('>>      RA = gexp(-I3*e1*pi/2/2);');
     RA = gexp(-I3*e1*pi/2/2);
     disp('>>      RB = gexp(-I3*e2*pi/2/2);');
     RB = gexp(-I3*e2*pi/2/2);
     disp('>>      Rtot =  RB/RA');
     Rtot =  RB/RA
     disp('>>      n = 8;                          % we rotate in 8 steps');
     n = 8;                          % we rotate in 8 steps
     disp('>>      R = gexp(sLog(Rtot)/n);');
     R = gexp(sLog(Rtot)/n);
  GAps = 'Refer to the tutorial before continuing >> ';
  disp(' ');  GAprompt;
  GAps = 'GAblock >> ';
     disp('>>      u = e1+e2-e3;');
     u = e1+e2-e3;
     disp('>>      v = e1+e3;');
     v = e1+e3;
     disp('>>      view = [-0.6  2.5  -1    1.16  -2  1.1];  % select the view');
     view = [-0.6  2.5  -1    1.16  -2  1.1];  % select the view
     disp('>>            % === initial orientation:');
           % === initial orientation:
     disp('>>      DrawBivector(RA*u/RA,RA*v/RA,''b'');  axis(view); GAview([30 30]); %%');
     DrawBivector(RA*u/RA,RA*v/RA,'b');  axis(view); GAview([30 30]); %%
  GAprompt;
     disp('>>            % === final orientation:');
           % === final orientation:
     disp('>>      DrawBivector(RB*u/RB,RB*v/RB,''g'');  axis(view);                  %%');
     DrawBivector(RB*u/RB,RB*v/RB,'g');  axis(view);                  %%
  GAprompt;
     disp('>>      axisR = unit(GAZ(-sLog(R)/I3));   % reorientation axis: ');
     axisR = unit(GAZ(-sLog(R)/I3));   % reorientation axis: 
     disp('>>      draw(axisR,''r'');                 %% displayed for visualization ');
     draw(axisR,'r');                 %% displayed for visualization 
  GAprompt;
     disp('>>            % === display of the 7 intermediate orientations ');
           % === display of the 7 intermediate orientations 
     disp('>>      Ri = RA;');
     Ri = RA;
     for i=1:n-1
     disp('>>      for i=1:n-1');
     disp(['i = ', num2str(i)])
     disp('>>          Ri = R*Ri;');
         Ri = R*Ri;
     disp('>>          ui = Ri*u/Ri;');
         ui = Ri*u/Ri;
     disp('>>          vi = Ri*v/Ri;');
         vi = Ri*v/Ri;
     disp('>>          DrawBivector(ui,vi);');
         DrawBivector(ui,vi);
     disp('>>          drawnow;');
         drawnow;
     disp('>>      end');
     end
     disp('>>      GAorbiter(125);');
     GAorbiter(125);
  disp(' ');    disp('End of GAblock sequence.  Returning to Matlab.');
elseif ( GAn == 4 ) 
  GAps = 'GAblock >> ';
     disp('>>      % LINE INTERSECTS LINE');
     % LINE INTERSECTS LINE
     disp('>>      p = e2; u = 0.2*e2 + e1;');
     p = e2; u = 0.2*e2 + e1;
     disp('>>      q = e1; v = e2-2*e1;');
     q = e1; v = e2-2*e1;
     disp('>>      clf; ');
     clf; 
     disp('>>      draw(p); GAview([0 90]);');
     draw(p); GAview([0 90]);
     disp('>>      DrawPolyline({p-2*u,p+2*u});');
     DrawPolyline({p-2*u,p+2*u});
     disp('>>      draw(q,''g''); DrawPolyline({q-v,q+2*v},''g'');');
     draw(q,'g'); DrawPolyline({q-v,q+2*v},'g');
  GAps = 'Refer to the tutorial before continuing >> ';
  disp(' ');  GAprompt;
  GAps = 'GAblock >> ';
     disp('>>      U = (q^v/(u^v)) * u');
     U = (q^v/(u^v)) * u
     disp('>>      V = (p^u/(v^u)) * v');
     V = (p^u/(v^u)) * v
     disp('>>      draw(U,''m'')          %% Draw U');
     draw(U,'m')          %% Draw U
  GAprompt;
     disp('>>      draw(V, ''m'')         %% Draw V');
     draw(V, 'm')         %% Draw V
  GAprompt;
     disp('>>      draw(U+V, ''r'' )      % Draw U+V');
     draw(U+V, 'r' )      % Draw U+V
  disp(' ');    disp('End of GAblock sequence.  Returning to Matlab.');
elseif ( GAn == 5 ) 
  GAps = 'GAblock >> ';
     disp('>>      % PROJECTION, REJECTION');
     % PROJECTION, REJECTION
     disp('>>      x = e1 + e2/2+e3;');
     x = e1 + e2/2+e3;
     disp('>>      A = e2 + e3/3;   % a linear subspace');
     A = e2 + e3/3;   % a linear subspace
     disp('>>      xA = geoall(x,A);');
     xA = geoall(x,A);
     disp('>>      clf; draw(x,''b''); draw(A,''g'');       %% Draw A and x');
     clf; draw(x,'b'); draw(A,'g');       %% Draw A and x
  GAprompt;
     disp('>>      DrawPolyline({xA.rej,xA.rej+xA.proj,xA.proj},''k'');');
     DrawPolyline({xA.rej,xA.rej+xA.proj,xA.proj},'k');
     disp('>>      draw(xA.rej,''m''); draw(xA.proj,''m''); %% Draw rej and proj');
     draw(xA.rej,'m'); draw(xA.proj,'m'); %% Draw rej and proj
  GAprompt;
     disp('>>      distxA = norm(xA.rej)');
     distxA = norm(xA.rej)
     disp('>>      tanglexA = xA.rej/xA.proj');
     tanglexA = xA.rej/xA.proj
     disp('>>      anglexA = atan(norm(tanglexA))*180/pi %%');
     anglexA = atan(norm(tanglexA))*180/pi %%
  GAprompt;
     disp('>>      B = e1^A;             % A planar subspace (containing A)');
     B = e1^A;             % A planar subspace (containing A)
     disp('>>      xB = geoall(x,B);		    ');
     xB = geoall(x,B);		    
     disp('>>      draw(B,''y'');          %%');
     draw(B,'y');          %%
  GAprompt;
     disp('>>      DrawPolyline({xB.rej,xB.rej+xB.proj,xB.proj},''k'');');
     DrawPolyline({xB.rej,xB.rej+xB.proj,xB.proj},'k');
     disp('>>      draw(xB.rej,''r''); draw(xB.proj,''r''); %% rej and proj for B');
     draw(xB.rej,'r'); draw(xB.proj,'r'); %% rej and proj for B
  GAprompt;
     disp('>>      distxB = norm(xB.rej)');
     distxB = norm(xB.rej)
     disp('>>      tanglexB = xB.rej/xB.proj');
     tanglexB = xB.rej/xB.proj
     disp('>>      anglexB = atan(norm(tanglexB))*180/pi');
     anglexB = atan(norm(tanglexB))*180/pi
     disp('>>      side = sign((x^B)/I3)');
     side = sign((x^B)/I3)
  disp(' ');    disp('End of GAblock sequence.  Returning to Matlab.');
elseif ( GAn == 6 ) 
  GAps = 'GAblock >> ';
     disp('>>      % MEET, JOIN ');
     % MEET, JOIN 
     disp('>>      A = e2^(e1+e3);');
     A = e2^(e1+e3);
     disp('>>      B = e1^(e2+e3/2);');
     B = e1^(e2+e3/2);
     disp('>>      cAB = inner(A,B);  pAB = cAB/B; rAB = A - pAB;');
     cAB = inner(A,B);  pAB = cAB/B; rAB = A - pAB;
     disp('>>      clf; draw(A,''b''); draw(B,''g''); %%');
     clf; draw(A,'b'); draw(B,'g'); %%
  GAprompt;
     disp('>>      draw(cAB,''r''); draw(pAB,''c''); draw(rAB,''m'');');
     draw(cAB,'r'); draw(pAB,'c'); draw(rAB,'m');
  GAps = 'Refer to the tutorial before continuing >> ';
  disp(' ');  GAprompt;
  GAps = 'GAblock >> ';
     disp('>>      mAB = meet(A,B)');
     mAB = meet(A,B)
     disp('>>      clf; draw(A,''b''); draw(B,''g''); draw(mAB,''y'');');
     clf; draw(A,'b'); draw(B,'g'); draw(mAB,'y');
  GAps = 'Refer to the tutorial before continuing >> ';
  disp(' ');  GAprompt;
  GAps = 'GAblock >> ';
     disp('>>       jAB = join(A,B)');
      jAB = join(A,B)
     disp('>>       draw(jAB,''k'');');
      draw(jAB,'k');
  disp(' ');    disp('End of GAblock sequence.  Returning to Matlab.');
elseif ( GAn == 7 ) 
  GAps = 'GAblock >> ';
     disp('>>      % MEET AND JOIN DECOMPOSED');