C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_MK_ORIENT    **
C**************************
C
CD-Routine: Calculate Orientation Matrix - LDM_MK_ORIENT
C This routine calculates the orientation matrix based on the current
C LDM parameters. This is normally used as a service routine to the
C subroutine LRL_GEN but may be called directly if required.
CD-end
CD-Fortran:
      SUBROUTINE LRL_MK_ORIENT (ORIENT)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Make the orientation matrix
C
C          (Used as a service routine for LRL_GEN but may be
C           called directly if desired)
C
C Author:  (Derived from make_orient in GENLAUE)
C
C Arguments:
C
      REAL ORIENT(3,3)
CD-Parameters:
C  ORIENT(3,3) r (W)    Returns the orientation matrix
CD-end
C Notes: Generalised matrix builder
C
C    The reciprocal cell parameters are calculated from the real cell
C    parameters.
C
C    ORIENT returns the orientation matrix
C      ... row 1 is components of reciprocal cell along X,
C      ... row 2 along Y, row 3 along Z.
C    Convention for crystal system:
C
C         1 - triclinic    2 - monoclinic    3 - orthorhombic
C         4 - tetragonal   5 - hexagonal     6 - rhombohedral
C         7 - cubic
C
C    Conventions for alignment are:
C       1..astar,   2..bstar,   3..cstar
C      -1..-astar, -2..-bstar, -3..-cstar
C
C    Thus if the astar axis lies along the X-axis in the direction of +x
C    then IBEAM = 1, and similarly for IROT.
C
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      CHARACTER*1 VALSTR
      INTEGER I, J, K
      INTEGER IFLAG, ND
      INTEGER SI, SJ, SK
      INTEGER ISYST, IBEAM, IROT
      REAL    CELLP(6)
      REAL    RCELL(6)
      REAL    RL(3)
      REAL    CC(3), SS(3)
      REAL    COS_PSI, COS_RHO, A, C, S, FUN
      REAL    DEG_TO_RAD
      REAL    ALPHA, BETA, GAMMA
      REAL    WSTD
C
C-------------------------------------------------------------------------------
C
C====== Initialise & get LDM parameters
C
      WSTD = 1.0
      CALL LDM_GET_A (CELLP(1), IFLAG, ND)
      CALL LDM_GET_B (CELLP(2), IFLAG, ND)
      CALL LDM_GET_C (CELLP(3), IFLAG, ND)
      CALL LDM_GET_ALPHA (CELLP(4), IFLAG, ND)
      CALL LDM_GET_BETA (CELLP(5), IFLAG, ND)
      CALL LDM_GET_GAMMA (CELLP(6), IFLAG, ND)
      CALL LDM_GET_SYST (VALSTR, ISYST, IFLAG)
      CALL LDM_GET_BEAM (VALSTR, IBEAM, IFLAG)
      CALL LDM_GET_ROT (VALSTR, IROT, IFLAG)
C
C====== Get reciprocal cell parameters
C
      CALL LRL_RECCEL (RCELL,CELLP,WSTD)
C
C====== Clear the orientation matrix
C
      DO 80 I = 1 , 3
         DO 70 J = 1 , 3
            ORIENT(I,J) = 0.
 70      CONTINUE
 80   CONTINUE
C
C====== i, j, k have special meanings in this code.
C====== i refers to the crystal rotation axis, Z
C
      I = IABS(IROT)
      SI = 1
      IF (IROT.LT.0) SI = -1
C
C====== j refers to IBEAM, the x-ray beam axis, X
C
      J = IABS(IBEAM)
      SJ = 1
      IF (IBEAM.LT.0) SJ = -1
C
C====== k is therefore the other axis, Y
C
      K = 6 - I - J
C
C====== The sign of k relies on i,j,k and X,Y,Z being right-handed
C
      IF (K.EQ.J+1 .OR. K.EQ.J-2) THEN
         SK = SI*SJ
      ELSE
         SK = -SI*SJ
      END IF
C
C====== Zero orientation matrix
C
      CALL MTX_ZERO(ORIENT,3,3)
