C  **********************************************
C SUBROUTINE TSTRJM ====*====3====*====4====*====5====*====6====*====7
C
C    SPHERICAL HARMONICS(J/2) ARE TRANSFORMED TO CUBIC(IL.GE.1)
C    HARMINICS OF DOUBLE REPRESENTATION
C           2000/01/11  Akira Yanase
C
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
      SUBROUTINE TSTRJM(JJ,U,KP,INS)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SN,J3HM
      COMMON/SPG2/IL,NGNG,IIGG(48),JV(2,3,48)
      COMMON/SPG6/SN(2,2,24),IDF(24,24)
      COMMON/JEQ3/J3HM(4,4,24)
      SAVE /SPG2/,/SPG6/,/JEQ3/,TC,TH,NDIMC,NDIMH
      COMPLEX*16 U(100),WD(24),WE(4,24),CW,WW
      COMPLEX*16 JJHM(22,22,24)
      DIMENSION  EULER(3,24)
      INTEGER ICREAL(2,3,24)
      INTEGER ICM(16,22,22),K1M(22,22),K2M(22,22),K3M(22,22),ISM(22,22)
      DIMENSION IC(16),KP(100),INS(4,24)
      INTEGER TC(24),TH(12),NDIMC(3),NDIMH(3)
      DATA TC/
     &  1,2,2,2,3,3,3,3,3,3,3,3,
     &  2,2,2,2,2,2,4,4,4,4,4,4/
      DATA TH/
     &  1,6,3,2,3,6,2,2,2,2,2,2/
      DATA NDIMC/
     &    2,2,4/
      DATA NDIMH/
     &    2,2,2/
      PAI=4.D0*DATAN(1.D0)
      CALL SPINM4
C
      IF(IL.GE.1) THEN
         CALL EULERC(ICREAL,EULER)
         JRCH=1
         NG=24
      ELSE
         WRITE(6,*) ' This TSTRJM is only for cubic group'
         STOP
      END IF
      IPR=0
      JDIM=JJ+1
      DO 31 J=1,JDIM
      DO 31 K=1,JDIM
         MA=JJ+2-2*J
         MB=JJ+2-2*K
         CALL TSRMHI(JJ,MA,MB,K1,K2,K3,IS,IC,IPR)
         DO 32 II=1,16
           ICM(II,J,K)=IC(II)
   32    CONTINUE
         K1M(J,K)=K1
         K2M(J,K)=K2
         K3M(J,K)=K3
         ISM(J,K)=IS
   31 CONTINUE
      DO 10 I=1,NG
C       WRITE(6,660) (EULER(K,I),K=1,3)
C 660   FORMAT(1H ,3F10.5)
        CB=COS(EULER(2,I))
        W1=COS(EULER(2,I)*0.5D0)
        W3=SIN(EULER(2,I)*0.5D0)
        DO 20 J=1,JDIM
        DO 21 K=1,JDIM
          MA=JJ+2-2*J
          MB=JJ+2-2*K
          WW=0.0
          DO 22 II=1,16
            IF(ICM(II,J,K).NE.0) THEN
            WW=WW+ICM(II,J,K)*(CB**(II-1))
            END IF
   22     CONTINUE
          WW=((WW*K1M(J,K))/K2M(J,K))*SQRT(DBLE(K3M(J,K)))
          IF(ISM(J,K).EQ.1) THEN
             WW=WW*W1
          ELSE IF(ISM(J,K).EQ.3) THEN
             WW=WW*W3
          END IF
          CA=COS(EULER(1,I)*0.5D0*MA)
          CC=COS(EULER(3,I)*0.5D0*MB)
          SA=SIN(EULER(1,I)*0.5D0*MA)
          SC=SIN(EULER(3,I)*0.5D0*MB)
          JJHM(J,K,I)=DCMPLX(CA,-SA)*DCMPLX(CC,-SC)*WW
   21   CONTINUE
   20   CONTINUE
   10 CONTINUE
      ND1=0
      IW=0
      NJR=3
      DO 801 JR=1,NJR
      NDIM=NDIMC(JR)
C
C   CHARCTER TEST
C
      WB=0.0D0
      DO 204 IG=1,NG
      CL=JJ+1
      IF(IG.EQ.1) GO TO 201
      IF(JRCH.EQ.1) XX=PAI/TC(IG)
      CL=SIN((JJ+1)*XX)/SIN(XX)
  201 CONTINUE
