                   Grupy punktowe w krystalografii: analiza macierzy 
                   © W.J. 2008-2017

Procedury i definicje systemowe:

load("nchrpl")$

load(draw)$

number: false $ ratprint: false $

[ sa60: sin(2*%pi/6), ca60: cos(2*%pi/6), sa120: sin(2*%pi/3), ca120: cos(2*%pi/3) ];

parametry1: [
line_width=1,
user_preamble="set size ratio 1",
color=black,
head_type = 'nofilled, head_length = .05, 
vector([-1,0],[2.1,0]), text:label(["y",1.1,-0.1]),
vector([0,1],[0,-2.1]), text:label(["x",-0.1,-1.1]),
implicit( x^2+y^2=1, x,-1.20,1.2, y,-1.2,1.2 ),

line_width=2,
implicit( (x-1)^2+(y-0)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x-0)^2+(y-1)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x-0)^2+(y+1)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x+1)^2+(y-0)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),

implicit( (x-sa60 )^2+(y-ca60 )^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x-sa120)^2+(y-ca120)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x-ca60 )^2+(y+sa60 )^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x-ca120)^2+(y+sa120)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),

implicit( (x+sa60 )^2+(y-ca60 )^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x+sa120)^2+(y-ca120)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x+ca60 )^2+(y-sa60 )^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),
implicit( (x+ca120)^2+(y-sa120)^2=0.0005, x,-1.2,1.2, y,-1.2,1.2 ),

line_width=1
] $

Definicje macierzy elementów symetrii:
(W poniższych definicjach znak _ oznacza minus lub kreskę ponad:
 czyli _3 trzy z kreską a _x to -x)
( L3Pz to oś sześć z kreską || Z )

E:     matrix( [ 1, 0, 0], [ 0, 1, 0], [ 0, 0, 1]) $
C:     matrix( [-1, 0, 0], [ 0,-1, 0], [ 0, 0,-1]) $
Pz:    matrix( [ 1, 0, 0], [ 0, 1, 0], [ 0, 0,-1]) $
Py:    matrix( [ 1, 0, 0], [ 0,-1, 0], [ 0, 0, 1]) $
Px:    matrix( [-1, 0, 0], [ 0, 1, 0], [ 0, 0, 1]) $
Pxy:   matrix( [ 0,-1, 0], [-1, 0, 0], [ 0, 0, 1]) $
Pyz:   matrix( [ 1, 0, 0], [ 0, 0,-1], [ 0,-1, 0]) $
Pxz:   matrix( [ 0, 0,-1], [ 0, 1, 0], [-1, 0, 0]) $
P_xy:  matrix( [ 0, 1, 0], [ 1, 0, 0], [ 0, 0, 1]) $
P_xz:  matrix( [ 0, 0, 1], [ 0, 1, 0], [ 1, 0, 0]) $
P_yz:  matrix( [ 1, 0, 0], [ 0, 0, 1], [ 0, 1, 0]) $
L2z:   Pz.C $
L2y:   Py.C $
L2x:   Px.C $
L2xy:  Pxy.C $
L2yz:  Pyz.C $
L2xz:  Pxz.C $
L2_xy: P_xy.C $
L2_xz: P_xz.C $
L2_yz: P_yz.C $
L4x:   matrix( [ 1, 0, 0], [ 0, 0,-1], [ 0, 1, 0]) $
L4y:   matrix( [ 0, 0, 1], [ 0, 1, 0], [-1, 0, 0]) $
L4z:   matrix( [ 0,-1, 0], [ 1, 0, 0], [ 0, 0, 1]) $
A4z:   L4z.C $
A4x:   L4x.C $
L3z:   matrix( [ ca120, -sa120, 0], [ sa120, ca120, 0], [ 0, 0, 1] ) $
L3Cz:  L3z.C $
L6z:   matrix( [ ca60,  -sa60,  0], [ sa60,   ca60, 0], [ 0, 0, 1] ) $
L3Pz:  L6z.C $
L3xyz: matrix( [ 0, 0, 1], [ 1, 0, 0], [ 0, 1, 0]) $

 
                                               Analiza macierzy
 

Definicja punktu:

Pkt1: matrix([0.7], [0.2], [0.4]) $

 
Tworzę macierz do analizy
(Macierz opisuje działanie nieznanego (na razie) elementu symetrii)
(Tu można wpisać macierze dowolnych elementów symetrii spośród zdefiniowanych powyżej) 
(Należy pamiętać, żeby w układzie heksagonalnym użyć tu poniżej tylko elementów symetrii związanych z Z i Y)
 

X: L4z.Py;

Sprawdzam ortogonalność
(Jeżeli wyjdzie inaczej niż 1 1 1 0 0 0 to STOP)

