src/IntRank2/CartMx2A.cc
/* This function will assume that the irreducible Cartesian spatial tensor
   array A (symmetric) has no asymmetry. That is, eta = 0 for A.  It will 
   return three Euler angles EA:{alpha,beta,gamma} that relate to the tensor
   PAS to the input array axes. In this case, i.e. eta=0, the Euler angle
   gamma is ill-defined and set to zero. The two other angles alpha and beta
   may be generated from the following:
                                                                           1
                                                   [ [   [   A         ] ] - ]
      1       [     2         ]                    | | 1 |    zz       | | 2 |
 A  = - del   | 3cos (beta)-1 |        beta = acos | | - | 2 ----- + 1 | |   |
  zz  2    zz [               ]                    | | 3 |   del       | |   |
                                                   [ [   [      zz     ] ]   ]

                                  [           A          ] 
                                  | 4          xz        |
                     alpha = acos | - ------------------ |
                                  | 3 del  * sin(2*beta) |
                                  [      zz              ]                   */
This function will assume that the irreducible Cartesian spatial tensor
array A (symmetric) has no asymmetry. That is, eta = 0 for A.  It will
return three Euler angles EA:{alpha,beta,gamma} that relate to the tensor
PAS to the input array axes.


This function will assume that the irreducible Cartesian spatial tensor array A ( symmetric ) has no asymmetry . That is eta = 0 for A . It will return three Euler angles EA : { alpha beta gamma } that relate to the tensor PAS to the input array axes .
void CartMx2A::SymCartEA()
  {
double TFCut = 1.e-14;
  string fct("SymCartEA");				
  ConvMeth = 2;						
  TrackConv(ConvMeth);					
  if(EtaVal)
    {
    C2Aerror(25, 1);					
    C2Aerror(26, fct, 1);				
    C2Afatal(30);					
    }
  if(A.is_zero())
    {
    C2Aerror(17, 1);					
    C2Aerror(26, fct, 1);				
    C2Afatal(30);					
    }
  if(!PASDelzz)
    {
    C2Aerror(18, 1);					
    C2Aerror(26, fct, 1);				
    C2Afatal(30);					
    }
  double Axx = A.getRe(0,0) - AisoVal;			
  double Axy = A.getRe(0,1) - AisoVal;			
  double Axz = A.getRe(0,2) - AisoVal;			
  double Ayy = A.getRe(1,1) - AisoVal;			
  double Ayz = A.getRe(1,2) - AisoVal;			
  double Azz = A.getRe(2,2) - AisoVal;			
  double beta;						
  double alpha;						
  beta  = acos(sqrt((1.0/3.0)*(2.0*Azz/PASDelzz+1.0)));	
  beta  = fabs(beta);					
  alpha = acos((4.0/3.0)*Axz/(PASDelzz*sin(2.0*beta)));	
  Sb    = sin(beta);					
  S2b   = sin(2.0*beta);				
  Ssqb  = Sb*Sb;					*sin(beta)
  Cb    = cos(beta);					
  Csqb  = Cb*Cb;					*cos(beta)
  Ca    = cos(alpha);					
  S2a   = sin(2.0*alpha); 				
  C2a   = cos(2.0*alpha); 				
  if((fabs(Axy - 0.75*PASDelzz*S2a*Ssqb)                 < TFCut)
  && (fabs(Axz - 0.75*PASDelzz*Ca*S2b)                   < TFCut)
  && (fabs(Axx - 0.25*PASDelzz*(3.*C2a*Ssqb-3.*Csqb+1.)) < TFCut)
  && (fabs(Ayy + 0.25*PASDelzz*(3.*C2a*Ssqb+3.*Csqb-1.)) < TFCut))
    { EA = EAngles(alpha, beta, 0.0); return; }
  if(beta < PI/2.0) beta = PI-beta;			
  else              beta = PI/2.0-beta;
  alpha = acos((4.0/3.0)*Axz/(PASDelzz*sin(2.0*beta)));	
  Sb    = sin(beta);					
  S2b   = sin(2.0*beta);				
  Ssqb  = Sb*Sb;					*sin(beta)
  Cb    = cos(beta);					
  Csqb  = Cb*Cb;					*cos(beta)
  Ca    = cos(alpha); 					
  S2a   = sin(2.0*alpha); 				
  C2a   = cos(2.0*alpha); 				
  if((fabs(Axy - 0.75*PASDelzz*S2a*Ssqb)                 < TFCut)
  && (fabs(Axz - 0.75*PASDelzz*Ca*S2b)                   < TFCut)
  && (fabs(Axx - 0.25*PASDelzz*(3.*C2a*Ssqb-3.*Csqb+1.)) < TFCut)
  && (fabs(Ayy + 0.25*PASDelzz*(3.*C2a*Ssqb+3.*Csqb-1.)) < TFCut))
    { EA = EAngles(alpha, beta, 0.0); return; }
  double beta2;
  if(beta < PI/2.0) beta2 = PI-beta;
  else              beta2 = PI/2.0-beta;
  beta2 = PI - beta + PI/2.0;				
  if(beta2 > PI) beta2 = beta - PI/2.0;		
  double b1alphaxz, b2alphaxz;
  b1alphaxz = acos((4.0/3.0)*Axz/(PASDelzz*sin(2.0*beta)));
  b2alphaxz = acos((4.0/3.0)*Axz/(PASDelzz*sin(2.0*beta2)));
  double b1alphayz, b2alphayz;
  b1alphayz = asin((4.0/3.0)*Ayz/(PASDelzz*sin(2.0*beta)));
  b2alphayz = asin((4.0/3.0)*Ayz/(PASDelzz*sin(2.0*beta2)));
  double diff1 = fabs(b1alphaxz - b1alphayz);
  double diff2 = fabs(b2alphaxz - b2alphayz);
  if(diff1 < diff2)
    { EA = EAngles(b1alphaxz, beta, 0.0); }
  else
    { EA = EAngles(b2alphaxz, beta2, 0.0); }
  TrackConv(16);					
  }
