45 theSpaceGroup(spacegroup),
59 cosa=std::cos(alpha), cosb=std::cos(beta), cosg=std::cos(gamma);
60 sina=std::sin(alpha), sinb=std::sin(beta), sing=std::sin(gamma);
62 cosar = (cosb*cosg-cosa)/(sinb*sing);
63 cosbr = (cosa*cosg-cosb)/(sina*sing);
64 cosgr = (cosa*cosb-cosg)/(sina*sinb);
67 theRecVolume = 1. / theVolume;
69 theRecSize[0] = sizeB * sizeC * sina / theVolume;
70 theRecSize[1] = sizeC * sizeA * sinb / theVolume;
71 theRecSize[2] = sizeA * sizeB * sing / theVolume;
92 x3=cosa, y3=(cosb-cosa*cosg)/sing, z3=std::sqrt(1.-x3*x3-y3*y3);
102 x3=cosa, y3=(cosb-cosa*cosg)/sing, z3=std::sqrt(1.-x3*x3-y3*y3);
151 if( aGroup >= 1 && aGroup <= 2 ) {
return Triclinic;}
152 else if(aGroup >= 3 && aGroup <= 15 ) {
return Monoclinic;}
153 else if(aGroup >= 16 && aGroup <= 74 ) {
return Orthorhombic;}
154 else if(aGroup >= 75 && aGroup <= 142) {
return Tetragonal;}
155 else if(aGroup == 146 || aGroup == 148 ||
156 aGroup == 155 || aGroup == 160 ||
157 aGroup == 161 || aGroup == 166 ||
159 else if(aGroup >= 143 && aGroup <= 167) {
return Hexagonal;}
160 else if(aGroup >= 168 && aGroup <= 194) {
return Hexagonal;}
161 else if(aGroup >= 195 && aGroup <= 230) {
return Cubic;}
201 G4double x3=cosa, y3=(cosb-cosa*cosg)/sing, z3=std::sqrt(1.-x3*x3-y3*y3);
210 vecout.push_back(aaa);
219 for(
auto &vec:vecout){
232 return FillAmorphous(Cij);
235 return FillCubic(Cij);
238 return FillTetragonal(Cij);
241 return FillOrthorhombic(Cij);
244 return FillRhombohedral(Cij);
247 return FillMonoclinic(Cij);
250 return FillTriclinic(Cij);
253 return FillHexagonal(Cij);
264 G4bool G4CrystalUnitCell::FillAmorphous(
G4double Cij[6][6])
const {
265 Cij[3][3] = 0.5*(Cij[0][0]-Cij[0][1]);
272 G4double C11=Cij[0][0], C12=Cij[0][1], C44=Cij[3][3];
274 for (
size_t i=0; i<6; i++) {
275 for (
size_t j=i; j<6; j++) {
276 if (i<3 && j<3) Cij[i][j] = (i==j) ? C11 : C12;
277 else if (i==j && i>=3) Cij[i][i] = C44;
282 ReflectElReduced(Cij);
284 return (C11!=0. && C12!=0. && C44!=0.);
289 G4bool G4CrystalUnitCell::FillTetragonal(
G4double Cij[6][6])
const {
290 G4double C11=Cij[0][0], C12=Cij[0][1], C13=Cij[0][2], C16=Cij[0][5];
291 G4double C33=Cij[2][2], C44=Cij[3][3], C66=Cij[5][5];
298 ReflectElReduced(Cij);
301 return (C11!=0. && C12!=0. && C13!=0. && C33!=0. && C44!=0. && C66!=0.);
306 G4bool G4CrystalUnitCell::FillOrthorhombic(
G4double Cij[6][6])
const {
308 ReflectElReduced(Cij);
311 for (
size_t i=0; i<6; i++) {
312 for (
size_t j=i+1; j<3; j++)
313 good &= (Cij[i][j] != 0);
321 G4bool G4CrystalUnitCell::FillRhombohedral(
G4double Cij[6][6])
const {
322 G4double C11=Cij[0][0], C12=Cij[0][1], C13=Cij[0][2], C14=Cij[0][3];
323 G4double C15=Cij[0][4], C33=Cij[2][2], C44=Cij[3][3], C66=0.5*(C11-C12);
334 return (C11!=0 && C12!=0 && C13!=0 && C14!=0. &&
335 C33!=0. && C44!=0. && C66!=0.);
340 G4bool G4CrystalUnitCell::FillMonoclinic(
G4double Cij[6][6])
const {
344 return (FillOrthorhombic(Cij) && Cij[0][5]!=0. && Cij[1][5]!=0. &&
345 Cij[2][5] != 0. && Cij[3][4]!=0.);
350 G4bool G4CrystalUnitCell::FillTriclinic(
G4double Cij[6][6])
const {
353 ReflectElReduced(Cij);
356 for (
size_t i=0; i<6; i++) {
357 for (
size_t j=i; j<6; j++) good &= (Cij[i][j] != 0);
366 G4bool G4CrystalUnitCell::FillHexagonal(
G4double Cij[6][6])
const {
368 Cij[4][5] = 0.5*(Cij[0][0] - Cij[0][1]);
375 G4bool G4CrystalUnitCell::ReflectElReduced(
G4double Cij[6][6])
const {
376 for (
size_t i=1; i<6; i++) {
377 for (
size_t j=i+1; j<6; j++) {
378 Cij[j][i] = Cij[i][j];
404 return a*a*a*std::sqrt(1.-3.*cosa*cosa+2.*cosa*cosa*cosa);
410 return a*b*c*std::sqrt(1.-cosa*cosa-cosb*cosb-cosg*cosg*2.*cosa*cosb*cosg);
413 return std::sqrt(3.0)/2.*a*a*c;
446 G4double a2 = a*a, b2 = b*b, c2 = c*c;
447 G4double h2 = h*h, k2 = k*k, l2 = l*l;
458 return a2 / ( h2+k2+l2 );
461 return 1.0 / ( (h2 + k2)/a2 + l2/c2 );
464 return 1.0 / ( h2/a2 + k2/b2 + l2/c2 );
467 cos2a=cosa*cosa; sin2a=sina*sina;
468 T = h2+k2+l2+2.*(h*k+k*l+h*l) * ((cos2a-cosa)/sin2a);
469 R = sin2a / (1. - 3*cos2a + 2.*cos2a*cosa);
474 return 1./(1./sin2b * (h2/a2+l2/c2-2*h*l*cosb/(a*c)) + k2/b2);
480 return 1. / ( (4.*(h2+k2+h*k) / (3.*a2)) + l2/c2 );
512 G4double a2 = a*a, b2 = b*b, c2 = c*c;
513 G4double h2 = h*h, k2 = k*k, l2 = l*l;
521 return a2 * (h2+k2+l2);
524 return (h2+k2)*a2 + l2*c2 ;
527 return h2*a2 + k2+b2 + h2*c2;
530 return (h2+k2+l2+2.*(h*k+k*l+h*l) * cosar)*a2;
533 return h2*a2+k2*b2+l2*c2+2.*h*l*a*c*cosbr;
536 return h2*a2+k2*b2+l2*c2+2.*k*l*b*c*cosar+2.*l*h*c*a*cosbr+2.*h*k*a*b*cosgr;
539 return (h2+k2+h*k)*a2 + l2*c2;
573 G4double a2 = a*a, b2 = b*b, c2 = c*c;
581 return (h1*h2 + k1*k2 + l1+l2) / (std::sqrt(h1*h1 + k1*k1 + l1*l1) * std::sqrt(h2*h2 + k2*k2 + l2*l2));
589 return dsp1dsp2 * (h1*h2*a2 + k1*k2*a2 + l1*l2*c2);
593 return dsp1dsp2 * (h1*h2*a2 + k1*k2*b2 + l1*l2*c2+
594 (k1*l2+k2*l1)*b*c*cosar+
595 (h1*l2+h2*l1)*a*c*cosbr+
596 (h1*k2+h2*k1)*a*b*cosgr);
600 return dsp1dsp2 * (h1*h2*a2 + k1*k2*b2 + l1*l2*c2+
601 (k1*l2+k2*l1)*b*c*cosar+
602 (h1*l2+h2*l1)*a*c*cosbr+
603 (h1*k2+h2*k1)*a*b*cosgr);
607 return dsp1dsp2 * (h1*h2*a2 + k1*k2*b2 + l1*l2*c2+
608 (k1*l2+k2*l1)*b*c*cosar+
609 (h1*l2+h2*l1)*a*c*cosbr+
610 (h1*k2+h2*k1)*a*b*cosgr);
614 return dsp1dsp2 *( (h1*h2 + k1*k2 + 0.5*(h1*k2+k1*h2))*a2 + l1*l2*c2);
const G4ThreeVector & GetRecBasis(G4int idx) const
DLL_API const Hep3Vector HepZHat
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector theRecAngle
G4double GetIntCosAng(G4int h1, G4int k1, G4int l1, G4int h2, G4int k2, G4int l2)
const G4ThreeVector & GetBasis(G4int idx) const
G4double GetIntSp2(G4int h, G4int k, G4int l)
G4ThreeVector theRecBasis[3]
virtual ~G4CrystalUnitCell()
G4double ComputeCellVolume()
Hep3Vector & rotateZ(double)
G4ThreeVector theUnitBasis[3]
theLatticeSystemType GetLatticeSystem()
const G4ThreeVector & GetUnitBasis(G4int idx) const
DLL_API const Hep3Vector HepYHat
G4bool FillElReduced(G4double Cij[6][6])
G4bool FillAtomicPos(G4ThreeVector &pos, std::vector< G4ThreeVector > &vecout)
static constexpr double deg
G4double GetRecIntSp2(G4int h, G4int k, G4int l)
DLL_API const Hep3Vector HepXHat
G4CrystalUnitCell(G4double sizeA, G4double sizeB, G4double sizeC, G4double alpha, G4double beta, G4double gamma, G4int spacegroup)
G4ThreeVector GetUnitBasisTrigonal()
G4ThreeVector theBasis[3]
static constexpr double halfpi
G4bool FillAtomicUnitPos(G4ThreeVector &pos, std::vector< G4ThreeVector > &vecout)
Hep3Vector & rotateX(double)
G4ThreeVector theRecUnitBasis[3]
static const G4double alpha
static const G4double pos
const G4ThreeVector & GetRecUnitBasis(G4int idx) const