X[1,1]^2 + X[1,2]^2 + X[1,3]^2;

X[2,1]^2 + X[2,2]^2 + X[2,3]^2;

X[3,1]^2 + X[3,2]^2 + X[3,3]^2;

X[1,1]*X[2,1] + X[1,2]*X[2,2] + X[1,3]*X[2,3];

X[1,1]*X[3,1] + X[1,2]*X[3,2] + X[1,3]*X[3,3];

X[3,1]*X[2,1] + X[3,2]*X[2,2] + X[3,3]*X[2,3];

Obliczam wyznacznik
(Jeżeli wyjdzie inaczej niż 1 lub -1 to STOP)

Det: determinant(X);

Teraz zmieniam macierz jeżeli była to oś inwersyjna:

if (Det=-1) then (Xnew1: copymatrix(X).C) else (Xnew1: copymatrix(X));

Wyznaczam krotność osi (używam zmienionej macierzy)
(Jeżeli wyjdzie krotność niedozwolona w krystalografii to STOP)

sladXnew1: mattrace(Xnew1);

Wynik: solve(2*cos(fi)+1=sladXnew1,fi);

krotnosc: 2*%pi/rhs(Wynik[1]);

Czyli podsumowując jest to oś:
(Jeżeli wyszło -2 to chodzi oczywiście o płaszczyznę)

Os: krotnosc*Det;

 
Teraz badam kierunek:

Najpierw tworzę nową macierz

Uwaga techniczna do Maximy:
 Nieprawidłowo by było tak:  Mnew2: Mnew1;  bo by się utożsamiły w pamięci. 
 (Ze zwykłymi zmiennymi tak można a z macierzami nie.)

Xnew2: copymatrix(Xnew1);

Xnew2[1,1]: Xnew1[1,1]-1 $  Xnew2[2,2]:Xnew1[2,2]-1 $  Xnew2[3,3]:Xnew1[3,3]-1 $

Xnew2;

Liczę minory (jeśli tylko się da to wynik biorę z pierwszego rzędu):

 - minory z trzeciego wiersza:

C1t: determinant(minor(Xnew2,3,1))*(-1)^(3+1);
C2t: determinant(minor(Xnew2,3,2))*(-1)^(3+2);
C3t: determinant(minor(Xnew2,3,3))*(-1)^(3+3);

naj1: max(abs(C1t), abs(C2t), abs(C3t))$ 
if (naj1>0) then (C1:C1t, C2:C2t, C3:C3t); [C1, C2, C3];

 - minory z drugiego wiersza:

C1t: determinant(minor(Xnew2,2,1))*(-1)^(2+1);
C2t: determinant(minor(Xnew2,2,2))*(-1)^(2+2);
C3t: determinant(minor(Xnew2,2,3))*(-1)^(2+3);

naj2: max(abs(C1t), abs(C2t), abs(C3t))$ 
if (naj1=0 and naj2>0) then (C1:C1t, C2:C2t, C3:C3t)$ [C1, C2, C3];

- minory z pierwszego wiersza:

C1t: determinant(minor(Xnew2,1,1))*(-1)^(1+1);
C2t: determinant(minor(Xnew2,1,2))*(-1)^(1+2);
C3t: determinant(minor(Xnew2,1,3))*(-1)^(1+3);

naj3: max(abs(C1t), abs(C2t), abs(C3t))$ 
if (naj1=0 and naj3>naj2) then (C1:C1t, C2:C2t, C3:C3t)$ [C1, C2, C3];

Jeżeli była to oś inwersyjna to znaleziony kierunek należy zamienić na przeciwny

C1: C1*Det $ C2: C2*Det $ C3: C3*Det $

Czyli znalezione wskaźniki wynoszą:

[C1,C2,C3];

Teraz przeliczenie na kąty (w/g pierwszego wiersza).

Nie da się przeliczyć na kąty w sposób automatyczny (czyli tu w Maximie) bo
 nie da się rozróżnić [-1 1 0] od [-1/2 1/2 0] a przecież wynikają z tego różne kąty itd. 
 
Więc dzielę na dwa główne rozwiązania
 
 1. W układzie prostokątnym:
     - wystarczają tu same C1:C2:C3 (=wskaźniki) 
       (więc dalsze tu obliczenia są niepotrzebne, czyli STOP, dalej można nie czytać)

 2. W układzie heksagonalnym:
     - gdy oba elementy zorientowane są względem osi Z (np. L3z i Pz)
        to odrazu wychodzą zrozumiałe wskaźniki, czyli STOP, dalej można nie czytać
     - dla typowych (pozostałych) przypadków poniżej jest metoda na narysowanie elementu symetrii na kole projekcji