void CartMx2A :: SymCartEA ( ) { double TFCut = 1.e-14 ; string fct ( " " ) ; ConvMeth = 2 ; TrackConv ( ConvMeth ) ; if ( EtaVal ) { C2Aerror ( 25 , 1 ) ; C2Aerror ( 26 , fct , 1 ) ; C2Afatal ( 30 ) ; } if ( A . is_zero ( ) ) { C2Aerror ( 17 , 1 ) ; C2Aerror ( 26 , fct , 1 ) ; C2Afatal ( 30 ) ; } if ( ! PASDelzz ) { C2Aerror ( 18 , 1 ) ; C2Aerror ( 26 , fct , 1 ) ; C2Afatal ( 30 ) ; } double Axx = A . getRe ( 0 , 0 ) - AisoVal ; double Axy = A . getRe ( 0 , 1 ) - AisoVal ; double Axz = A . getRe ( 0 , 2 ) - AisoVal ; double Ayy = A . getRe ( 1 , 1 ) - AisoVal ; double Ayz = A . getRe ( 1 , 2 ) - AisoVal ; double Azz = A . getRe ( 2 , 2 ) - AisoVal ; double beta ; double alpha ; beta = acos ( sqrt ( ( 1.0 / 3.0 ) * ( 2.0 * Azz / PASDelzz + 1.0 ) ) ) ; beta = fabs ( beta ) ; alpha = acos ( ( 4.0 / 3.0 ) * Axz / ( PASDelzz * sin ( 2.0 * beta ) ) ) ; Sb = sin ( beta ) ; S2b = sin ( 2.0 * beta ) ; Ssqb = Sb * Sb ; Cb = cos ( beta ) ; Csqb = Cb * Cb ; Ca = cos ( alpha ) ; S2a = sin ( 2.0 * alpha ) ; C2a = cos ( 2.0 * alpha ) ; if ( ( fabs ( Axy - 0.75 * PASDelzz * S2a * Ssqb ) < TFCut ) && ( fabs ( Axz - 0.75 * PASDelzz * Ca * S2b ) < TFCut ) && ( fabs ( Axx - 0.25 * PASDelzz * ( 3. * C2a * Ssqb - 3. * Csqb + 1. ) ) < TFCut ) && ( fabs ( Ayy + 0.25 * PASDelzz * ( 3. * C2a * Ssqb + 3. * Csqb - 1. ) ) < TFCut ) ) { EA = EAngles ( alpha , beta , 0.0 ) ; return ; } if ( beta < PI / 2.0 ) beta = PI - beta ; else beta = PI / 2.0 - beta ; alpha = acos ( ( 4.0 / 3.0 ) * Axz / ( PASDelzz * sin ( 2.0 * beta ) ) ) ; Sb = sin ( beta ) ; S2b = sin ( 2.0 * beta ) ; Ssqb = Sb * Sb ; Cb = cos ( beta ) ; Csqb = Cb * Cb ; Ca = cos ( alpha ) ; S2a = sin ( 2.0 * alpha ) ; C2a = cos ( 2.0 * alpha ) ; if ( ( fabs ( Axy - 0.75 * PASDelzz * S2a * Ssqb ) < TFCut ) && ( fabs ( Axz - 0.75 * PASDelzz * Ca * S2b ) < TFCut ) && ( fabs ( Axx - 0.25 * PASDelzz * ( 3. * C2a * Ssqb - 3. * Csqb + 1. ) ) < TFCut ) && ( fabs ( Ayy + 0.25 * PASDelzz * ( 3. * C2a * Ssqb + 3. * Csqb - 1. ) ) < TFCut ) ) { EA = EAngles ( alpha , beta , 0.0 ) ; return ; } double beta2 ; if ( beta < PI / 2.0 ) beta2 = PI - beta ; else beta2 = PI / 2.0 - beta ; beta2 = PI - beta + PI / 2.0 ; if ( beta2 > PI ) beta2 = beta - PI / 2.0 ; double b1alphaxz , b2alphaxz ; b1alphaxz = acos ( ( 4.0 / 3.0 ) * Axz / ( PASDelzz * sin ( 2.0 * beta ) ) ) ; b2alphaxz = acos ( ( 4.0 / 3.0 ) * Axz / ( PASDelzz * sin ( 2.0 * beta2 ) ) ) ; double b1alphayz , b2alphayz ; b1alphayz = asin ( ( 4.0 / 3.0 ) * Ayz / ( PASDelzz * sin ( 2.0 * beta ) ) ) ; b2alphayz = asin ( ( 4.0 / 3.0 ) * Ayz / ( PASDelzz * sin ( 2.0 * beta2 ) ) ) ; double diff1 = fabs ( b1alphaxz - b1alphayz ) ; double diff2 = fabs ( b2alphaxz - b2alphayz ) ; if ( diff1 < diff2 ) { EA = EAngles ( b1alphaxz , beta , 0.0 ) ; } else { EA = EAngles ( b2alphaxz , beta2 , 0.0 ) ; } TrackConv ( 16 ) ; }
This function will assume that the irreducible Cartesian spatial tensor
array A (symmetric) has no asymmetry.
This function will assume that the irreducible Cartesian spatial tensor array A ( symmetric ) has no asymmetry .
// Our function name // Our conversion type // Track conv. process //              Insure We Are Dealing With Asymmetric, Non-Zero Tensor // Eta is not zero // Bad function call // Cant get EA // Cartesian array 0 // Bad function call // Cant get EA // Anisotropy is 0 // Bad function call // Cant get EA // Oriented Axx value // Oriented Axy value // Oriented Axz value // Oriented Ayy value // Oriented Ayz value // Oriented Azz value //              Determine  The Possible Values Of Beta From Azz //         Use Of cos(beta)*cos(beta) Makes The Result MultiValued //              Use Of acos Produces A Result Between [-90,90] // Here is angle beta // Here is angle alpha // Beta based on Azz // [-90,90] ==> [0,PI]  // Alpha based on Axz // sin(beta) // sin(2.0*beta) // sin(beta)*sin(beta) // cos(beta) // cos(beta)*cos(beta) // cos(alpha) // sin(2.0*alpha) // sin(2.0*alpha) // Try alternate beta // Alpha based on Axz // sin(beta) // sin(2.0*beta) // sin(beta)*sin(beta) // cos(beta) // cos(beta)*cos(beta) // cos(alpha) // sin(2.0*alpha) // sin(2.0*alpha) // sosiz - don't think (or know if) the rest of this is used anymore? // Here is other beta // Both are possible // Track conv. process