find_inv_G

This function computes the invariants of a one turn map matrix.
The resulting invariant matrices G1,G2,G3 satisfy
 M^T G_a M = G_a for a=1,2,3
 Algorithm from PhD thesis of B. Nash

0001 function [G1 G2 G3]=find_inv_G(m66)
0002 %This function computes the invariants of a one turn map matrix.
0003 %The resulting invariant matrices G1,G2,G3 satisfy
0004 % M^T G_a M = G_a for a=1,2,3
0005 % Algorithm from PhD thesis of B. Nash
0006 
0007 [V,d]=eigs(m66);
0008 v1=V(:,1);
0009 v2=V(:,2);
0010 v3=V(:,3);
0011 v4=V(:,4);
0012 v5=V(:,5);
0013 v6=V(:,6);
0014 
0015 %compute norms of each:
0016 for j=1:6
0017 n(j)=-(1/2)*i*V(:,j)'*jmat(3)*V(:,j);
0018 end
0019 
0020 %put positive modes before negative modes (swap columns if first n in pair
0021 %is negative)
0022 for j=1:3
0023     if n(2*j-1) <0 
0024         V(:,[2*j-1,2*j]) = V(:,[2*j,2*j-1]);
0025         n([2*j-1,2*j]) = n([2*j,2*j-1]);
0026     end
0027 end
0028 %now normalize each vector
0029 
0030 for j=1:6
0031     Vn(:,j)=V(:,j)/sqrt(abs(n(j)));
0032 end
0033 
0034 %find the vecs that project most onto x,y,z, and reorder
0035 %call these vectors Vxyz, ordered via (vx,vmx,vy,vmy,vz,vmz)
0036 Vxyz(:,1)=(1/sqrt(2))*[i 1 0 0 0 0]';
0037 Vxyz(:,2)=i*conj(Vxyz(:,1));
0038 Vxyz(:,3)=(1/sqrt(2))*[0 0 i 1 0 0]';
0039 Vxyz(:,4)=i*conj(Vxyz(:,3));
0040 Vxyz(:,5)=(1/sqrt(2))*[0 0 0 0 i 1]';
0041 Vxyz(:,6)=i*conj(Vxyz(:,5));
0042 
0043 for j=1:6
0044     for k=1:6
0045         nn(j,k) = (1/2)*abs(sqrt(-i*Vn(:,j)'*jmat(3)*Vxyz(:,k)));
0046     end
0047 end
0048 
0049 %now use nn to order pairs
0050 n1 = nn(1,[1,3,5]);
0051 n3 = nn(3,[1,3,5]);
0052 n5 = nn(5,[1,3,5]);
0053 
0054 [vals1,sind1]=sort(n1,'descend');
0055 [vals3,sind3]=sort(n3,'descend');
0056 [vals5,sind5]=sort(n5,'descend');
0057 
0058 orderp = [2*sind1(1)-1,2*sind3(1)-1,2*sind5(1)-1];
0059 ordern = [2*sind1(1),2*sind3(1),2*sind5(1)];
0060 
0061 V_ordered = Vn;
0062 V_ordered(:,[1,3,5])=V_ordered(:,orderp); %reorder positive modes
0063 V_ordered(:,[2,4,6])=V_ordered(:,ordern); %reorder negative modes
0064 
0065 
0066 v1=V_ordered(:,1);
0067 v2=V_ordered(:,3);
0068 v3=V_ordered(:,5);
0069 
0070 
0071 G1 = -(1/2)*jmat(3)*(conj(v1)*v1.'+v1*v1')*jmat(3);
0072 G2 = -(1/2)*jmat(3)*(conj(v2)*v2.'+v2*v2')*jmat(3);
0073 G3 = -(1/2)*jmat(3)*(conj(v3)*v3.'+v3*v3')*jmat(3);