C     WRITE(6,693) L,IG,CL
C 693 FORMAT(2I5,F8.4)
      IF(ABS(CL).LT.0.001D0) GO TO 204
      WA=0.0D0
      DO 205 IN=1,NDIM
      IF(JRCH.EQ.1) THEN
        IF(JR.EQ.1.OR.(JR.EQ.2.AND.IG.LE.12)) THEN
          WA=WA+DCONJG(SN(IN,IN,IG))
        ELSE IF(JR.EQ.2.AND.IG.GE.13) THEN
          WA=WA-DCONJG(SN(IN,IN,IG))
        ELSE IF(JR.EQ.3) THEN
          WA=WA+DCONJG(J3HM(IN,IN,IG))
        ELSE
          WRITE(6,*) ' SOMETING WRONG IN TSTRJM'
          STOP
        END IF
      ELSE 
          WRITE(6,*) ' SOMETING WRONG IN TSTRJM'
          STOP
      END IF
  205 CONTINUE
C     WRITE(6,691) WA,CL
C 691 FORMAT(8F10.4)
      IF(ABS(WA).LT.0.001D0) GO TO 204
      WB=WB+CL*WA
  204 CONTINUE
      WB=WB/NG
      NI=WB+0.5D0
C     WRITE(6,690) L,JR,NI
C 690 FORMAT(16I5)
      IF(NI.EQ.0) GO TO 801
C
C  PROJECTION OPERATOR
C
      NC=1
      NB=1
      IB=1
  150 CONTINUE
      MB=JJ-2*(IB-1)
C      WRITE(6,695) JJ,JR,NI,MB,NC,NB
C  695 FORMAT(' J=',I4,'  JR=',I4,'  NI=',6I4)
      DO 101 IA=1,JJ+1
      WD(IA)=0.0D0
      DO 104 IG=1,NG
      IF(ABS(JJHM(IA,IB,IG)).LE.0.000001D0) GO TO 104
      IF(JRCH.EQ.1) THEN
         IF(JR.EQ.1.OR.(JR.EQ.2.AND.IG.LE.12)) THEN
           WD(IA)=WD(IA)+JJHM(IA,IB,IG)*DCONJG(SN(NC,NB,IG))
         ELSE IF(JR.EQ.2.AND.IG.GE.13) THEN
           WD(IA)=WD(IA)-JJHM(IA,IB,IG)*DCONJG(SN(NC,NB,IG))
         ELSE IF(JR.EQ.3) THEN
           WD(IA)=WD(IA)+JJHM(IA,IB,IG)*DCONJG(J3HM(NC,NB,IG))
         END IF
      END IF
  104 CONTINUE
  101 CONTINUE
C
C    ORTHOGONALITY TEST
C
      IF(ND1.EQ.0) GO TO 111
      DO 112 ID=1,ND1
      ID1=INS(1,ID)
      ID2=INS(2,ID)
      WW=0.0D0
      DO 113 I=ID1,ID2
      JU=KP(I)
      WW=WW+CONJG(U(I))*WD(JU)
  113 CONTINUE
      IF(ABS(WW).LE.1.0D-4) GO TO 112
C      write(6,*) ' NON ORTHOGONALITY TO',ID,WW
C
C    ORTHOGONALIZATION
C
      DO 121 I=INS(1,ID),INS(2,ID)
      JU=KP(I)
      WD(JU)=WD(JU)-U(I)*WW
  121 CONTINUE
C
C
  112 CONTINUE
C
C   REGISTRATION
C
  111 ND1=ND1+1
      INS(1,ND1)=IW+1
      WA=0.0D0
      IIWW=IW
      DO 114 IA=1,J+1
      IF(ABS(WD(IA)).LT.1.0D-4) GO TO 114
      IW=IW+1
      KP(IW)=IA
      U(IW)=WD(IA)
      WA=WA+DBLE(U(IW)*DCONJG(U(IW)))
  114 CONTINUE
      WA=SQRT(WA)
      IF(WA.GT.1.0D-4) GO TO 116
      IW=IIWW
      ND1=ND1-1
      GO TO 303
  116 INS(2,ND1)=IW
      INS(3,ND1)=NC
      INS(4,ND1)=JR
      CW=1.0/WA
      DO 117 I=INS(1,ND1),INS(2,ND1)
  117 U(I)=U(I)*CW