C
C====== Convert cell angles to radians
C
      DEG_TO_RAD = 3.1415926536/180.0
      ALPHA = RCELL(4)*DEG_TO_RAD
      BETA = RCELL(5)*DEG_TO_RAD
      GAMMA = RCELL(6)*DEG_TO_RAD
C
C====== Set reciprocal cell lengths
C
      RL(1) = RCELL(1)
      RL(2) = RCELL(2)
      RL(3) = RCELL(3)
C
C====== Now go and do specific stuff for each crystal system
C
      IF (ISYST.EQ.2) THEN
C
C====== Monoclinic
C
         CC(1) = 0.0
         CC(2) = COS(BETA)
         CC(3) = 0.0
C
         SS(1) = 1.0
         SS(2) = SIN(BETA)
         SS(3) = 1.0
C
         ORIENT(1,J) = SJ*RL(J)*SS(K)
         ORIENT(1,K) = SJ*RL(K)*CC(I)
         ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
         ORIENT(3,I) = SI*RL(I)
         ORIENT(3,J) = SI*RL(J)*CC(K)
         ORIENT(3,K) = SI*RL(K)*CC(J)

      ELSE IF (ISYST.EQ.3) THEN
C
C====== Orthorhombic
C
         ORIENT(1,J) = SJ*RL(J)
         ORIENT(2,K) = SK*RL(K)
         ORIENT(3,I) = SI*RL(I)

      ELSE IF (ISYST.EQ.4) THEN
C
C====== Tetragonal
C
         ORIENT(1,J) = SJ*RL(J)
         ORIENT(2,K) = SK*RL(K)
         ORIENT(3,I) = SI*RL(I)

      ELSE IF (ISYST.EQ.5) THEN
C
C====== Hexagonal
C
         CC(1) = 0.0
         CC(2) = 0.0
         CC(3) = 0.5
C
         SS(1) = 1.0
         SS(2) = 1.0
         SS(3) = SQRT(3.0)/2.0
C
         ORIENT(1,J) = SJ*RL(J)*SS(K)
         ORIENT(1,K) = SJ*RL(K)*CC(I)
         ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
         ORIENT(3,I) = SI*RL(I)
         ORIENT(3,J) = SI*RL(J)*CC(K)
         ORIENT(3,K) = SI*RL(K)*CC(J)

      ELSE IF (ISYST.EQ.6) THEN
C
C====== Rhombohedral
C
         A = RCELL(1)
         C = COS(ALPHA)
         S = SIN(ALPHA)
         FUN = SQRT(1.-3.*C**2+2.*C**3)
C
         ORIENT(1,J) = SJ*A*S
         ORIENT(1,K) = SJ*A*C*(1.-C)/S
         ORIENT(2,K) = SK*A*FUN/S
         ORIENT(3,I) = SI*A
         ORIENT(3,J) = SI*A*C
         ORIENT(3,K) = SI*A*C

      ELSE IF (ISYST.EQ.7) THEN
C
C====== Cubic
C
         A = RCELL(1)
C
         ORIENT(1,J) = SJ*A
         ORIENT(2,K) = SK*A
         ORIENT(3,I) = SI*A

      ELSE
C
C====== Triclinic
C
         CC(1) = COS(ALPHA)
         CC(2) = COS(BETA)
         CC(3) = COS(GAMMA)
C
         SS(1) = SIN(ALPHA)
         SS(2) = SIN(BETA)
         SS(3) = SIN(GAMMA)
C
C====== See Buerger's Xray Crystallography for the next formulae, p348
C
         COS_PSI = (CC(I)-CC(K)*CC(J))/SS(K)
         COS_RHO = SQRT(1.0-CC(I)**2-CC(J)**2-CC(K)**2+2.0*CC(I)*CC(J)
     +             *CC(K))/SS(K)
C
         ORIENT(1,J) = SJ*RL(J)*SS(K)
         ORIENT(1,K) = SJ*RL(K)*COS_PSI
         ORIENT(2,K) = SK*RL(K)*COS_RHO
         ORIENT(3,I) = SI*RL(I)
         ORIENT(3,J) = SI*RL(J)*CC(K)
         ORIENT(3,K) = SI*RL(K)*CC(J)
      END IF
      RETURN
      END
