C----- The rotation matrix is stored in the array ROTATE.
C----- Use single precision on cray C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION PROP(16), ROTATE(3,3), TERM1(3,3), TERM2(3,3), INDX(3)
C----- Subroutines: C
C CROSS -- cross product of two vectors C
C LUDCMP -- LU decomposition C
C LUBKSB -- linear equation solver based on LU decomposition C method (must call LUDCMP first)
C----- PROP -- constants characterizing the crystal orientation C (INPUT) C
C PROP(1) - PROP(3) -- direction of the first vector in C local cubic crystal system C PROP(4) - PROP(6) -- direction of the first vector in C global system C
C PROP(9) - PROP(11)-- direction of the second vector in C local cubic crystal system
C PROP(12)- PROP(14)-- direction of the second vector in C global system C
C----- ROTATE -- rotation matrix (OUTPUT): C
C ROTATE(i,1) -- direction cosines of direction [1 0 0] in C local cubic crystal system
C ROTATE(i,2) -- direction cosines of direction [0 1 0] in C local cubic crystal system
C ROTATE(i,3) -- direction cosines of direction [0 0 1] in C local cubic crystal system
C----- local matrix: TERM1
CALL CROSS (PROP(1), PROP(9), TERM1, ANGLE1)
C----- LU decomposition of TERM1
CALL LUDCMP (TERM1, 3, 3, INDX, DCMP)
C----- inverse matrix of TERM1: TERM2 DO J=1,3 DO I=1,3
IF (I.EQ.J) THEN TERM2(I,J)=1. ELSE
TERM2(I,J)=0. END IF END DO END DO
DO J=1,3
CALL LUBKSB (TERM1, 3, 3, INDX, TERM2(1,J)) END DO
C----- global matrix: TERM1
CALL CROSS (PROP(4), PROP(12), TERM1, ANGLE2)
C----- Check: the angle between first and second vector in local and C global systems must be the same. The relative difference must be C less than 0.1%. C
IF (ABS(ANGLE1/ANGLE2-1.).GT.0.001) THEN WRITE (6,*)
2 '***ERROR - angles between two vectors are not the same' STOP END IF
C----- rotation matrix: ROTATE DO J=1,3 DO I=1,3
ROTATE(I,J)=0. DO K=1,3
ROTATE(I,J)=ROTATE(I,J)+TERM1(I,K)*TERM2(K,J) END DO END DO END DO
RETURN END
C-----------------------------------
SUBROUTINE CROSS (A, B, C, ANGLE)
C----- (1) normalize vectors A and B to unit vectors C (2) store A, B and A*B (cross product) in C
C----- Use single precision on cray C
IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(3), B(3), C(3,3)
SUM1=SQRT(A(1)**2+A(2)**2+A(3)**2) SUM2=SQRT(B(1)**2+B(2)**2+B(3)**2)
IF (SUM1.EQ.0.) THEN
WRITE (6,*) '***ERROR - first vector is zero' STOP ELSE
DO I=1,3
C(I,1)=A(I)/SUM1 END DO END IF
IF (SUM2.EQ.0.) THEN
WRITE (6,*) '***ERROR - second vector is zero' STOP ELSE
DO I=1,3
C(I,2)=B(I)/SUM2 END DO END IF
ANGLE=0. DO I=1,3
ANGLE=ANGLE+C(I,1)*C(I,2) END DO
ANGLE=ACOS(ANGLE)
C(1,3)=C(2,1)*C(3,2)-C(3,1)*C(2,2) C(2,3)=C(3,1)*C(1,2)-C(1,1)*C(3,2) C(3,3)=C(1,1)*C(2,2)-C(2,1)*C(1,2)
SUM3=SQRT(C(1,3)**2+C(2,3)**2+C(3,3)**2)
IF (SUM3.LT.1.E-8) THEN WRITE (6,*)
2 '***ERROR - first and second vectors are parallel' STOP END IF
RETURN END
C----------------------------------------------------------------------
SUBROUTINE SLIPSYS (ISPDIR, ISPNOR, NSLIP, SLPDIR, SLPNOR, 2 ROTATE)
C----- This subroutine generates all slip systems in the same set for C a CUBIC crystal. For other crystals (e.g., HCP, Tetragonal, C Orthotropic, ...), it has to be modified to include the effect of C crystal aspect ratio.
C----- Denote s as a slip direction and m as normal to a slip plane. C In a cubic crystal, (s,-m), (-s,m) and (-s,-m) are NOT considered C independent of (s,m).
C----- Subroutines: LINE1 and LINE
C----- Variables: C
C ISPDIR -- a typical slip direction in this set of slip systems C (integer) (INPUT)
C ISPNOR -- a typical normal to slip plane in this set of slip C systems (integer) (INPUT)
C NSLIP -- number of independent slip systems in this set C (OUTPUT)
C SLPDIR -- unit vectors of all slip directions (OUTPUT) C SLPNOR -- unit normals to all slip planes (OUTPUT) C ROTATE -- rotation matrix (INPUT)
C ROTATE(i,1) -- direction cosines of [100] in global system C ROTATE(i,2) -- direction cosines of [010] in global system C ROTATE(i,3) -- direction cosines of [001] in global system C
C NSPDIR -- number of all possible slip directions in this set C NSPNOR -- number of all possible slip planes in this set