C
C    PARTNER
C
      DO 131 IA=1,JJ+1
      WE(2,IA)=0
      WE(3,IA)=0
      WE(4,IA)=0
      DO 137 I=INS(1,ND1),INS(2,ND1)
      DO 134 IG=1,NG
      DO 136 N=2,NDIM
        IF(JRCH.EQ.1) THEN
         IF(JR.EQ.1.OR.(JR.EQ.2.AND.IG.LE.12)) THEN
           WE(N,IA)=WE(N,IA)+JJHM(IA,KP(I),IG)*DCONJG(SN(N,1,IG))*U(I)
         ELSE IF(JR.EQ.2.AND.IG.GT.12) THEN
           WE(N,IA)=WE(N,IA)-JJHM(IA,KP(I),IG)*DCONJG(SN(N,1,IG))*U(I)
         ELSE IF(JR.EQ.3) THEN
           WE(N,IA)=WE(N,IA)+JJHM(IA,KP(I),IG)*DCONJG(J3HM(N,1,IG))*U(I)
         END IF
        END IF
  136 CONTINUE   
  134 CONTINUE
  137 CONTINUE
  131 CONTINUE
C
C     ORTHOGONARITY TEST FOR PARTNER
C
      DO 141 ID=1,ND1
         ID1=INS(1,ID)
         ID2=INS(2,ID)
         DO 142 N=2,NDIM
            WW=0.0D0
            DO 143 I=ID1,ID2
            JU=KP(I)
            WW=WW+CONJG(U(I))*WE(N,JU)
 143        CONTINUE
            IF(ABS(WW).GE.1.0D-10) THEN
            WRITE(6,666) JR,N,ID,WW
 666           FORMAT(' JR=',I3,'  N=',I3
     & ,' THIS PARTNER IS NONORTHOGONAL TO',I5,2D15.3)
            END IF
 142     CONTINUE
 141  CONTINUE   
C
C    REGISTRATION OF PARTNER
C
      DO 138 N=2,NDIM
      ND1=ND1+1
      INS(1,ND1)=IW+1
      WA=0.0D0
      DO 139 IA=1,JJ+1
         IF(ABS(WE(N,IA)).LT.1.D-4) GO TO 139
         IW=IW+1
         KP(IW)=IA
         U(IW)=WE(N,IA)
         WA=WA+DBLE(U(IW)*DCONJG(U(IW)))
 139  CONTINUE
      IF(WA.LT.1.0D-4) THEN
         WRITE(6,699) JJ,JR,N
 699     FORMAT(' PARTNER BECOMES ZERO VECTOR',3I5)
         STOP
      END IF
      INS(2,ND1)=IW
      INS(3,ND1)=N
      INS(4,ND1)=JR
         WA=SQRT(WA)
      DO 140 I=INS(1,ND1),INS(2,ND1)
         U(I)=U(I)/WA
 140  CONTINUE   
 138  CONTINUE
      NI=NI-1
      IF(NI.EQ.0) GO TO 801
      GO TO 302
C
C    FIND A NEW PIVOT
C
  303 CONTINUE
C      WRITE(6,693)
C  693 FORMAT(' ZERO VECTOR')
  302 CONTINUE
      IB=IB+1
      IF(IB.LE.JDIM) GO TO 150
      WRITE(6,601)
  601 FORMAT(' STOP NEAR 302 IN TSTRLM')
      STOP
  801 CONTINUE
      NW=IW
      RETURN
      END
      SUBROUTINE TSJMDS(JJ,U,KP,INS)
      COMPLEX*16 U(100)
      DIMENSION KP(100),INS(4,24)
       DO 2 I=1,JJ+1
         WRITE(6,601) (INS(K,I),K=1,4),I
 601     FORMAT(2I5,'  N=',I2,'  JR=',I2,I5)
         DO 3 I1=INS(1,I),INS(2,I)
         WRITE(6,602) JJ-2*(KP(I1)-1),U(I1)
 602     FORMAT(' M=',I5,'/2 (',2F15.12,')')
 3       CONTINUE
 2       CONTINUE
      RETURN
      END
C SUBROUTINE SPINM4 ====*====3====*====4====*====5====*====6====*====7
C
C      ROTATION MATRIIS FOR J=3/2 SPACE ARE PREPARED
C      FOR CUBIC OR HEXAGONAL OPERATIONS
C
C      1989.12.12 :  A. YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
C
      SUBROUTINE SPINM4
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 J3HM
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      COMMON/JEQ3/J3HM(4,4,24)
      SAVE /SPG2/,/JEQ3/
      DIMENSION  EULER(3,24)
      INTEGER ICREAL(2,3,24)
      INTEGER IC(16),ICM(16,4,4),K1M(4,4),K2M(4,4),K3M(4,4),ISM(4,4)
C
      IF(IL.GE.1) THEN
         CALL EULERC(ICREAL,EULER)
         NN=24
      ELSE
         CALL EULERH(ICREAL,EULER)
         NN=12
      END IF
      IPR=0
      DO 31 J=1,4
      DO 31 K=1,4
         MA=5-2*J
         MB=5-2*K
         CALL TSRMHI(3,MA,MB,K1,K2,K3,IS,IC,IPR)
         DO 32 II=1,16
           ICM(II,J,K)=IC(II)
   32    CONTINUE
         K1M(J,K)=K1
         K2M(J,K)=K2
         K3M(J,K)=K3
         ISM(J,K)=IS
   31 CONTINUE
      DO 10 I=1,NN