Najpierw przeliczam te wskaźniki na kąty

naj: max(abs(C1), abs(C2), abs(C3));

if (abs(naj)=3) then ([C1red:C1/(2*sqrt(3)), C2red:C2/(2*sqrt(3)), C3red:C3/(2*sqrt(3))]) else ([C1red:C1/2, C2red:C2/2, C3red:C3/2]);

if (abs(naj)=1) then ([C1red:C1, C2red:C2, C3red:C3]);

[C1red, C2red, C3red];

[ alfax:ev(realpart(acos(C1red))),  alfay:ev(realpart(acos(C2red))),  alfaz:ev(realpart(acos(C3red))) ];

Te kąty odmierzamy od osi prostokątnego układu współrzędnych:

%*360/(2*%pi);

 
Summa Summarum
(Zbieram tu wszystkie dotychczasowe wyniki)
 

[X, Xnew1, Xnew2, Os, [C1, C2, C3], [alfax*360/(2*%pi), alfay*360/(2*%pi), alfaz*360/(2*%pi)]];

Niestety Maxima nie potrafi narysować znalezionego elementu symetrii na kole projekcji
(a przynanmniej nie da się tego zrobić łatwo)...
więc na rysunku znajdzie się rozkład położeń odpowiadający działaniu znalezionego elementu symetrii
(małe niebieskie kółko oznacza położenie z ujemnym Z)
 
A element symetrii pasujący do tych położeń trzeba sobie wyobrazić

[ NPkt1:            ev(Pkt1,numer),
  NPkt2:         ev(X.NPkt1,numer),
  NPkt3:       ev(X.X.NPkt1,numer),
  NPkt4:     ev(X.X.X.NPkt1,numer), 
  NPkt5:   ev(X.X.X.X.NPkt1,numer),
  NPkt6: ev(X.X.X.X.X.NPkt1,numer)
];

wxdraw2d(

parametry1,

if NPkt1[3,1]>0 then color= red else color=blue,
if NPkt1[3,1]>0 then 
  implicit( (x-NPkt1[2,1])^2 + (y-(-NPkt1[1,1]))^2=0.005, x,-1.2,1.2, y,-1.2,1.2 ) 
 else 
  implicit( (x-NPkt1[2,1])^2 + (y-(-NPkt1[1,1]))^2=0.002, x,-1.2,1.2, y,-1.2,1.2 ),

if NPkt2[3,1]>0 then color= red else color=blue,
if NPkt2[3,1]>0 then 
  implicit( (x-NPkt2[2,1])^2 + (y-(-NPkt2[1,1]))^2=0.005, x,-1.2,1.2, y,-1.2,1.2 ) 
 else 
  implicit( (x-NPkt2[2,1])^2 + (y-(-NPkt2[1,1]))^2=0.002, x,-1.2,1.2, y,-1.2,1.2 ),

if NPkt3[3,1]>0 then color= red else color=blue,
if NPkt3[3,1]>0 then 
  implicit( (x-NPkt3[2,1])^2 + (y-(-NPkt3[1,1]))^2=0.005, x,-1.2,1.2, y,-1.2,1.2 ) 
 else 
  implicit( (x-NPkt3[2,1])^2 + (y-(-NPkt3[1,1]))^2=0.002, x,-1.2,1.2, y,-1.2,1.2 ),

if NPkt4[3,1]>0 then color= red else color=blue,
if NPkt4[3,1]>0 then 
  implicit( (x-NPkt4[2,1])^2 + (y-(-NPkt4[1,1]))^2=0.005, x,-1.2,1.2, y,-1.2,1.2 ) 
 else 
  implicit( (x-NPkt4[2,1])^2 + (y-(-NPkt4[1,1]))^2=0.002, x,-1.2,1.2, y,-1.2,1.2 ),

if NPkt5[3,1]>0 then color= red else color=blue,
if NPkt5[3,1]>0 then 
  implicit( (x-NPkt5[2,1])^2 + (y-(-NPkt5[1,1]))^2=0.005, x,-1.2,1.2, y,-1.2,1.2 ) 
 else 
  implicit( (x-NPkt5[2,1])^2 + (y-(-NPkt5[1,1]))^2=0.002, x,-1.2,1.2, y,-1.2,1.2 ),

if NPkt6[3,1]>0 then color= red else color=blue,
if NPkt6[3,1]>0 then 
  implicit( (x-NPkt6[2,1])^2 + (y-(-NPkt6[1,1]))^2=0.005, x,-1.2,1.2, y,-1.2,1.2 ) 
 else 
  implicit( (x-NPkt6[2,1])^2 + (y-(-NPkt6[1,1]))^2=0.002, x,-1.2,1.2, y,-1.2,1.2 )

);