/*
Space Group Info's (c) 1994-96 Ralf W. Grosse-Kunstleve
*/
/*
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
The GNU GPL can also be found at http://www.gnu.org
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define SGCLIB_C__
#include "sginfo.h"
static const char *Err_Ill_SMx_in_List =
"Error: Illegal SeitzMx in list";
void SetSgError(const char *msg)
{
if (SgError == NULL) SgError = msg;
}
int iModPositive(int ix, int iy)
{
if (iy > 0)
{
ix %= iy;
if (ix < 0) ix += iy;
}
return ix;
}
static void SwapRTMx(T_RTMx *Mx_a, T_RTMx *Mx_b)
{
int i;
T_RTMx BufMx;
for (i = 0; i < 12; i++) BufMx.a[i] = Mx_a->a[i];
for (i = 0; i < 12; i++) Mx_a->a[i] = Mx_b->a[i];
for (i = 0; i < 12; i++) Mx_b->a[i] = BufMx.a[i];
}
static void CopyRotMxInfo(T_RotMxInfo *target, const T_RotMxInfo *source)
{
memcpy(target, source, sizeof (*target));
}
static void SwapRotMxInfo(T_RotMxInfo *RMx_a, T_RotMxInfo *RMx_b)
{
T_RotMxInfo Buffer;
memcpy(&Buffer, RMx_a, sizeof (Buffer));
memcpy(RMx_a, RMx_b, sizeof (Buffer));
memcpy(RMx_b, &Buffer, sizeof (Buffer));
}
int traceRotMx(const int *RotMx)
{
return RotMx[0] + RotMx[4] + RotMx[8];
}
int deterRotMx(const int *RotMx)
{
int det;
det = RotMx[0] * (RotMx[4] * RotMx[8] - RotMx[5] * RotMx[7]);
det -= RotMx[1] * (RotMx[3] * RotMx[8] - RotMx[5] * RotMx[6]);
det += RotMx[2] * (RotMx[3] * RotMx[7] - RotMx[4] * RotMx[6]);
return det;
}
void RotMx_t_Vector(int *R_t_V, const int *RotMx, const int *Vector, int FacTr)
{
const int *vec;
if (FacTr > 0)
{
vec = Vector;
*R_t_V = *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec;
*R_t_V %= FacTr; if (*R_t_V < 0) *R_t_V += FacTr;
R_t_V++;
vec = Vector;
*R_t_V = *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec;
*R_t_V %= FacTr; if (*R_t_V < 0) *R_t_V += FacTr;
R_t_V++;
vec = Vector;
*R_t_V = *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec++;
*R_t_V += *RotMx * *vec;
*R_t_V %= FacTr; if (*R_t_V < 0) *R_t_V += FacTr;
}
else
{
vec = Vector;
*R_t_V = *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec++;
*R_t_V++ += *RotMx++ * *vec;
vec = Vector;
*R_t_V = *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec++;
*R_t_V++ += *RotMx++ * *vec;
vec = Vector;
*R_t_V = *RotMx++ * *vec++;
*R_t_V += *RotMx++ * *vec++;
*R_t_V += *RotMx * *vec;
}
}
void RotMxMultiply(int *rmxab, const int *rmxa, const int *rmxb)
{
const int *a, *b;
/* no loops to be as fast as possible */
/* -funroll-loops in gcc would have had the same effect ;-) (sxm)*/
a = rmxa;
b = rmxb;
*rmxab = *a++ * *b; b += 3; /* r11 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a * *b; b -= 5;
rmxab++;
a = rmxa;
*rmxab = *a++ * *b; b += 3; /* r12 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a * *b; b -= 5;
rmxab++;
a = rmxa;
*rmxab = *a++ * *b; b += 3; /* r13 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a++ * *b; b = rmxb;
rmxab++;
rmxa = a;
*rmxab = *a++ * *b; b += 3; /* r21 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a * *b; b -= 5;
rmxab++;
a = rmxa;
*rmxab = *a++ * *b; b += 3; /* r22 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a * *b; b -= 5;
rmxab++;
a = rmxa;
*rmxab = *a++ * *b; b += 3; /* r23 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a++ * *b; b = rmxb;
rmxab++;
rmxa = a;
*rmxab = *a++ * *b; b += 3; /* r31 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a * *b; b -= 5;
rmxab++;
a = rmxa;
*rmxab = *a++ * *b; b += 3; /* r32 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a * *b; b -= 5;
rmxab++;
a = rmxa;
*rmxab = *a++ * *b; b += 3; /* r33 */
*rmxab += *a++ * *b; b += 3;
*rmxab += *a * *b;
}
void RotateRotMx(int *RotMx, const int *RMx, const int *InvRMx)
{
int BufMx[9];
RotMxMultiply(BufMx, RotMx, InvRMx);
RotMxMultiply(RotMx, RMx, BufMx);
}
void SeitzMxMultiply(T_RTMx *smxab, const T_RTMx *smxa, const T_RTMx *smxb)
{
const int *ar, *a, *b, *bt;
int *ab;
/* no loops to be as fast as possible */
ar = smxa->a;
a = smxa->a;
b = smxb->a;
ab = smxab->a;
*ab = *a++ * *b; b += 3; /* r11 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r12 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r13 */
*ab += *a++ * *b; b += 3;
*ab += *a++ * *b; b = smxb->a;
ab++;
ar = a;
*ab = *a++ * *b; b += 3; /* r21 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r22 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r23 */
*ab += *a++ * *b; b += 3;
*ab += *a++ * *b; b = smxb->a;
ab++;
ar = a;
*ab = *a++ * *b; b += 3; /* r31 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r32 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r33 */
*ab += *a++ * *b; b += 3;
*ab += *a++ * *b++; bt = b;
ab++;
ar = smxa->a;
*ab = *ar++ * *b++; /* t1 */
*ab += *ar++ * *b++;
*ab += *ar++ * *b; b = bt;
*ab += *a++;
*ab %= STBF; if (*ab < 0) *ab += STBF;
ab++;
*ab = *ar++ * *b++; /* t2 */
*ab += *ar++ * *b++;
*ab += *ar++ * *b; b = bt;
*ab += *a++;
*ab %= STBF; if (*ab < 0) *ab += STBF;
ab++;
*ab = *ar++ * *b++; /* t3 */
*ab += *ar++ * *b++;
*ab += *ar * *b;
*ab += *a;
*ab %= STBF; if (*ab < 0) *ab += STBF;
}
void RTMxMultiply(T_RTMx *rtmxab, const T_RTMx *rtmxa, const T_RTMx *rtmxb,
int FacAug, int FacTr)
{
const int *ar, *a, *b, *bt;
int *ab;
/* no loops to be as fast as possible */
ar = rtmxa->a;
a = rtmxa->a;
b = rtmxb->a;
ab = rtmxab->a;
*ab = *a++ * *b; b += 3; /* r11 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r12 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r13 */
*ab += *a++ * *b; b += 3;
*ab += *a++ * *b; b = rtmxb->a;
ab++;
ar = a;
*ab = *a++ * *b; b += 3; /* r21 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r22 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r23 */
*ab += *a++ * *b; b += 3;
*ab += *a++ * *b; b = rtmxb->a;
ab++;
ar = a;
*ab = *a++ * *b; b += 3; /* r31 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r32 */
*ab += *a++ * *b; b += 3;
*ab += *a * *b; b -= 5;
ab++;
a = ar;
*ab = *a++ * *b; b += 3; /* r33 */
*ab += *a++ * *b; b += 3;
*ab += *a++ * *b++; bt = b;
ab++;
if (FacTr > 0)
{
ar = rtmxa->a;
*ab = *ar++ * *b++; /* t1 */
*ab += *ar++ * *b++;
*ab += *ar++ * *b; b = bt;
*ab += *a++ * FacAug;
*ab %= FacTr; if (*ab < 0) *ab += FacTr;
ab++;
*ab = *ar++ * *b++; /* t2 */
*ab += *ar++ * *b++;
*ab += *ar++ * *b; b = bt;
*ab += *a++ * FacAug;
*ab %= FacTr; if (*ab < 0) *ab += FacTr;
ab++;
*ab = *ar++ * *b++; /* t3 */
*ab += *ar++ * *b++;
*ab += *ar * *b;
*ab += *a * FacAug;
*ab %= FacTr; if (*ab < 0) *ab += FacTr;
}
else
{
ar = rtmxa->a;
*ab = *ar++ * *b++; /* t1 */
*ab += *ar++ * *b++;
*ab += *ar++ * *b; b = bt;
*ab += *a++ * FacAug;
ab++;
*ab = *ar++ * *b++; /* t2 */
*ab += *ar++ * *b++;
*ab += *ar++ * *b; b = bt;
*ab += *a++ * FacAug;
ab++;
*ab = *ar++ * *b++; /* t3 */
*ab += *ar++ * *b++;
*ab += *ar * *b;
*ab += *a * FacAug;
}
}
void InverseRotMx(const int *RotMx, int *InvRotMx)
{
InvRotMx[0] = RotMx[4] * RotMx[8] - RotMx[5] * RotMx[7];
InvRotMx[1] = - RotMx[1] * RotMx[8] + RotMx[2] * RotMx[7];
InvRotMx[2] = RotMx[1] * RotMx[5] - RotMx[2] * RotMx[4];
InvRotMx[3] = - RotMx[3] * RotMx[8] + RotMx[5] * RotMx[6];
InvRotMx[4] = RotMx[0] * RotMx[8] - RotMx[2] * RotMx[6];
InvRotMx[5] = - RotMx[0] * RotMx[5] + RotMx[2] * RotMx[3];
InvRotMx[6] = RotMx[3] * RotMx[7] - RotMx[4] * RotMx[6];
InvRotMx[7] = - RotMx[0] * RotMx[7] + RotMx[1] * RotMx[6];
InvRotMx[8] = RotMx[0] * RotMx[4] - RotMx[1] * RotMx[3];
}
void InverseRTMx(const T_RTMx *RTMx, T_RTMx *InvRTMx)
{
int *iR;
const int *T;
iR = InvRTMx->s.R;
InverseRotMx(RTMx->s.R, iR);
T = RTMx->s.T;
InvRTMx->s.T[0] = - iR[0] * T[0] - iR[1] * T[1] - iR[2] * T[2];
InvRTMx->s.T[1] = - iR[3] * T[0] - iR[4] * T[1] - iR[5] * T[2];
InvRTMx->s.T[2] = - iR[6] * T[0] - iR[7] * T[1] - iR[8] * T[2];
}
int IsSMxTransl0(const T_LatticeInfo *LatticeInfo, const int *SeitzMxT)
{
int iTrV, nTrV, t;
const int *TrV;
nTrV = LatticeInfo->nTrVector;
TrV = LatticeInfo->TrVector;
for (iTrV = 0; iTrV < nTrV; iTrV++)
{
t = (SeitzMxT[0] + TrV[0]) % STBF;
if (t == 0) {
t = (SeitzMxT[1] + TrV[1]) % STBF;
if (t == 0) {
t = (SeitzMxT[2] + TrV[2]) % STBF;
if (t == 0)
return 1;
}}
TrV += 3;
}
return 0;
}
static int IsSpecialSeitzMx(T_SgInfo *SgInfo, const T_RTMx *SMx, int ExpandLT)
{
int i, special, smx11;
const T_LatticeInfo *ExpLT;
/* check rotation part for identity or inversion operation */
smx11 = SMx->s.R[0];
if ( smx11 != 1
&& smx11 != -1) return 0;
for (i = 1; i < 9; i++)
{
if (i % 4) {
if (SMx->s.R[i] != 0) return 0; }
else {
if (SMx->s.R[i] != smx11) return 0; }
}
if (smx11 == 1) special = SpecialSMx_Identity;
else special = SpecialSMx_Inversion;
/* rotation part is identity or inversion
check translation part now
*/
if (IsSMxTransl0(SgInfo->LatticeInfo, SMx->s.T) == 1)
return (special | SpecialSMx_Transl0);
if (ExpandLT && (smx11 == 1 || SgInfo->Centric))
{
/* try to expand lattice type */
ExpLT = NULL;
switch (SgInfo->LatticeInfo->Code)
{
case 'P':
if (IsSMxTransl0(LI_A, SMx->s.T) == 1)
{ ExpLT = LI_A; break; }
if (IsSMxTransl0(LI_B, SMx->s.T) == 1)
{ ExpLT = LI_B; break; }
if (IsSMxTransl0(LI_C, SMx->s.T) == 1)
{ ExpLT = LI_C; break; }
if (IsSMxTransl0(LI_I, SMx->s.T) == 1)
{ ExpLT = LI_I; break; }
if (IsSMxTransl0(LI_R, SMx->s.T) == 1)
{ ExpLT = LI_R; break; }
if (IsSMxTransl0(LI_S, SMx->s.T) == 1)
{ ExpLT = LI_S; break; }
if (IsSMxTransl0(LI_T, SMx->s.T) == 1)
{ ExpLT = LI_T; break; }
case 'A':
case 'B':
case 'C':
if (IsSMxTransl0(LI_F, SMx->s.T) == 1)
ExpLT = LI_F;
}
if (ExpLT != NULL)
{
SgInfo->LatticeInfo = ExpLT;
SgInfo->StatusLatticeTr = -1;
return (special | SpecialSMx_Transl0);
}
}
return special;
}
int CompareSeitzMx(const T_LatticeInfo *LatticeInfo,
const T_RTMx *SeitzMxA, const T_RTMx *SeitzMxB)
{
int i, iTrV, nTrV, t;
const int *TrV;
/* compare rotation part */
for (i = 0; i < 9; i++)
if (SeitzMxA->s.R[i] != SeitzMxB->s.R[i]) return 1;
/* rotation part is same
check translation
*/
nTrV = LatticeInfo->nTrVector;
TrV = LatticeInfo->TrVector;
for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3)
{
t = (SeitzMxA->s.T[0] + TrV[0]) % STBF;
if (t == SeitzMxB->s.T[0]) {
t = (SeitzMxA->s.T[1] + TrV[1]) % STBF;
if (t == SeitzMxB->s.T[1]) {
t = (SeitzMxA->s.T[2] + TrV[2]) % STBF;
if (t == SeitzMxB->s.T[2])
return 0;
}}
}
return -1;
}
int GetRotMxOrder(const int *RotMx)
{
int deter = deterRotMx(RotMx);
if (deter == -1 || deter == 1)
{
switch (traceRotMx(RotMx))
{
case -3: return -1;
case -2: return -6;
case -1: if (deter == -1) return -4;
else return 2;
case 0: if (deter == -1) return -3;
else return 3;
case 1: if (deter == -1) return -2;
else return 4;
case 2: return 6;
case 3: return 1;
}
}
return 0;
}
static int nNextBasis_of_DirCode(const int DirCode,
const int **RMx, const int **InvRMx)
{
switch (DirCode)
{
case '.': *RMx = *InvRMx = NULL; return 1;
case '=':
case '"':
case '\'':
case '|':
case '\\': *RMx = RMx_3_111; *InvRMx = RMx_3i111; return 3;
case '*': *RMx = RMx_4_001; *InvRMx = RMx_4i001; return 4;
default:
break;
}
SetSgError("Internal Error: Corrupt DirCode");
return -1;
}
int GetRotMxInfo(const int *RotMx, T_RotMxInfo *RotMxInfo)
{
int i;
int nNextBasis, iNextBasis;
int nLoopInv, iLoopInv;
int Order, AbsOrder;
int RMxCopy[9], MatchMx[9], InvMatchMx[9], REV[3];
int *mmx;
const int *NBRMx, *InvNBRMx;
const T_TabXtalRotMx *txrmx;
Order = GetRotMxOrder(RotMx);
if (RotMxInfo)
RotMxInfo->Order = Order;
if (Order)
{
AbsOrder = abs(Order);
if (Order > 0)
for (i = 0; i < 9; i++) RMxCopy[i] = RotMx[i];
else
for (i = 0; i < 9; i++) RMxCopy[i] = -RotMx[i];
for (txrmx = TabXtalRotMx; txrmx->Order; txrmx++)
if (txrmx->Order == AbsOrder) break;
while (txrmx->Order == AbsOrder)
{
nNextBasis = nNextBasis_of_DirCode(txrmx->DirCode, &NBRMx, &InvNBRMx);
if (nNextBasis < 0)
return 0;
if (AbsOrder > 2) nLoopInv = 2;
else nLoopInv = 1;
for (i = 0; i < 9; i++) MatchMx[i] = txrmx->RMx[i];
for (iNextBasis = 0; iNextBasis < nNextBasis; iNextBasis++)
{
if (iNextBasis)
RotateRotMx(MatchMx, NBRMx, InvNBRMx);
mmx = MatchMx;
for (iLoopInv = 0; iLoopInv < nLoopInv; iLoopInv++)
{
if (iLoopInv)
{
InverseRotMx(MatchMx, InvMatchMx);
mmx = InvMatchMx;
}
for (i = 0; i < 9; i++)
if (mmx[i] != RMxCopy[i]) break;
if (i == 9) /* matrix has been found */
{
if (RotMxInfo)
{
RotMxInfo->Inverse = iLoopInv;
if (nNextBasis == 3)
{
switch(iNextBasis)
{
case 0: RotMxInfo->RefAxis = 'z'; break;
case 1: RotMxInfo->RefAxis = 'x'; break;
case 2: RotMxInfo->RefAxis = 'y'; break;
}
}
else
RotMxInfo->RefAxis = 'o';
RotMxInfo->DirCode = txrmx->DirCode;
for (i = 0; i < 3; i++)
RotMxInfo->EigenVector[i] = txrmx->EigenVector[i];
for (; iNextBasis--;)
{
RotMx_t_Vector(REV, NBRMx, RotMxInfo->EigenVector, 0);
if (iNextBasis-- == 0)
{
for (i = 0; i < 3; i++)
RotMxInfo->EigenVector[i] = REV[i];
break;
}
RotMx_t_Vector(RotMxInfo->EigenVector, NBRMx, REV, 0);
}
}
return Order;
}
}
}
txrmx++;
}
}
return 0;
}
const T_RotMxInfo *ListOrBufRotMxInfo(const T_SgInfo *SgInfo, int iList,
T_RotMxInfo *BufRotMxInfo)
{
T_RotMxInfo *RMxI;
RMxI = SgInfo->ListRotMxInfo;
if (RMxI)
RMxI += iList;
else
{
RMxI = BufRotMxInfo;
if (GetRotMxInfo(SgInfo->ListSeitzMx[iList].s.R, RMxI) == 0) {
SetSgError(Err_Ill_SMx_in_List);
return NULL;
}
}
return RMxI;
}
static int CoreAdd2ListSeitzMx(T_SgInfo *SgInfo, const T_RTMx *NewSMx)
{
int i, iList;
T_RTMx *lsmx;
T_RotMxInfo RotMxInfo;
const T_LatticeInfo *LI;
static const char *Err_NonXtalOp =
"Error: Generators produce non-crystallographic operation";
if (SgInfo->GenOption) LI = SgInfo->LatticeInfo;
else LI = LI_P;
lsmx = SgInfo->ListSeitzMx;
for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
if (CompareSeitzMx(LI, NewSMx, lsmx) == 0)
return 0; /* matrix is not unique */
if (GetRotMxInfo(NewSMx->s.R, &RotMxInfo) == 0) {
SetSgError(Err_NonXtalOp);
return -1;
}
if (SgInfo->nList >= SgInfo->MaxList)
{
if (SgInfo->nList >= 192)
SetSgError(Err_NonXtalOp);
else
SetSgError("Internal Error: Allocated space for ListSeitzMx too small");
return -1;
}
for (i = 0; i < 12; i++) lsmx->a[i] = NewSMx->a[i];
if (SgInfo->ListRotMxInfo != NULL)
CopyRotMxInfo(&SgInfo->ListRotMxInfo[SgInfo->nList], &RotMxInfo);
SgInfo->nList++;
return 1;
}
int Add2ListSeitzMx(T_SgInfo *SgInfo, const T_RTMx *NewSMx)
{
int i, special, Enter;
int iMult, jMult;
T_RTMx TrialSMx, *iMultSMx, *jMultSMx;
static const char *Err_Ill_lattice_translation =
"Error: Illegal lattice translation";
if (SgInfo->nList == 0)
{
/* make sure identity matrix is first in list */
InitSeitzMx(&TrialSMx, 1);
if (CoreAdd2ListSeitzMx(SgInfo, &TrialSMx) < 0)
return -1;;
}
for (i = 0; i < 9; i++)
TrialSMx.s.R[i] = NewSMx->s.R[i];
for (i = 0; i < 3; i++) {
TrialSMx.s.T[i] = NewSMx->s.T[i] % STBF;
if (TrialSMx.s.T[i] < 0)
TrialSMx.s.T[i] += STBF;
}
iMult = SgInfo->nList;
iMultSMx = &SgInfo->ListSeitzMx[iMult];
jMult = 1;
jMultSMx = &SgInfo->ListSeitzMx[1]; /* skip first = identity matrix */
for (;;)
{
Enter = 1;
special = IsSpecialSeitzMx(SgInfo, &TrialSMx, 1);
if (special & SpecialSMx_Identity)
{
if (! (special & SpecialSMx_Transl0)) {
SetSgError(Err_Ill_lattice_translation);
return -1;
}
if (SgInfo->GenOption)
Enter = 0;
}
else if (special & SpecialSMx_Inversion)
{
if (! (special & SpecialSMx_Transl0))
{
if (SgInfo->Centric) {
SetSgError(Err_Ill_lattice_translation);
return -1;
}
SgInfo->InversionOffOrigin = 1;
}
else
{
if (SgInfo->InversionOffOrigin)
SgInfo->Centric = 1;
SgInfo->InversionOffOrigin = 0;
if (SgInfo->GenOption)
{
if (SgInfo->Centric == 0)
SgInfo->Centric = -1;
Enter = 0;
}
else
SgInfo->Centric = 1;
}
}
if (Enter && CoreAdd2ListSeitzMx(SgInfo, &TrialSMx) < 0)
return -1;
if (SgInfo->GenOption < 0)
break;
if (jMult > iMult)
{
iMult++;
iMultSMx++;
jMult = 1;
jMultSMx = &SgInfo->ListSeitzMx[1]; /* skip first = identity matrix */
}
if (iMult == SgInfo->nList)
break;
SeitzMxMultiply(&TrialSMx, jMultSMx, iMultSMx);
jMult++;
jMultSMx++;
}
return 0;
}
int AddInversion2ListSeitzMx(T_SgInfo *SgInfo)
{
T_RTMx SeitzMx;
InitSeitzMx(&SeitzMx, -1);
return Add2ListSeitzMx(SgInfo, &SeitzMx);
}
int AddLatticeTr2ListSeitzMx(T_SgInfo *SgInfo,
const T_LatticeInfo *LatticeInfo)
{
int iTrV, nTrV;
const int *TrV;
T_RTMx SeitzMx;
InitRotMx(SeitzMx.s.R, 1);
nTrV = LatticeInfo->nTrVector;
TrV = &LatticeInfo->TrVector[3]; /* skip first vector which is always 000 */
for (iTrV = 1; iTrV < nTrV; iTrV++)
{
SeitzMx.s.T[0] = *TrV++;
SeitzMx.s.T[1] = *TrV++;
SeitzMx.s.T[2] = *TrV++;
if (Add2ListSeitzMx(SgInfo, &SeitzMx) < 0)
return -1;
}
if (SgInfo->GenOption)
SgInfo->StatusLatticeTr = 0;
else
SgInfo->StatusLatticeTr = 1;
return 0;
}
static int RemoveLatticeTr(T_SgInfo *SgInfo)
{
int iList, jList, i;
T_RTMx *smxi, *smxj, *lastsmx;
T_RotMxInfo *lrmxi, *lastrmxi;
if (SgInfo->LatticeInfo->Code == 'P')
return 0;
if (SgInfo->StatusLatticeTr == -1)
{
if (AddLatticeTr2ListSeitzMx(SgInfo, SgInfo->LatticeInfo) < 0)
return -1;
}
iList = 0;
while (iList < SgInfo->nList)
{
smxi = &SgInfo->ListSeitzMx[iList];
jList = iList + 1;
while (jList < SgInfo->nList)
{
smxj = &SgInfo->ListSeitzMx[jList];
if (CompareSeitzMx(SgInfo->LatticeInfo, smxi, smxj) == 0)
{
/* copy last element to this place */
SgInfo->nList--;
lastsmx = &SgInfo->ListSeitzMx[SgInfo->nList];
for (i = 0; i < 12; i++) smxj->a[i] = lastsmx->a[i];
if (SgInfo->ListRotMxInfo != NULL)
{
lrmxi = &SgInfo->ListRotMxInfo[jList];
lastrmxi = &SgInfo->ListRotMxInfo[SgInfo->nList];
CopyRotMxInfo(lrmxi, lastrmxi);
}
}
else
jList++;
}
iList++;
}
SgInfo->StatusLatticeTr = 0;
return 0;
}
static int RemoveInversion(T_SgInfo *SgInfo)
{
int i, iList, special, deter, Centric, InversionOffOrigin;
T_RTMx *lsmx, *smx, ProperSMx;
T_RotMxInfo *lrmxi, *rmxi;
Centric = 0;
InversionOffOrigin = 0;
lsmx = SgInfo->ListSeitzMx;
for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
{
special = IsSpecialSeitzMx(SgInfo, lsmx, 0);
if (special & SpecialSMx_Inversion)
{
if (special & SpecialSMx_Transl0)
Centric = 1;
else
InversionOffOrigin = 1;
break;
}
}
if (InversionOffOrigin && Centric) {
SetSgError("Internal Error: Corrupt SgInfo");
return -1;
}
if (Centric == 0)
{
if (InversionOffOrigin) {
SgInfo->Centric = 0;
SgInfo->InversionOffOrigin = 1;
}
else
{
if (SgInfo->Centric) SgInfo->Centric = -1;
SgInfo->InversionOffOrigin = 0;
}
}
else
{
SgInfo->InversionOffOrigin = 0;
lsmx = SgInfo->ListSeitzMx;
lrmxi = SgInfo->ListRotMxInfo;
iList = 0;
while (iList < SgInfo->nList)
{
deter = deterRotMx(lsmx->s.R);
if (deter == -1 && SgInfo->Centric == -1)
{
for (i = 0; i < 9; i++)
ProperSMx.s.R[i] = -lsmx->s.R[i];
for (i = 0; i < 3; i++) {
ProperSMx.s.T[i] = -lsmx->s.T[i] % STBF;
if (ProperSMx.s.T[i] < 0)
ProperSMx.s.T[i] += STBF;
}
smx = SgInfo->ListSeitzMx;
for (i = 0; i < SgInfo->nList; i++, smx++)
if (CompareSeitzMx(LI_P, &ProperSMx, smx) == 0) break;
if (i == SgInfo->nList)
{
for (i = 0; i < 12; i++) lsmx->a[i] = ProperSMx.a[i];
deter = deterRotMx(lsmx->s.R);
if (deter != 1 || (lrmxi && GetRotMxInfo(lsmx->s.R, lrmxi) == 0)) {
SetSgError(Err_Ill_SMx_in_List);
return -1;
}
}
}
if (deter == -1)
{
/* copy last element to this place */
SgInfo->nList--;
if (SgInfo->nList != iList)
{
smx = &SgInfo->ListSeitzMx[SgInfo->nList];
for (i = 0; i < 12; i++) lsmx->a[i] = smx->a[i];
if (lrmxi)
{
rmxi = &SgInfo->ListRotMxInfo[SgInfo->nList];
CopyRotMxInfo(lrmxi, rmxi);
}
}
}
else if (deter == 1)
{
lsmx++;
if (lrmxi != NULL) lrmxi++;
iList++;
}
else {
SetSgError(Err_Ill_SMx_in_List);
return -1;
}
}
SgInfo->Centric = -1;
}
return 0;
}
int ApplyOriginShift(T_SgInfo *SgInfo)
{
int HaveOrSh, iList, i;
T_RTMx *lsmx, SMx;
int OrSh[3], lo[3];
HaveOrSh = 0;
for (i = 0; i < 3; i++) {
OrSh[i] = SgInfo->OriginShift[i] * (STBF / 12);
if (OrSh[i]) HaveOrSh = 1;
}
if (HaveOrSh == 0)
return 0;
lsmx = SgInfo->ListSeitzMx;
for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
{
RotMx_t_Vector(lo, lsmx->s.R, OrSh, STBF);
for (i = 0; i < 3; i++)
lsmx->s.T[i] = iModPositive(lsmx->s.T[i] - lo[i] + OrSh[i], STBF);
}
if (SgInfo->Centric == -1)
{
lsmx = &SMx;
InitSeitzMx(lsmx, -1);
RotMx_t_Vector(lo, lsmx->s.R, OrSh, STBF);
for (i = 0; i < 3; i++)
lsmx->s.T[i] = iModPositive(lsmx->s.T[i] - lo[i] + OrSh[i], STBF);
if (CoreAdd2ListSeitzMx(SgInfo, lsmx) < 0)
return -1;
SgInfo->Centric = 0;
SgInfo->InversionOffOrigin = 1;
}
return 1;
}
static void TidyTranslation(T_SgInfo *SgInfo)
{
int iList;
int iTrV, nTrV;
const int *TrV;
T_RTMx *lsmx;
int t1, t2, t3, *lt1, *lt2, *lt3, mint1, mint2, mint3;
int n0t, n0mint;
nTrV = SgInfo->LatticeInfo->nTrVector;
lsmx = SgInfo->ListSeitzMx;
for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
{
mint1 = *(lt1 = &lsmx->s.T[0]);
mint2 = *(lt2 = &lsmx->s.T[1]);
mint3 = *(lt3 = &lsmx->s.T[2]);
TrV = SgInfo->LatticeInfo->TrVector;
for (iTrV = 0; iTrV < nTrV; iTrV++)
{
t1 = ((*lt1) + *TrV++) % STBF;
t2 = ((*lt2) + *TrV++) % STBF;
t3 = ((*lt3) + *TrV++) % STBF;
n0t = (t1 == 0) + (t2 == 0) + (t3 == 0);
n0mint = (mint1 == 0) + (mint2 == 0) + (mint3 == 0);
if ( n0t > n0mint
|| ( n0t == n0mint
&& ( t1 + t2 + t3 < mint1 + mint2 + mint3
|| ( t1 + t2 + t3 == mint1 + mint2 + mint3
&& (t1 < mint1 || (t1 == mint1 && t2 < mint2)))))) {
mint1 = t1; mint2 = t2; mint3 = t3;
}
}
*lt1 = mint1;
*lt2 = mint2;
*lt3 = mint3;
}
}
static T_SgInfo *Pt_SgInfo_ListSortFunction = NULL;
static int SgInfoListSortFunction(const int *iList_a, const int *iList_b)
{
int val_a, val_b, n0_a, n0_b, i;
T_RTMx *smx_a, *smx_b;
T_RotMxInfo RotMxInfo_a, RotMxInfo_b, *rmxi_a, *rmxi_b;
T_SgInfo *SgInfo = Pt_SgInfo_ListSortFunction;
if (SgError != NULL) return 0;
if (SgInfo->ListRotMxInfo == NULL)
{
rmxi_a = &RotMxInfo_a;
rmxi_b = &RotMxInfo_b;
smx_a = &SgInfo->ListSeitzMx[*iList_a];
smx_b = &SgInfo->ListSeitzMx[*iList_b];
if ( GetRotMxInfo(smx_a->s.R, rmxi_a) == 0
|| GetRotMxInfo(smx_b->s.R, rmxi_b) == 0)
{
SetSgError(Err_Ill_SMx_in_List);
return 0;
}
}
else
{
rmxi_a = &SgInfo->ListRotMxInfo[*iList_a];
rmxi_b = &SgInfo->ListRotMxInfo[*iList_b];
}
val_a = abs(rmxi_a->Order);
val_b = abs(rmxi_b->Order);
if (val_a == 1 && val_b != 1) return -1;
if (val_a != 1 && val_b == 1) return 1;
if (rmxi_a->Order == 1 && rmxi_b->Order != 1) return -1;
if (rmxi_a->Order != 1 && rmxi_b->Order == 1) return 1;
if (val_a != 1)
{
if (val_a > val_b) return -1;
if (val_a < val_b) return 1;
if (rmxi_a->Order > rmxi_b->Order) return -1;
if (rmxi_a->Order < rmxi_b->Order) return 1;
}
n0_a = n0_b = 0;
for (i = 0; i < 3; i++)
{
if (rmxi_a->EigenVector[i] == 0) n0_a++;
if (rmxi_b->EigenVector[i] == 0) n0_b++;
}
if (n0_a > n0_b) return -1;
if (n0_a < n0_b) return 1;
val_a = val_b = 0;
for (i = 0; i < 3; i++)
{
if (val_a < abs(rmxi_a->EigenVector[i]))
val_a = abs(rmxi_a->EigenVector[i]);
if (val_b < abs(rmxi_b->EigenVector[i]))
val_b = abs(rmxi_b->EigenVector[i]);
}
if (val_a < val_b) return -1;
if (val_a > val_b) return 1;
val_a = 100 * abs(rmxi_a->EigenVector[2]);
val_a += 10 * abs(rmxi_a->EigenVector[0]);
val_a += abs(rmxi_a->EigenVector[1]);
val_b = 100 * abs(rmxi_b->EigenVector[2]);
val_b += 10 * abs(rmxi_b->EigenVector[0]);
val_b += abs(rmxi_b->EigenVector[1]);
if (n0_a < 2)
{
if (val_a < val_b) return -1;
if (val_a > val_b) return 1;
}
else
{
if (val_a > val_b) return -1;
if (val_a < val_b) return 1;
}
for (i = 0; i < 3; i++)
{
if (rmxi_a->EigenVector[i] > rmxi_b->EigenVector[i]) return -1;
if (rmxi_a->EigenVector[i] < rmxi_b->EigenVector[i]) return 1;
}
if (rmxi_a->Inverse < rmxi_b->Inverse) return -1;
if (rmxi_a->Inverse > rmxi_b->Inverse) return 1;
smx_a = &SgInfo->ListSeitzMx[*iList_a];
smx_b = &SgInfo->ListSeitzMx[*iList_b];
for (i = 0; i < 3; i++)
{
if (smx_a->s.T[i] < smx_b->s.T[i]) return -1;
if (smx_a->s.T[i] > smx_b->s.T[i]) return 1;
}
return 0;
}
static int PostSortSgInfoList(const T_SgInfo *SgInfo, int *List_iList)
{
int nList, iL_iL, jL_iL;
T_RTMx BufMx1, BufMx2, *smxab;
const T_RTMx *lsmx, *smxb;
T_RotMxInfo RotMxInfo;
const T_RotMxInfo *rmxi;
int nO_1, iO, save, i, i_;
nList = SgInfo->nList;
iL_iL = 0;
while (iL_iL < nList)
{
lsmx = &SgInfo->ListSeitzMx[List_iList[iL_iL]];
rmxi = ListOrBufRotMxInfo(SgInfo, List_iList[iL_iL], &RotMxInfo);
if (rmxi == NULL)
return -1;
iL_iL++;
iO = rmxi->Order;
if (iO < 0 && iO % 2) iO *= 2;
nO_1 = abs(iO) - 1;
smxab = &BufMx2;
smxb = lsmx;
for (iO = 1; iO < nO_1; iO++)
{
SeitzMxMultiply(smxab, lsmx, smxb);
for (jL_iL = iL_iL; jL_iL < nList; jL_iL++)
{
smxb = &SgInfo->ListSeitzMx[List_iList[jL_iL]];
if (CompareSeitzMx(SgInfo->LatticeInfo, smxab, smxb) == 0)
break;
}
if (jL_iL < nList)
{
save = List_iList[jL_iL];
for (i = i_ = jL_iL; i > iL_iL; i = i_)
List_iList[i] = List_iList[--i_];
List_iList[iL_iL++] = save;
}
smxb = smxab;
if (iO % 2) smxab = &BufMx1;
else smxab = &BufMx2;
}
}
return 0;
}
static void SortSgInfoList(T_SgInfo *SgInfo, int *List_iList)
{
int i, j, refi;
int nList;
T_RTMx *lsmx;
T_RotMxInfo *lrmxi;
if (SgError != NULL) return;
nList = SgInfo->nList;
lsmx = SgInfo->ListSeitzMx;
lrmxi = SgInfo->ListRotMxInfo;
Pt_SgInfo_ListSortFunction = SgInfo;
for (i = 0; i < nList; i++) List_iList[i] = i;
qsort((void *) List_iList, nList, sizeof (*List_iList),
(int (*)(const void *, const void *)) SgInfoListSortFunction);
Pt_SgInfo_ListSortFunction = NULL;
if (SgError != NULL) return;
if (PostSortSgInfoList(SgInfo, List_iList) != 0)
return;
for (i = 0; i < nList; i++)
{
j = List_iList[i];
if (j != i)
{
for (refi = i + 1; refi < nList; refi++)
if (List_iList[refi] == i) break;
if (refi >= nList) {
SetSgError("Internal Error: SortSgInfoList(): Corrupt List_iList");
return;
}
SwapRTMx(&lsmx[i], &lsmx[j]);
if (lrmxi != NULL)
SwapRotMxInfo(&lrmxi[i], &lrmxi[j]);
List_iList[refi] = j;
}
}
}
int FindSeitzMx(const T_SgInfo *SgInfo,
int Order, int HonorSign, int RefAxis, int DirCode)
{
int iList;
int MatchingOrder;
T_RTMx *lsmx;
T_RotMxInfo *lrmxi, RotMxInfo;
lsmx = SgInfo->ListSeitzMx;
lrmxi = SgInfo->ListRotMxInfo;
if (lrmxi == NULL) lrmxi = &RotMxInfo;
for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
{
if (lrmxi == &RotMxInfo)
{
if (GetRotMxInfo(lsmx->s.R, lrmxi) == 0) {
SetSgError(Err_Ill_SMx_in_List);
return -1;
}
}
if (HonorSign == 0)
MatchingOrder = (abs(Order) == abs(lrmxi->Order));
else
MatchingOrder = (Order == lrmxi->Order);
if ( MatchingOrder
&& lrmxi->Inverse == 0
&& (RefAxis == 0 || RefAxis == lrmxi->RefAxis)
&& (DirCode == 0 || DirCode == lrmxi->DirCode))
{
if (DirCode != '*') return iList;
if ( lrmxi->EigenVector[0] == 1
&& lrmxi->EigenVector[1] == 1
&& lrmxi->EigenVector[2] == 1) return iList;
}
if (lrmxi != &RotMxInfo) lrmxi++;
}
return -1;
}
static int FindXtalSystem(T_SgInfo *SgInfo)
{
int iList, i;
int HonorSign = 0, CheckEnantiomorph;
const T_RTMx *lsmx;
T_RotMxInfo RotMxInfo;
const T_RotMxInfo *lrmxi;
enum Axis { N1 = 0, N2, N3, N4, N6, EndOfAxis };
int N_count[EndOfAxis];
SgInfo->XtalSystem = XS_Unknown;
SgInfo->UniqueRefAxis = 0;
SgInfo->UniqueDirCode = 0;
SgInfo->ExtraInfo = EI_Unknown;
CheckEnantiomorph = 0;
for (i = 0; i < EndOfAxis; i++) N_count[i] = 0;
lsmx = SgInfo->ListSeitzMx;
lrmxi = SgInfo->ListRotMxInfo;
if (lrmxi == NULL) lrmxi = &RotMxInfo;
for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
{
if (lrmxi == &RotMxInfo)
{
if (GetRotMxInfo(lsmx->s.R, &RotMxInfo) == 0) {
SetSgError(Err_Ill_SMx_in_List);
return XS_Unknown;
}
}
switch(abs(lrmxi->Order))
{
case 1: i = N1; break;
case 2: i = N2; break;
case 3: i = N3; break;
case 4: i = N4; break;
case 6: i = N6; break;
default:
SetSgError("Internal Error: FindXtalSystem(): Corrupt ListRotMxInfo");
return XS_Unknown;
}
if (lrmxi->Inverse == 0) /* skip N^-1 */
N_count[i]++;
if (lrmxi != &RotMxInfo) lrmxi++;
}
i = EndOfAxis;
if (SgInfo->InversionOffOrigin == 1)
{
for (i = 0; i < EndOfAxis; i++)
{
if (N_count[i] % 2) break;
N_count[i] /= 2;
}
}
if (i == EndOfAxis)
{
if (N_count[N3] == 4) SgInfo->XtalSystem = XS_Cubic;
else if (N_count[N3] > 1) SgInfo->XtalSystem = XS_Unknown;
else if (N_count[N6] == 1) SgInfo->XtalSystem = XS_Hexagonal;
else if (N_count[N6] > 0) SgInfo->XtalSystem = XS_Unknown;
else if (N_count[N3] == 1) SgInfo->XtalSystem = XS_Trigonal;
else if (N_count[N4] == 1) SgInfo->XtalSystem = XS_Tetragonal;
else if (N_count[N4] > 0) SgInfo->XtalSystem = XS_Unknown;
else if (N_count[N2] > 2) SgInfo->XtalSystem = XS_Orthorhombic;
else if (N_count[N2] > 0) SgInfo->XtalSystem = XS_Monoclinic;
else if (N_count[N1] > 0) SgInfo->XtalSystem = XS_Triclinic;
HonorSign = 1;
iList = FindSeitzMx(SgInfo, -1, HonorSign, 'o', '.');
if (iList < 0) HonorSign = 0;
switch(SgInfo->XtalSystem)
{
case XS_Monoclinic:
iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '=');
if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
break;
case XS_Tetragonal:
CheckEnantiomorph = 1;
iList = FindSeitzMx(SgInfo, 4, HonorSign, 0, '=');
if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
break;
case XS_Trigonal:
CheckEnantiomorph = 1;
iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '=');
if (iList < 0)
iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '*');
if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
break;
case XS_Hexagonal:
CheckEnantiomorph = 1;
iList = FindSeitzMx(SgInfo, 6, HonorSign, 0, '=');
if (iList < 0) SgInfo->XtalSystem = XS_Unknown;
break;
case XS_Cubic:
iList = FindSeitzMx(SgInfo, 4, HonorSign, 0, '=');
if (iList >= 0) CheckEnantiomorph = 1;
break;
default:
iList = -1;
break;
}
}
if (SgInfo->XtalSystem == XS_Unknown) {
SetSgError("Error: Can't determine crystal system");
}
else if (iList >= 0)
{
lrmxi = ListOrBufRotMxInfo(SgInfo, iList, &RotMxInfo);
if (lrmxi == NULL) {
SgInfo->XtalSystem = XS_Unknown;
return XS_Unknown;
}
if (SgInfo->XtalSystem != XS_Cubic)
{
SgInfo->UniqueRefAxis = lrmxi->RefAxis;
SgInfo->UniqueDirCode = lrmxi->DirCode;
if (SgInfo->XtalSystem == XS_Trigonal && lrmxi->DirCode == '=')
{
switch (lrmxi->RefAxis)
{
case 'z':
switch (SgInfo->LatticeInfo->Code)
{
case 'R': SgInfo->ExtraInfo = EI_Obverse; break;
case 'T': SgInfo->ExtraInfo = EI_Reverse; break;
}
break;
case 'y':
switch (SgInfo->LatticeInfo->Code)
{
case 'S': SgInfo->ExtraInfo = EI_Obverse; break;
case 'R': SgInfo->ExtraInfo = EI_Reverse; break;
}
break;
case 'x':
switch (SgInfo->LatticeInfo->Code)
{
case 'T': SgInfo->ExtraInfo = EI_Obverse; break;
case 'S': SgInfo->ExtraInfo = EI_Reverse; break;
}
break;
}
}
}
if ( HonorSign == 0 /* no inversion matrix */
&& SgInfo->LatticeInfo->Code == 'P'
&& CheckEnantiomorph == 1)
{
lsmx = &SgInfo->ListSeitzMx[iList];
if (GetRotMxOrder(lsmx->s.R) > 1)
{
i = lsmx->s.T[0] * lrmxi->EigenVector[0];
i += lsmx->s.T[1] * lrmxi->EigenVector[1];
i += lsmx->s.T[2] * lrmxi->EigenVector[2];
if (i % (STBF / 2)) SgInfo->ExtraInfo = EI_Enantiomorphic;
}
}
}
return SgInfo->XtalSystem;
}
static int BuildGenerator_iList(T_SgInfo *SgInfo)
{
int iList, iList_1, nG;
int SgInfo_CI, HonorSign, Flag3asterisk, FlagPG;
SgInfo_CI = (SgInfo->Centric || SgInfo->InversionOffOrigin);
SgInfo->PointGroup = PG_Unknown;
nG = SgInfo->nGenerator = 0;
#define G_iL SgInfo->Generator_iList
HonorSign = 1;
iList_1 = FindSeitzMx(SgInfo, -1, HonorSign, 'o', '.');
if (iList_1 < 0) HonorSign = 0;
if (SgInfo->XtalSystem == XS_Unknown)
FindXtalSystem(SgInfo);
switch(SgInfo->XtalSystem)
{
case XS_Triclinic:
if (iList_1 < 0)
iList_1 = FindSeitzMx(SgInfo, 1, HonorSign, 'o', '.');
if (iList_1 >= 0) G_iL[nG++] = iList_1;
if (SgInfo_CI)
SgInfo->PointGroup = PG_1b;
else
SgInfo->PointGroup = PG_1;
SgInfo->nGenerator = nG;
return 0;
case XS_Monoclinic:
iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '=');
if (iList < 0) break;
G_iL[nG++] = iList;
if (SgInfo_CI)
SgInfo->PointGroup = PG_2_m;
else if (deterRotMx(SgInfo->ListSeitzMx[iList].s.R) == -1)
SgInfo->PointGroup = PG_m;
else
SgInfo->PointGroup = PG_2;
if (iList_1 >= 0) G_iL[nG++] = iList_1;
SgInfo->nGenerator = nG;
return 0;
case XS_Orthorhombic:
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'z', '=');
if (iList >= 0) G_iL[nG++] = iList;
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'x', '=');
if (iList >= 0) G_iL[nG++] = iList;
if (nG < 2)
{
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'y', '=');
if (iList >= 0) G_iL[nG++] = iList;
}
if (nG != 2) break;
if (SgInfo_CI)
SgInfo->PointGroup = PG_mmm;
else if ( deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1
|| deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
SgInfo->PointGroup = PG_mm2;
else
SgInfo->PointGroup = PG_222;
if (iList_1 >= 0) G_iL[nG++] = iList_1;
SgInfo->nGenerator = nG;
return 0;
case XS_Tetragonal:
iList = FindSeitzMx(SgInfo, 4, HonorSign, 0, '=');
if (iList < 0) break;
G_iL[nG++] = iList;
if ( SgInfo->UniqueRefAxis != 'x')
{
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'x', '=');
if (iList >= 0) G_iL[nG++] = iList;
}
if (nG < 2 && SgInfo->UniqueRefAxis != 'z')
{
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'z', '=');
if (iList >= 0) G_iL[nG++] = iList;
}
if (nG < 2 && SgInfo->UniqueRefAxis != 'y')
{
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'y', '=');
if (iList >= 0) G_iL[nG++] = iList;
}
if (nG < 2)
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_4_m;
else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
SgInfo->PointGroup = PG_4b;
else
SgInfo->PointGroup = PG_4;
}
else
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_4_mmm;
else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
{
if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
SgInfo->PointGroup = PG_4bm2;
else
SgInfo->PointGroup = PG_4b2m;
}
else
{
if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
SgInfo->PointGroup = PG_4mm;
else
SgInfo->PointGroup = PG_422;
}
}
if (iList_1 >= 0) G_iL[nG++] = iList_1;
SgInfo->nGenerator = nG;
return 0;
case XS_Trigonal:
case XS_Hexagonal:
Flag3asterisk = 0;
if (SgInfo->XtalSystem == XS_Trigonal)
{
iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '=');
if (iList < 0)
{
iList = FindSeitzMx(SgInfo, 3, HonorSign, 0, '*');
Flag3asterisk = 1;
}
}
else
iList = FindSeitzMx(SgInfo, 6, HonorSign, 0, '=');
if (iList < 0) break;
G_iL[nG++] = iList;
iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '\'');
if (iList >= 0) G_iL[nG++] = iList;
FlagPG = -1;
if (nG < 2)
{
iList = FindSeitzMx(SgInfo, 2, HonorSign, 0, '"');
if (iList >= 0) G_iL[nG++] = iList;
FlagPG = 1;
}
if (SgInfo->XtalSystem == XS_Trigonal)
{
if (nG < 2)
{
if (SgInfo_CI) SgInfo->PointGroup = PG_3b;
else SgInfo->PointGroup = PG_3;
}
else
{
if (Flag3asterisk == 1)
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_3bm;
else
{ FlagPG = deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R);
if (FlagPG == -1) SgInfo->PointGroup = PG_3m;
else SgInfo->PointGroup = PG_32;
}
}
else if (FlagPG == -1)
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_3b1m;
else
{ FlagPG = deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R);
if (FlagPG == -1) SgInfo->PointGroup = PG_31m;
else SgInfo->PointGroup = PG_312;
}
}
else
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_3bm1;
else
{ FlagPG = deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R);
if (FlagPG == -1) SgInfo->PointGroup = PG_3m1;
else SgInfo->PointGroup = PG_321;
}
}
}
}
else
{
if (nG < 2)
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_6_m;
else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
SgInfo->PointGroup = PG_6b;
else
SgInfo->PointGroup = PG_6;
}
else
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_6_mmm;
else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
{
if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == FlagPG)
SgInfo->PointGroup = PG_6b2m;
else
SgInfo->PointGroup = PG_6bm2;
}
else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[1]].s.R) == -1)
SgInfo->PointGroup = PG_6mm;
else
SgInfo->PointGroup = PG_622;
}
}
if (iList_1 >= 0) G_iL[nG++] = iList_1;
SgInfo->nGenerator = nG;
return 0;
case XS_Cubic:
FlagPG = 0;
iList = FindSeitzMx(SgInfo, 4, HonorSign, 'z', '=');
if (iList < 0) {
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'z', '=');
FlagPG = 1;
}
if (iList < 0) break;
G_iL[nG++] = iList;
iList = FindSeitzMx(SgInfo, 2, HonorSign, 'x', '=');
if (iList < 0) break;
G_iL[nG++] = iList;
iList = FindSeitzMx(SgInfo, 3, HonorSign, 'o', '*');
if (iList < 0) break;
G_iL[nG++] = iList;
if (FlagPG)
{
if (SgInfo_CI) SgInfo->PointGroup = PG_m3b;
else SgInfo->PointGroup = PG_23;
}
else
{
if (SgInfo_CI)
SgInfo->PointGroup = PG_m3bm;
else if (deterRotMx(SgInfo->ListSeitzMx[G_iL[0]].s.R) == -1)
SgInfo->PointGroup = PG_4b3m;
else
SgInfo->PointGroup = PG_432;
}
if (iList_1 >= 0) G_iL[nG++] = iList_1;
SgInfo->nGenerator = nG;
return 0;
default:
break;
}
#undef G_iL
return -1;
}
static int BuildHSym(T_SgInfo *SgInfo)
{
int NeedDash, HaveOrSh, nGRT, iList, iG, ip, os, i;
int AbsOrder, RefAxis, DirCode, ScrewPrimitive, Screw;
int PreviousRotation, PreviousRefAxis, PreviousDirCode;
int nTrV, iTrV, OrSh[3], RO[3], Transl[3];
const int *TrV, *ht;
T_RTMx SMx_1;
const T_RTMx *SMx;
const T_RotMxInfo *rmxi;
char *hsym, *hsym_mark;
struct
{
T_RotMxInfo RMxI_Buf;
const T_RotMxInfo *RMxI;
int Transl[3];
}
GRT[sizeof SgInfo->Generator_iList / sizeof (*SgInfo->Generator_iList) + 1];
const char *Digits = "0123456";
if (SgInfo->nGenerator == 0) {
SetSgError("Internal Error: BuildHSym(): Empty generator list");
return -1;
}
HaveOrSh = 0;
for (i = 0; i < 3; i++) {
OrSh[i] = SgInfo->OriginShift[i] * (STBF / 12);
if (OrSh[i]) HaveOrSh = 1;
}
NeedDash = 0;
nGRT = 0;
for (iG = 0; iG < SgInfo->nGenerator; iG++)
{
iList = SgInfo->Generator_iList[iG];
GRT[nGRT].RMxI = ListOrBufRotMxInfo(SgInfo, iList, &GRT[nGRT].RMxI_Buf);
if (GRT[nGRT].RMxI == NULL)
return -1;
SMx = &SgInfo->ListSeitzMx[iList];
RotMx_t_Vector(RO, SMx->s.R, OrSh, STBF);
for (i = 0; i < 3; i++)
GRT[nGRT].Transl[i] = iModPositive(SMx->s.T[i] + RO[i] - OrSh[i], STBF);
if (GRT[nGRT].RMxI->Order == -1)
{
for (i = 0; i < 3; i++)
if (GRT[nGRT].Transl[i] != 0) break;
if (i == 3) NeedDash = 1;
else nGRT++;
}
else
nGRT++;
}
if (SgInfo->Centric)
{
if (HaveOrSh == 0)
NeedDash = 1;
else
{
for (iG = 0; iG < nGRT; iG++)
if (GRT[iG].RMxI->Order == 1) break;
InitSeitzMx(&SMx_1, -1);
if (GetRotMxInfo(SMx_1.s.R, &GRT[iG].RMxI_Buf) != -1) {
SetSgError("Internal Error: BuildHSym(): Corrupt GetRotMxInfo()");
return -1;
}
GRT[iG].RMxI = &GRT[iG].RMxI_Buf;
for (i = 0; i < 3; i++)
GRT[iG].Transl[i] = iModPositive(-2 * OrSh[i], STBF);
if (iG == nGRT)
nGRT++;
}
}
hsym = SgInfo->HallSymbol;
for (i = 0; i <= MaxLenHallSymbol; i++)
*hsym++ = '\0';
PreviousRotation = 0;
PreviousRefAxis = 0;
PreviousDirCode = 0;
hsym = SgInfo->HallSymbol;
if (NeedDash)
*hsym++ = '-';
else
*hsym++ = ' ';
*hsym++ = SgInfo->LatticeInfo->Code;
nTrV = SgInfo->LatticeInfo->nTrVector;
for (iG = 0; iG < nGRT; iG++)
{
rmxi = GRT[iG].RMxI;
AbsOrder = abs(rmxi->Order);
RefAxis = rmxi->RefAxis;
DirCode = rmxi->DirCode;
if (RefAxis == 'o') RefAxis = 0;
if (DirCode == '=' || DirCode == '.') DirCode = 0;
if (iG == 0)
{
if (RefAxis == 'z') RefAxis = 0;
}
else
{
if (AbsOrder == 2)
{
if (PreviousRotation == 2 || PreviousRotation == 4)
{
if (RefAxis == 'x') RefAxis = 0;
}
else if (PreviousRotation == 3 || PreviousRotation == 6)
{
if ( PreviousDirCode == '*'
|| RefAxis == PreviousRefAxis) RefAxis = 0;
if (DirCode == '\'') DirCode = 0;
}
}
else if (AbsOrder == 3)
{
if (DirCode == '*') DirCode = 0;
}
}
PreviousRotation = AbsOrder;
PreviousRefAxis = rmxi->RefAxis;
PreviousDirCode = rmxi->DirCode;
*hsym++ = ' ';
if (rmxi->Order < 0) *hsym++ = '-';
*hsym++ = Digits[AbsOrder];
if (RefAxis) *hsym++ = RefAxis;
if (DirCode) *hsym++ = DirCode;
TrV = SgInfo->LatticeInfo->TrVector;
for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3)
{
for (i = 0; i < 3; i++)
if ((GRT[iG].Transl[i] + TrV[i]) % STBF != 0)
break;
if (i == 3)
break;
}
if (iTrV < nTrV)
continue; /* next iG */
hsym_mark = hsym;
TrV = SgInfo->LatticeInfo->TrVector;
for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3, hsym = hsym_mark)
{
for (i = 0; i < 3; i++)
Transl[i] = iModPositive(GRT[iG].Transl[i] + TrV[i], STBF);
Screw = 0;
for (ip = 0; ip < 3; ip++)
if (rmxi->EigenVector[ip] != 0) break;
if (ip < 3 && rmxi->EigenVector[ip] == 1)
{
for (i = ip + 1; i < 3; i++)
if (rmxi->EigenVector[i] != 0) break;
if (i == 3)
{
ScrewPrimitive = STBF / AbsOrder;
Screw = Transl[ip] / ScrewPrimitive;
i = Screw * ScrewPrimitive;
if (i % 3)
{
*hsym++ = Digits[Screw];
Transl[ip] -= i;
}
}
}
ht = HallTranslations;
while (*ht)
{
for (i = 0; i < 3; i++)
if (Transl[i] < ht[i + 1]) break;
if (i == 3)
{
for (i = 0; i < 3; i++)
Transl[i] -= ht[i + 1];
*hsym++ = (char) *ht;
}
ht += 4;
}
for (i = 0; i < 3; i++)
if (Transl[i] != 0)
break;
if (i == 3)
break;
}
if (iTrV == nTrV)
return 0;
}
if (nGRT == 0)
{
*hsym++ = ' ';
*hsym++ = '1';
}
if (HaveOrSh)
{
*hsym++ = ' ';
*hsym++ = '(';
for (i = 0; i < 3; i++)
{
if (i) *hsym++ = ' ';
os = iModPositive(SgInfo->OriginShift[i], 12);
if (os > 6)
{
*hsym++ = '-';
os = 12 - os;
}
*hsym++ = Digits[os];
}
*hsym++ = ')';
}
*hsym = '\0';
if (SgInfo->HallSymbol[MaxLenHallSymbol] != '\0') {
SetSgError("Internal Error: BuildHSym(): MaxLenHallSymbol too small");
return -1;
}
return 1;
}
static int BuildHallSymbol(T_SgInfo *SgInfo, int FixedOriginShift)
{
int ix, iy, iz;
int status;
static const int ShiftTable[] = { 0, 1, -1, 2, -2, 3 };
if (SgError != NULL) return -1;
if (SgInfo->nGenerator == 0)
{
if (BuildGenerator_iList(SgInfo) != 0)
{
SetSgError("Error: Can't build generator list");
return -1;
}
}
if (FixedOriginShift)
{
status = BuildHSym(SgInfo);
if (status == 1)
return 0;
}
else
{
for (ix = 0; ix < 6; ix++)
{
SgInfo->OriginShift[0] = ShiftTable[ix];
for (iy = 0; iy < 6; iy++)
{
SgInfo->OriginShift[1] = ShiftTable[iy];
for (iz = 0; iz < 6; iz++)
{
SgInfo->OriginShift[2] = ShiftTable[iz];
status = BuildHSym(SgInfo);
if (status < 0)
return -1;
if (status == 1)
return 0;
}
}
}
}
SetSgError("Error: Can't build Hall Symbol");
return -1;
}
void InitSgInfo(T_SgInfo *SgInfo)
{
int i;
SgInfo->GenOption = 0;
SgInfo->Centric = 0;
SgInfo->InversionOffOrigin = 0;
SgInfo->LatticeInfo = LI_P;
SgInfo->StatusLatticeTr = 0;
for (i = 0; i < 3; i++)
SgInfo->OriginShift[i] = 0;
SgInfo->nList = 0;
SgInfo->OrderL = 0;
SgInfo->OrderP = 0;
SgInfo->XtalSystem = XS_Unknown;
SgInfo->UniqueRefAxis = 0;
SgInfo->UniqueDirCode = 0;
SgInfo->ExtraInfo = EI_Unknown;
SgInfo->PointGroup = PG_Unknown;
SgInfo->nGenerator = 0;
SgInfo->HallSymbol[0] = '\0';
SgInfo->TabSgName = NULL;
SgInfo->CCMx_LP = NULL;
SgInfo->n_si_Vector = -1;
}
int CompleteSgInfo(T_SgInfo *SgInfo)
{
int List_iList[192];
const T_TabSgName *tsgn;
if (SgInfo->StatusLatticeTr == -1)
{
if (AddLatticeTr2ListSeitzMx(SgInfo, SgInfo->LatticeInfo) < 0)
return -1;
}
if (ApplyOriginShift(SgInfo) < 0)
return -1;
if (SgInfo->nList > sizeof List_iList / sizeof (*List_iList)) {
SetSgError("Internal Error: CompleteSgInfo()");
return -1;
}
if (SgInfo->nList > 1)
{
SortSgInfoList(SgInfo, List_iList);
if (SgError != NULL) return -1;
}
if (RemoveLatticeTr(SgInfo) != 0)
return -1;
if (RemoveInversion(SgInfo) != 0)
return -1;
TidyTranslation(SgInfo);
if (SgInfo->nList > 1)
{
SortSgInfoList(SgInfo, List_iList);
if (SgError != NULL) return -1;
}
SgInfo->OrderP = SgInfo->nList;
if (SgInfo->Centric == -1) SgInfo->OrderP *= 2;
SgInfo->OrderL = SgInfo->OrderP * SgInfo->LatticeInfo->nTrVector;
if (BuildHallSymbol(SgInfo, 0) != 0)
return -1;
for (tsgn = TabSgName; tsgn->HallSymbol; tsgn++)
if ( strcmp(tsgn->HallSymbol, SgInfo->HallSymbol) == 0
&& ( SgInfo->TabSgName == NULL
|| SgInfo->TabSgName == tsgn))
break;
if (SgInfo->TabSgName != NULL && tsgn->HallSymbol == NULL)
{
if (SgError) return -1;
printf("Internal Error: Input/Output HallSymbol mismatch: %s <> %s",
SgInfo->TabSgName->HallSymbol, SgInfo->HallSymbol);
SetSgError(SgErrorBuffer);
return -1;
}
if (tsgn->HallSymbol)
SgInfo->TabSgName = tsgn;
SgInfo->CCMx_LP = NULL;
switch (SgInfo->LatticeInfo->Code)
{
case 'P': SgInfo->CCMx_LP = CCMx_PP; break;
case 'A': SgInfo->CCMx_LP = CCMx_AP; break;
case 'B': SgInfo->CCMx_LP = CCMx_BP; break;
case 'C': SgInfo->CCMx_LP = CCMx_CP; break;
case 'I': SgInfo->CCMx_LP = CCMx_IP; break;
case 'R':
switch (SgInfo->UniqueRefAxis) {
case 0:
case 'z': SgInfo->CCMx_LP = CCMx_RP_z; break;
case 'y': SgInfo->CCMx_LP = CCMx_RP_y; break;
default: break;
}
break;
case 'S':
switch (SgInfo->UniqueRefAxis) {
case 0:
case 'y': SgInfo->CCMx_LP = CCMx_SP_y; break;
case 'x': SgInfo->CCMx_LP = CCMx_SP_x; break;
default: break;
}
break;
case 'T':
switch (SgInfo->UniqueRefAxis) {
case 0:
case 'x': SgInfo->CCMx_LP = CCMx_TP_x; break;
case 'z': SgInfo->CCMx_LP = CCMx_TP_z; break;
default: break;
}
break;
case 'F': SgInfo->CCMx_LP = CCMx_FP; break;
default:
break;
}
if (SgInfo->CCMx_LP == NULL) {
SetSgError("Internal Error: Illegal lattice code");
return -1;
}
return 0;
}
int CB_SMx(T_RTMx *CSiC,
const T_RTMx *CBMx, const T_RTMx *SMx, const T_RTMx *InvCBMx)
{
int i;
T_RTMx BufMx;
RTMxMultiply(&BufMx, SMx, InvCBMx, CTBF / STBF, CTBF);
RTMxMultiply(CSiC, CBMx, &BufMx, CRBF, CRBF * CTBF);
for (i = 0; i < 9; i++)
{
if (CSiC->s.R[i] % (CRBF * CRBF)) {
SetSgError("Internal Error: Corrupt CBMx/SMx/InvCBMx");
return -1;
}
CSiC->s.R[i] /= (CRBF * CRBF);
}
for (i = 0; i < 3; i++)
{
if (CSiC->s.T[i] % (CRBF * (CTBF / STBF))) {
SetSgError("Internal Error: Out of STBF range");
return -1;
}
CSiC->s.T[i] /= (CRBF * (CTBF / STBF));
}
return 0;
}
int TransformSgInfo(const T_SgInfo *SgInfo,
const T_RTMx *CBMx, const T_RTMx *InvCBMx,
T_SgInfo *BC_SgInfo)
{
int iList, f, i;
int nTrV, iTrV, nLoopInv, iLoopInv;
const int *TrV;
T_RTMx SMx, BC_SMx;
const T_RTMx *lsmx;
nLoopInv = Sg_nLoopInv(SgInfo);
nTrV = SgInfo->LatticeInfo->nTrVector;
TrV = SgInfo->LatticeInfo->TrVector;
for (iTrV = 0; iTrV < nTrV; iTrV++, TrV += 3)
{
for (iLoopInv = 0; iLoopInv < nLoopInv; iLoopInv++)
{
if (iLoopInv == 0) f = 1;
else f = -1;
lsmx = SgInfo->ListSeitzMx;
for (iList = 0; iList < SgInfo->nList; iList++, lsmx++)
{
for (i = 0; i < 9; i++)
SMx.s.R[i] = f * lsmx->s.R[i];
for (i = 0; i < 3; i++)
SMx.s.T[i] = f * lsmx->s.T[i] + TrV[i];
if (CB_SMx(&BC_SMx, CBMx, &SMx, InvCBMx) != 0)
return -1;
if (Add2ListSeitzMx(BC_SgInfo, &BC_SMx) < 0)
return -1;
}
}
}
return 0;
}
#undef SGCLIB_C__
syntax highlighted by Code2HTML, v. 0.9.1