C       WRITE(6,660) (EULER(K,I),K=1,3)
C 660   FORMAT(1H ,3F10.5)
        CB=COS(EULER(2,I))
        W1=COS(EULER(2,I)*0.5D0)
        W3=SIN(EULER(2,I)*0.5D0)
        DO 20 J=1,4
        DO 21 K=1,4
          MA=5-2*J
          MB=5-2*K
          WW=0.0
          DO 22 II=1,16
            IF(ICM(II,J,K).NE.0) THEN
            WW=WW+ICM(II,J,K)*(CB**(II-1))
            END IF
   22     CONTINUE
          WW=((WW*K1M(J,K))/K2M(J,K))*SQRT(DBLE(K3M(J,K)))
          IF(ISM(J,K).EQ.1) THEN
             WW=WW*W1
          ELSE
             WW=WW*W3
          END IF
          CA=COS(EULER(1,I)*0.5D0*MA)
          CC=COS(EULER(3,I)*0.5D0*MB)
          SA=SIN(EULER(1,I)*0.5D0*MA)
          SC=SIN(EULER(3,I)*0.5D0*MB)
          J3HM(J,K,I)=DCMPLX(CA,-SA)*DCMPLX(CC,-SC)*WW
          IF(IL.LE.0) THEN
             J3HM(J,K,I+12)=J3HM(J,K,I)
          END IF
   21   CONTINUE
C        WRITE(6,600) I,MA,(J3HM(J,K,I),K=1,4)
C  600   FORMAT(2I3,'/2',4(' (',2F10.7,')'))
C        WRITE(11) (J3HM(J,K,I),K=1,4)
   20   CONTINUE
C      WRITE(6,600)
   10 CONTINUE
      RETURN
      END
C SUBROUTINE SPNML4 ====*====3====*====4====*====5====*====6====*====7
C
C      MULITIPLICATIO TABLE OF SPIN HALF MATRIX
C      FOR CUBIC OR HEXAGONAL OPERATIONS
C
C      1988.10.10 :  A. YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
      SUBROUTINE SPNML4
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 J3HM
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      COMMON/JEQ3/J3HM(4,4,24)
      SAVE /SPG2/,/JEQ3/
      COMPLEX*16 WW,WM(4,4)
      INTEGER IDF(24,24)
      IF(IL.GE.1) NN=24
      IF(IL.LE.0) NN=12
      DO 10 I=1,NN
      DO 20 J=1,NN
C---
         DO 31 K1=1,4
         DO 33 K2=1,4
            WW=0.0
            DO 32 K=1,4
            WW=WW+J3HM(K1,K,I)*J3HM(K,K2,J)
   32       CONTINUE
         WM(K1,K2)=WW
   33    CONTINUE
C        WRITE(6,600) I,J,3-2*K1,(WM(K1,K2),K2=1,2)
C 600   FORMAT(2I3,I2,'/2',2(' (',2F10.7,')'))
   31    CONTINUE
C        WRITE(6,600)
C---
      DO 41 K=1,NN
         INIT=0
      DO 42 K1=1,4
      DO 43 K2=1,4
         IF(ABS(WM(K1,K2)).GT.0.000001D0) THEN
           IF(ABS(WM(K1,K2)-J3HM(K1,K2,K)).LT.0.000001D0) THEN
              INDS=1
           ELSE IF(ABS(WM(K1,K2)+J3HM(K1,K2,K)).LT.0.000001D0) THEN
              INDS=-1
           ELSE
              INDS=0
           END IF
           IF(INDS.EQ.0) GO TO 41
           IF(INIT.EQ.1.AND.INDS.NE.INDO) GO TO 41
           INDO=INDS
           INIT=1
         END IF
   43 CONTINUE
   42 CONTINUE
         IDF(I,J)=INDS*K
         IF(IL.LE.0) THEN
             IDF(I+12,J)=INDS*(K+12)
             IDF(I+12,J+12)=IDF(I,J)
             IDF(I,J+12)=INDS*(K+12)
         END IF
         GO TO 20
   41 CONTINUE
         WRITE(6,900) I,J
  900    FORMAT(' STOP AT 41 IN SPNMLT FOR',2I5)
         STOP
   20 CONTINUE
      WRITE(6,610) (IDF(I,J),J=1,NN)
  610 FORMAT(24I3)
   10 CONTINUE
      RETURN
      END