      FUNCTION AIJ(R,I,J)                                                       
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R(3)                                                            
C                                                                               
      GO TO (10,20,30),I                                                        
   10 GO TO (100,200,300),J                                                     
   20 GO TO (300,100,200),J                                                     
   30 GO TO (200,300,100),J                                                     
  100 AIJ=0.                                                                    
      RETURN                                                                    
C                                                                               
  200 MM=MOD(I+1,3)+1                                                           
      AIJ=R(MM)                                                                 
      RETURN                                                                    
C                                                                               
  300 MM=MOD(I,3)+1                                                             
      AIJ=-R(MM)                                                                
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION AXAXB(R,I,J,K,L,M,N,IP)                                          
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      AXAXB=AXA(R,I,J,R,K,L)*BIJK(R,M,N,IP)                                     
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION BIJK(R,I,JJ,KK)                                                  
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R(3)                                                            
C                                                                               
      DIV=2.                                                                    
      IF(JJ.EQ.KK)DIV=1.                                                        
      J=MIN0(JJ,KK)                                                             
      K=MAX0(JJ,KK)                                                             
      GO TO (10,20,30),I                                                        
   10 GO TO (12,14,400),J                                                       
   12 GO TO (100,200,300),K                                                     
   14 IF(K-2)400,400,100                                                        
   20 GO TO (22,24,400),J                                                       
   22 GO TO (400,300,100),K                                                     
   24 IF(K-2)100,100,200                                                        
   30 GO TO (32,34,100),J                                                       
   32 GO TO (400,100,200),K                                                     
   34 IF(K-2)400,400,300                                                        
  100 BIJK=0.                                                                   
      RETURN                                                                    
C                                                                               
  200 MM=MOD(I,3)+1                                                             
      BIJK=.5*R(MM)/DIV                                                         
      RETURN                                                                    
C                                                                               
  300 MM=MOD(I+1,3)+1                                                           
      BIJK=.5*R(MM)/DIV                                                         
      RETURN                                                                    
C                                                                               
  400 BIJK=-.5*R(I)/DIV                                                         
      RETURN                                                                    
      END                                                                       
      FUNCTION CIJKL(R,I,JJ,KK,LL)                                              
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R(3)                                                            
C                                                                               
      IF(JJ.EQ.KK.AND.JJ.EQ.LL)GO TO 7                                          
      IF(JJ.NE.KK.AND.JJ.NE.LL.AND.KK.NE.LL)GO TO 8                             
      IF(JJ-KK)1,2,3                                                            
    1 J=JJ                                                                      
      L=KK                                                                      
      K=LL                                                                      
      GO TO 4                                                                   
    2 J=MIN0(JJ,LL)                                                             
      K=KK                                                                      
      L=MAX0(JJ,LL)                                                             
      GO TO 4                                                                   
    3 J=KK                                                                      
      K=LL                                                                      
      L=JJ                                                                      
    4 DIV=3.                                                                    
      GO TO 9                                                                   
    7 J=JJ                                                                      
      K=JJ                                                                      
      L=JJ                                                                      
      DIV=1.                                                                    
      GO TO 9                                                                   
    8 J=1                                                                       
      K=2                                                                       
      L=3                                                                       
      DIV=6.                                                                    
    9 GO TO (10,20,30),I                                                        
   10 GO TO (12,15,300),J                                                       
   12 GOTO (13,100,100),K                                                       
   13 GO TO (100,200,300),L                                                     
   15 IF(K-2)16,16,200                                                          
   16 IF(L-2)200,200,300                                                        
   20 GO TO (22,25,200),J                                                       
   22 GOTO (23,24,300),K                                                        
   23 GO TO (300,100,200),L                                                     
   24 IF(L-2)300,300,100                                                        
   25 IF(K-2)26,26,100                                                          
   26 IF(L-2)100,100,200                                                        
   30 GO TO (32,35,100),J                                                       
   32 GO TO (33,34,200),K                                                       
   33 GO TO (200,300,100),L                                                     
   34 IF(L-2)200,200,100                                                        
   35 IF(K-2)36,36,300                                                          
   36 IF(L-2)300,300,100                                                        
  100 CIJKL=0.                                                                  
      RETURN                                                                    
C                                                                               
  200 MM=MOD(I+1,3)+1                                                           
      CIJKL=-R(MM)/(6.*DIV)                                                     
      RETURN                                                                    
C                                                                               
  300 MM=MOD(I,3)+1                                                             
      CIJKL=R(MM)/(6.*DIV)                                                      
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      SUBROUTINE CORECT(R3,R1,R2,ELT,THETA,RC,BIGR,BIGRP,RP2)                   
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R3(3),R1(3),R2(3),RC(3),BIGR(3),BIGRP(3),ELT(3,3),              
     & RP(3),RPP(3),ACON(3,3),BCON(3)                                           
C                                                                               
      RP2=0.                                                                    
      DO 10 I=1,3                                                               
      RP(I)=R2(I)-R1(I)                                                         
   10 RP2=RP2+RP(I)**2                                                          
      ISIG=1                                                                    
      DO 20 I53K=1,3                                                            
      I=I53K                                                                    
      BCON(I)=RC(I)                                                             
      DO 20 J53K=1,3                                                            
      J=J53K                                                                    
      ACON(I,J)=FIJ(THETA,ELT,RP,RP2,I,J)                                       
   20 BCON(I)=BCON(I)+ACON(I,J)*R2(J)                                           
      GO TO 40                                                                  
C                                                                               
C---- ENTRY POINT ---- ENTRY POINT ----                                         
C                                                                               
      ENTRY CORLIB(R3,ELT,RC)                                                   
C                                                                               
C---- ENTRY POINT ---- ENTRY POINT ----                                         
C                                                                               
      DO 30 I=1,3                                                               
      BCON(I)=RC(I)                                                             
      DO 30 J=1,3                                                               
   30 ACON(I,J)=0.                                                              
      ISIG=0                                                                    
   40 DO 50 I1=1,3                                                              
      DO 50 J1=1,3                                                              
      I=I1                                                                      
      J=J1                                                                      
   50 ACON(I,J)=ACON(I,J)+EIJ(ELT,I,J)                                          
      CALL SOLVE(ACON,BCON,R3)                                                  
      IF(ISIG.GT.0)GO TO 60                                                     
      RETURN                                                                    
C                                                                               
   60 DO 70 I=1,3                                                               
   70 RPP(I)=R3(I)-R2(I)                                                        
      CALL VXV(RP,RPP,BIGR)                                                     
      CALL VXV(RP,BIGR,BIGRP)                                                   
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION DETERM(C1,C2,C3)                                                 
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION C1(3),C2(3),C3(3)                                               
C                                                                               
      DETERM=C1(1)*C2(2)*C3(3)+C1(2)*C2(3)*C3(1)+C1(3)*C2(1)*C3(2)              
     & -C1(3)*C2(2)*C3(1)-C1(2)*C2(1)*C3(3)-C1(1)*C2(3)*C3(2)                   
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION DIJKLM(R,I,JI,KI,LI,MI)                                          
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R(3)                                                            
C                                                                               
      IN=JI*KI*LI*MI                                                            
      IF(MOD(IN,6).EQ.0.AND.IN.LE.18)GO TO 1080                                 
      IF(IN.EQ.4.OR.IN.EQ.9.OR.IN.EQ.36)GO TO 1060                              
      IF(JI.NE.KI.OR.KI.NE.LI.OR.LI.NE.MI)GO TO 1040                            
      DIV=1.                                                                    
      GO TO 1100                                                                
 1040 DIV=4.                                                                    
      GO TO 1100                                                                
 1060 DIV=6.                                                                    
      GO TO 1100                                                                
 1080 DIV=12.                                                                   
 1100 IF(MOD(IN,3).NE.0)GO TO 1500                                              
      M=3                                                                       
      IN=IN/3                                                                   
      IF(MOD(IN,3).NE.0)GO TO 1250                                              
      L=3                                                                       
      IN=IN/3                                                                   
      IF(MOD(IN,3).NE.0)GO TO 1125                                              
      K=3                                                                       
      J=IN/3                                                                    
      GO TO 1900                                                                
 1125 IF(MOD(IN,2).NE.0)GO TO 1150                                              
      K=2                                                                       
      J=IN/2                                                                    
      GO TO 1900                                                                
 1150 K=1                                                                       
      J=1                                                                       
      GO TO 1900                                                                
 1250 IF(MOD(IN,2).NE.0)GO TO 1375                                              
      L=2                                                                       
      IN=IN/2                                                                   
      IF(MOD(IN,2).NE.0)GO TO 1150                                              
      K=2                                                                       
      J=IN/2                                                                    
      GO TO 1900                                                                
 1375 L=1                                                                       
      GO TO 1150                                                                
 1500 IF(MOD(IN,2).NE.0)GO TO 1750                                              
      M=2                                                                       
      IN=IN/2                                                                   
      GO TO 1250                                                                
 1750 M=1                                                                       
      GO TO 1375                                                                
 1900 GO TO (100,200,300),I                                                     
  100 GO TO (110,150,2000),J                                                    
  110 GO TO (120,140,3000),K                                                    
  120 GO TO (125,130,2000),L                                                    
  125 GO TO (1000,4000,3000),M                                                  
  130 IF(M-2)1000,2000,1000                                                     
  140 IF(L-2)1000,145,4000                                                      
  145 IF(M-2)1000,4000,3000                                                     
  150 IF(K-2)1000,160,1000                                                      
  160 IF(L-2)1000,170,5000                                                      
  170 IF(M-2)1000,2000,1000                                                     
  200 GO TO (210,250,2000),J                                                    
  210 GO TO (220,240,1000),K                                                    
  220 GO TO (225,230,5000),L                                                    
  225 GO TO (2000,3000,1000),M                                                  
  230 IF(M-2)1000,2000,4000                                                     
  240 IF(L-2)1000,245,3000                                                      
  245 IF(M-2)1000,3000,1000                                                     
  250 IF(K-2)1000,260,4000                                                      
  260 IF(L-2)1000,270,2000                                                      
  270 IF(M-2)1000,1000,4000                                                     
  300 GO TO (310,350,1000),J                                                    
  310 GO TO (320,340,4000),K                                                    
  320 GO TO (325,330,2000),L                                                    
  325 GO TO (2000,1000,4000),M                                                  
  330 IF(M-2)1000,5000,3000                                                     
  340 IF(L-2)1000,345,1000                                                      
  345 IF(M-2)1000,1000,4000                                                     
  350 IF(K-2)1000,360,3000                                                      
  360 IF(L-2)1000,370,2000                                                      
  370 IF(M-2)1000,2000,3000                                                     
 1000 DIJKLM=0.                                                                 
      RETURN                                                                    
C                                                                               
 2000 FAC=1.                                                                    
 2010 DIJKLM=FAC*R(I)/(24.*DIV)                                                 
      RETURN                                                                    
C                                                                               
 3000 MM=MOD(I+1,3)+1                                                           
 3010 DIJKLM=-R(MM)/(24.*DIV)                                                   
      RETURN                                                                    
C                                                                               
 4000 MM=MOD(I,3)+1                                                             
      GO TO 3010                                                                
 5000 FAC=2.                                                                    
      GO TO 2010                                                                
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION D2KDL(R,EL,S,I,J,K,L)                                            
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R(3),EL(3,3),S(3,3)                                             
C                                                                               
      C(I,J,K,L)=CIJKL(R,I,J,K,L)                                               
      C6(I,K,L,M)=6.*C(I,K,L,M)                                                 
      ACC(I,J,K,L,M,N)=3.*AIJ(R,I,K)*C(J,L,M,N)                                 
C                                                                               
      SUM=AIJ(R,I,K)*AIJ(R,J,L)+AIJ(R,I,L)*AIJ(R,J,K)                           
      DO 1000 MXT=1,3                                                           
      M=MXT                                                                     
      SUM=SUM+S(M,I)*C6(J,K,L,M)+S(M,J)*C6(I,K,L,M)                             
      DO 1000 NXY=1,3                                                           
      N=NXY                                                                     
      SUM=SUM+EL(M,N)*(ACC(I,J,K,L,M,N)+ACC(I,J,L,K,M,N)+ACC(J,I,K,L,M,N        
     &      )+ACC(J,I,L,K,M,N)+ACC(I,J,M,N,K,L)+ACC(J,I,M,N,K,L)+ACC(I,J        
     &      ,M,N,L,K)+ACC(J,I,M,N,L,K)+4.*(BIJK(R,I,K,M)*BIJK(R,J,L,N)+         
     &      BIJK(R,I,L,M)*BIJK(R,J,K,N)))                                       
 1000 CONTINUE                                                                  
      D2KDL=SUM                                                                 
      IF(K.EQ.L)D2KDL=SUM/2.                                                    
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION D2KDS(R,EL,I,J,K,L)                                              
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R(3),EL(3,3)                                                    
C                                                                               
      CCC(J,K,L,M)=3.*CIJKL(R,J,K,L,M)                                          
C                                                                               
      SUM=0.                                                                    
      IF(I.NE.L)GO TO 1000                                                      
      SUM=SUM+AIJ(R,J,K)                                                        
      DO 500 MXV=1,3                                                            
      DO 500 NXV=1,3                                                            
      M=MXV                                                                     
      N=NXV                                                                     
  500 SUM=SUM+EL(M,N)*CCC(J,K,M,N)                                              
 1000 IF(J.NE.L)GO TO 2000                                                      
      SUM=SUM+AIJ(R,I,K)                                                        
      DO 1500 MTY=1,3                                                           
      DO 1500 NTY=1,3                                                           
      M=MTY                                                                     
      N=NTY                                                                     
 1500 SUM=SUM+EL(M,N)*CCC(I,K,M,N)                                              
 2000 D2KDS=SUM                                                                 
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION D2KDT(I,J,K,L)                                                   
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      D2KDT=0.                                                                  
      IF(I.EQ.K.AND.J.EQ.L)D2KDT=1.                                             
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION D3KDL(R,EL,S,I,J,K,L,M)                                          
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      INTEGER P,PX                                                              
C                                                                               
      DIMENSION R(3),EL(3,3),S(3,3)                                             
C                                                                               
      ABB(J,K,L,M,N)=2.*AIJ(R,J,L)*BIJK(R,K,M,N)                                
      AAB(I,J,K,L,M,N,NN)=AIJ(R,I,L)*AIJ(R,J,N)*BIJK(R,K,M,NN)                  
      PP(I,J,K,L,M)=PIJKLM(R,I,J,K,L,M)+PIJKLM(R,I,J,K,M,L)                     
C                                                                               
      SUM=0.                                                                    
      DO 1000 NHJ=1,3                                                           
      N=NHJ                                                                     
      SUM=SUM+S(N,I)*(ABB(J,K,L,M,N)+ABB(K,J,L,M,N)+ABB(J,K,M,L,N)+ABB(K        
     &    ,J,M,L,N))+S(N,J)*(ABB(I,K,L,M,N)+ABB(K,I,L,M,N)+ABB(I,K,M,L,N        
     &    )+ABB(K,I,M,L,N))+S(N,K)*(ABB(I,J,L,M,N)+ABB(J,I,L,M,N)+ABB(I,        
     &    J,M,L,N)+ABB(J,I,M,L,N))                                              
      DO 1000 PX=1,3                                                            
      P=PX                                                                      
      SUM=SUM+S(N,J)*S(P,K)*PP(I,N,P,L,M)+S(N,I)*S(P,K)*PP(J,N,P,L,M)           
     &    +S(N,I)*S(P,J)*PP(K,N,P,L,M)+2.*EL(N,P)*(AAB(I,J,K,L,M,N,P)           
     &    +AAB(I,K,J,L,M,N,P)+AAB(J,K,I,L,M,N,P)+AAB(I,J,K,N,P,L,M)             
     &    +AAB(I,K,J,N,P,L,M)+AAB(J,K,I,N,P,L,M))                               
 1000 CONTINUE                                                                  
      D3KDL=SUM                                                                 
      IF(L.EQ.M)D3KDL=SUM/2.                                                    
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION D3KDS(R,EL,S,I,J,K,LP,MP)                                        
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      INTEGER P,PHY                                                             
C                                                                               
      DIMENSION R(3),EL(3,3),S(3,3)                                             
C                                                                               
      BB(K,L,M)=2.*BIJK(R,K,L,M)                                                
      ABB(J,K,L,M,N)=2.*AIJ(R,J,L)*BIJK(R,K,M,N)                                
C                                                                               
      SUM=0.                                                                    
      DO 1000 LHY=1,3                                                           
      L=LHY                                                                     
      DO 1000 MHY=1,3                                                           
      M=MHY                                                                     
      IF(L.EQ.LP.AND.I.EQ.MP)SUM=SUM+S(M,J)*BB(K,L,M)+S(M,K)*BB(J,L,M)          
      IF(M.EQ.LP.AND.J.EQ.MP)SUM=SUM+S(L,I)*BB(K,L,M)                           
      IF(M.EQ.LP.AND.K.EQ.MP)SUM=SUM+S(L,I)*BB(J,L,M)+S(L,J)*BB(I,L,M)          
      IF(L.EQ.LP.AND.J.EQ.MP)SUM=SUM+S(M,K)*BB(I,L,M)                           
      DO 1000 NHY=1,3                                                           
      N=NHY                                                                     
      IF(L.NE.LP)GO TO 500                                                      
      IF(I.EQ.MP)SUM=SUM+EL(M,N)*(ABB(J,K,M,L,N)+ABB(K,J,M,L,N))                
      IF(J.EQ.MP)SUM=SUM+EL(M,N)*(ABB(I,K,M,L,N)+ABB(K,I,M,L,N))                
      IF(K.EQ.MP)SUM=SUM+EL(M,N)*(ABB(I,J,M,L,N)+ABB(J,I,M,L,N))                
 500  DO 1000 PHY=1,3                                                           
      P=PHY                                                                     
      IF(L.NE.LP)GO TO 750                                                      
      IF(J.EQ.MP)SUM=SUM+EL(N,P)*S(M,K)*PIJKLM(R,I,L,M,N,P)                     
      IF(I.EQ.MP)SUM=SUM+EL(N,P)*(S(M,K)*PIJKLM(R,J,L,M,N,P)+S(M,J)*            
     &               PIJKLM(R,K,L,M,N,P))                                       
 750  IF(M.NE.LP)GO TO 1000                                                     
      IF(K.EQ.MP)SUM=SUM+EL(N,P)*(S(L,J)*PIJKLM(R,I,L,M,N,P)+S(L,I)*            
     &               PIJKLM(R,J,L,M,N,P))                                       
      IF(J.EQ.MP)SUM=SUM+EL(N,P)*S(L,I)*PIJKLM(R,K,L,M,N,P)                     
 1000 CONTINUE                                                                  
      D3KDS=SUM                                                                 
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION EIJ(EL,I,J)                                                      
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION EL(3,3)                                                         
C                                                                               
      IF(I.EQ.J)GO TO 20                                                        
      EIJ=.5*EL(I,J)-(1./24.)*(EL4(EL,I,1,1,J)+EL4(EL,I,2,2,J)+EL4(EL,          
     & I,3,3,J))                                                                
      RETURN                                                                    
C                                                                               
   20 K=MOD(I,3)+1                                                              
      L=MOD(I+1,3)+1                                                            
      EIJ=1.-.5*(EL(K,K)+EL(L,L))+(1./24.)*(EL4(EL,I,I,K,K)+EL4(EL,I,I,         
     & L,L)+3.*(EL(K,K)**2+EL(L,L)**2)+2.*EL4(EL,K,K,L,L))                      
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION EL4(EL,I,J,K,L)                                                  
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION EL(3,3)                                                         
      EL4=EL(I,J)*EL(K,L)+EL(I,K)*EL(J,L)+EL(I,L)*EL(J,K)                       
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
                  SUBROUTINE EULER(PHI,CHI,OMEGA,TRANS)                         
C                                                                               
C     SUBROUTINE TO GENERATE THE MATRIX B-INVERSE FROM THE EULERIAN             
C                      ANGLES PHI,CHI,AND OMEGA.                                
C                                                                               
C                            BY L. W. FINGER                                    
C                                                                               
      DIMENSION TRANS(3,3)                                                      
C                                                                               
      DATA RAD/57.295779/                                                       
      ASSIGN 10 TO KGO                                                          
      GO TO 100                                                                 
   10 TRANS(1,1)=CP*CX*CO-SP*SO                                                 
      TRANS(1,2)=-CP*CX*SO-SP*CO                                                
      TRANS(1,3)=CP*SX                                                          
      TRANS(2,1)=SP*CX*CO+CP*SO                                                 
      TRANS(2,2)=-SP*CX*SO+CP*CO                                                
      TRANS(2,3)=SP*SX                                                          
      TRANS(3,1)=-SX*CO                                                         
      TRANS(3,2)=SX*SO                                                          
      TRANS(3,3)=CX                                                             
      RETURN                                                                    
C                                                                               
C                                                                               
      ENTRY DTDPHI(PHI,CHI,OMEGA,TRANS)                                         
C                                                                               
C                                                                               
      ASSIGN 20 TO KGO                                                          
      GO TO 100                                                                 
   20 TRANS(1,1)=-SP*CX*CO-CP*SO                                                
      TRANS(1,2)=SP*CX*SO-CP*CO                                                 
      TRANS(1,3)=-SP*SX                                                         
      TRANS(2,1)=CP*CX*CO-SP*SO                                                 
      TRANS(2,1)=-CP*CX*SO-SP*CO                                                
      TRANS(2,3)=CP*SX                                                          
      DO 25 I=1,3                                                               
      DO 23 J=1,2                                                               
   23 TRANS(J,I)=TRANS(J,I)/RAD                                                 
   25 TRANS(3,I)=0.                                                             
      RETURN                                                                    
C                                                                               
C                                                                               
      ENTRY DTDCHI(PHI,CHI,OMEGA,TRANS)                                         
C                                                                               
C                                                                               
      ASSIGN 30 TO KGO                                                          
      GO TO 100                                                                 
   30 TRANS(1,1)=-CP*SX*CO                                                      
      TRANS(1,2)=CP*SX*SO                                                       
      TRANS(1,3)=CP*CX                                                          
      TRANS(2,1)=-SP*SX*CO                                                      
      TRANS(2,2)=SP*SX*SO                                                       
      TRANS(2,3)=SP*CX                                                          
      TRANS(3,1)=-CX*CO                                                         
      TRANS(3,2)=CX*SO                                                          
      TRANS(3,3)=-SX                                                            
      DO 35 I=1,3                                                               
      DO 35 J=1,3                                                               
   35 TRANS(I,J)=TRANS(I,J)/RAD                                                 
      RETURN                                                                    
C                                                                               
C                                                                               
      ENTRY DTDOMG(PHI,CHI,OMEGA,TRANS)                                         
C                                                                               
C                                                                               
      ASSIGN 40 TO KGO                                                          
      GO TO 100                                                                 
   40 TRANS(1,1)=-CP*CX*SO-SP*CO                                                
      TRANS(1,2)=-CP*CX*CO+SP*SO                                                
      TRANS(2,1)=-SP*CX*SO+CP*CO                                                
      TRANS(2,2)=-SP*CX*CO-CP*SO                                                
      TRANS(3,1)=SX*SO                                                          
      TRANS(3,2)=SX*CO                                                          
      DO 45 I=1,3                                                               
      DO 43 J=1,2                                                               
   43 TRANS(I,J)=TRANS(I,J)/RAD                                                 
   45 TRANS(I,3)=0.                                                             
      RETURN                                                                    
C                                                                               
  100 CP=COS(PHI/RAD)                                                           
      CX=COS(CHI/RAD)                                                           
      CO=COS(OMEGA/RAD)                                                         
      SP=SIN(PHI/RAD)                                                           
      SX=SIN(CHI/RAD)                                                           
      SO=SIN(OMEGA/RAD)                                                         
      GO TO KGO,(10,20,30,40)                                                   
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION FIJ(THETA,EL,RP,RP2,I,J)                                         
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION EL(3,3),RP(3),CON(3)                                            
C                                                                               
      IF(RP2.LT.0.01.OR.THETA.LT.0.000001)GO TO 100                             
      K=MOD(I,3)+1                                                              
      L=MOD(I+1,3)+1                                                            
      KK=MOD(J,3)+1                                                             
      LL=MOD(J+1,3)+1                                                           
      M2=MOD(J-I+2,3)+1                                                         
      M3=MOD(J-I+3,3)+1                                                         
      M4=MOD(J-I+4,3)+1                                                         
      CON1=1.-THETA/4.-.5*(EL(K,K)+EL(L,L))                                     
      CON2=.5*EL(I,K)                                                           
      CON3=.5*EL(I,L)                                                           
      CON(M3)=-(RP(KK)**2+RP(LL)**2)                                            
      CON(M4)=RP(J)*RP(KK)                                                      
      CON(M2)=RP(J)*RP(LL)                                                      
      FIJ=(CON1*CON(1)+CON2*CON(2)+CON3*CON(3))*THETA/(2.*RP2)                  
      RETURN                                                                    
C                                                                               
  100 FIJ=0.                                                                    
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
                      SUBROUTINE FILLIN(IK,IP,P,NA)                             
C                                                                               
C         SUBROUTINE TO FILL IN THE SPECIAL PARAMETERS FROM THE LIST.           
C     IP=TYPE OF PARAMETER, P=CURRENT VALUE, IK IS POSITION IN PARAMETER        
C                                  LIST.                                        
C                                                                               
C                             BY L. W. FINGER                                   
C                                                                               
C                                                                               
      DIMENSION SAVE(60),KDEGN(2640),AA(2640)                                   
C                                                                               
      COMMON/XPAR/NPAR,ATOX(3,3),XTOA(3,3),PAR(60),PARNAM(2,60),IPAR(60)        
     & ,IRPAR(60),OXYZ(3),SXYZ(3),OXP(3),T(3,3),EL(3,3),S(3,3),TRNS(3,3)        
     &  ,TTRNS(3,3), TRANSI(3,3), TTRNSI(3,3), RTEM(3), BTEM(6),GTEM(10)        
      COMMON/LSMAT/ B(200),IEXT,IPARA(1000),MRANK,KCYL,LABEL(3),JSCRF(10        
     & ), NRSF, SIGMA, MDEGN, MAXIMA, IALFA(60), RDAMP(200), IDAMP(200),        
     & A(2200)                                                                  
      COMMON/IND23/IND2(2,6),IND3(3,10)                                         
      COMMON/D/ACALC,AIND(60),ANOM(2,20),BCALC,BIND(60),DADB(6,60),DADO(        
     & 60),DADX(3,60),DBDB(6,60),DBDO(60),DBDX(3,60),H(3),IH(3),ISCAT(2,        
     & 60), ISCL, MFSEL, OCCA(60), OCCUP(60),SCALE(10),SCAT(20),SITE(60)        
     & ,DADC(10,60), DBDC(10,60), DADD(15,60), DBDD(15,60),AMINUS,BMINUS        
     & ,DRDE(7),NDEGN,VORFAK,IPAREX(7),DAALFA(60),DBALFA(60)                    
C                                                                               
      EQUIVALENCE (SAVE,AIND),(KDEGN,AA,DADC(1,37))                             
C                                                                               
      I1(K)=MOD(K/256**3,256)                                                   
      I2(K)=MOD(K/256**2,256)                                                   
      I3(K)=MOD(K/256,256)                                                      
      I4(K)=MOD(K,256)                                                          
C                                                                               
      JD=(IP+2)/3                                                               
      GO TO(100,200,300,400,400,500,500,600,600,600),JD                         
C---- SPECIAL ORIGIN                                                            
  100 OXYZ(IP)=P                                                                
      RETURN                                                                    
C                                                                               
C---- ORIENTATION                                                               
  200 OXP(IP-3)=P                                                               
      RETURN                                                                    
C                                                                               
C---- SPECIAL XYZ. WE MUST DETERMINE IF IT APPLIES TO ATOM NA.                  
  300 LDEGN=2*NDEGN                                                             
      DO 310 I=1,LDEGN,2                                                        
      IF(I2(KDEGN(I)).NE.IK)GO TO 310                                           
      IF(I1(KDEGN(I)).NE.0)GO TO 310                                            
      IF(I3(KDEGN(I)).NE.NA) GO TO 310                                          
      GO TO 320                                                                 
  310 CONTINUE                                                                  
      RETURN                                                                    
C                                                                               
  320 SXYZ(IP-6)=P                                                              
      RETURN                                                                    
C                                                                               
C---- T TENSOR ELEMENT                                                          
  400 I=IP-9                                                                    
      J=IND2(1,I)                                                               
      K=IND2(2,I)                                                               
      T(J,K)=P                                                                  
      IF(J.NE.K)T(K,J)=P                                                        
      RETURN                                                                    
C                                                                               
C---- L TENSOR ELEMENT                                                          
  500 I=IP-15                                                                   
      J=IND2(1,I)                                                               
      K=IND2(2,I)                                                               
      EL(J,K)=P                                                                 
      IF(J.NE.K)EL(K,J)=P                                                       
      RETURN                                                                    
C                                                                               
C---- ELEMENT OF S TENSOR                                                       
  600 I=IP-21                                                                   
      IF(I.GT.6)GO TO 605                                                       
      J=IND2(1,I)                                                               
      K=IND2(2,I)                                                               
      GO TO 610                                                                 
  605 K=IND2(1,I-3)                                                             
      J=IND2(2,I-3)                                                             
  610 S(J,K)=P                                                                  
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION SKTLS(R,I,J,T,EL,S)                                              
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      DIMENSION R(3),T(3,3),EL(3,3),S(3,3)                                      
C                                                                               
      SUM=T(I,J)                                                                
      DO 1000 K23=1,3                                                           
      K=K23                                                                     
      SUM=SUM+AIJ(R,J,K)*S(K,I)+AIJ(R,I,K)*S(K,J)                               
      DO 1000 L23=1,3                                                           
      L=L23                                                                     
      SUM=SUM+AIJ(R,I,K)*AIJ(R,J,L)*EL(K,L)                                     
      DO 1000 M23=1,3                                                           
      M=M23                                                                     
      SUM=SUM+3.*EL(L,M)*(S(K,I)*CIJKL(R,J,K,L,M)+S(K,J)*CIJKL(R,I,K,L,M        
     & ))                                                                       
      DO 1000 N23=1,3                                                           
      N=N23                                                                     
      SUM=SUM+EL(K,L)*EL(M,N)*(3.*(AIJ(R,I,K)*CIJKL(R,J,L,M,N)+AIJ(R,J,K        
     & )*CIJKL(R,I,L,M,N))+2.*BIJK(R,I,K,M)*BIJK(R,J,L,N))                      
 1000 CONTINUE                                                                  
      SKTLS=SUM                                                                 
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
                          SUBROUTINE SOLVE(A,B,C)                               
C                                                                               
C                             BY L. W. FINGER                                   
C                                                                               
      DIMENSION A(3,3),B(3),C(3)                                                
C                                                                               
      DELTA=DETERM(A(1,1),A(1,2),A(1,3))                                        
      C(1)=DETERM(B(1),A(1,2),A(1,3))/DELTA                                     
      C(2)=DETERM(A(1,1),B(1),A(1,3))/DELTA                                     
      C(3)=DETERM(A(1,1),A(1,2),B(1))/DELTA                                     
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
                          SUBROUTINE SPVAL                                      
C                                                                               
C               DUMMY SUBROUTINES SPVAL AND SPDERI                              
      RETURN                                                                    
C                                                                               
C                                                                               
      ENTRY SPDERI                                                              
C                                                                               
C                                                                               
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
                       SUBROUTINE TFM2(T,A,B)                                   
C                                                                               
C                          BY L. W. FINGER                                      
C                                                                               
      DIMENSION A(6),B(6),C(10),D(10),T(3,3),BTEM(3,3),GTEM(3,3,3),J(3)         
C                                                                               
      COMMON/IND23/IND2(2,6),IND3(3,10)                                         
C                                                                               
      EQUIVALENCE(BTEM(1,1),GTEM(1,1,1)),(J1,J(1)),(J2,J(2)),(J3,J(3))          
C                                                                               
      DATA PI/3.1415926/                                                        
C                                                                               
      DO 10 I=1,6                                                               
      DO 5 K=1,2                                                                
    5 J(K)=IND2(K,I)                                                            
      BTEM(J1,J2)=A(I)                                                          
      BTEM(J2,J1)=A(I)                                                          
   10 CONTINUE                                                                  
      DO 20 I=1,6                                                               
      B(I)=0.                                                                   
      DO 15 K=1,2                                                               
   15 J(K)=IND2(K,I)                                                            
      DO 18 K=1,3                                                               
      DO 18 L=1,3                                                               
   18 B(I)=B(I)+T(J1,K)*T(J2,L)*BTEM(K,L)                                       
   20 B(I)=B(I)*2.*PI**2                                                        
      RETURN                                                                    
C                                                                               
C                                                                               
      ENTRY TFM3(T,C,D)                                                         
C                                                                               
C                                                                               
      DO 30 I=1,10                                                              
      DO 25 K=1,3                                                               
   25 J(K)=IND3(K,I)                                                            
      GTEM(J1,J2,J3)=C(I)                                                       
      GTEM(J2,J3,J1)=C(I)                                                       
      GTEM(J3,J1,J2)=C(I)                                                       
      GTEM(J3,J2,J1)=C(I)                                                       
      GTEM(J2,J1,J3)=C(I)                                                       
   30 GTEM(J1,J3,J2)=C(I)                                                       
      DO 40 I=1,10                                                              
      D(I)=0.                                                                   
      DO 35 K=1,3                                                               
   35 J(K)=IND3(K,I)                                                            
      DO 38 K=1,3                                                               
      DO 38 L=1,3                                                               
      DO 38 M=1,3                                                               
   38 D(I)=D(I)+T(J1,K)*T(J2,L)*T(J3,M)*GTEM(K,L,M)                             
   40 D(I)=D(I)*4.*PI**3/3.                                                     
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION TKTLS(R,I,J,K,EL,S)                                              
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      INTEGER P,P12                                                             
C                                                                               
      DIMENSION R(3),EL(3,3),S(3,3)                                             
C                                                                               
      B(K,L,M)=BIJK(R,K,L,M)                                                    
      A(I,J)=AIJ(R,I,J)                                                         
C                                                                               
      SUM=0.                                                                    
      DO 1000 L12=1,3                                                           
      L=L12                                                                     
      DO 1000 M12=1,3                                                           
      M=M12                                                                     
      SUM=SUM+2.*(S(L,I)*S(M,J)*B(K,L,M)+S(L,I)*S(M,K)*B(J,L,M)+S(L,J)*         
     & S(M,K)*B(I,L,M))                                                         
      DO 1000 N12=1,3                                                           
      N=N12                                                                     
      SUM=SUM+2.*EL(M,N)*(S(L,I)*(A(J,M)*B(K,L,N)+A(K,M)*B(J,L,N))+S(L,J        
     & )*(A(I,M)*B(K,L,N)+A(K,M)*B(I,L,N))+S(L,K)*(A(I,M)*B(J,L,N)+             
     & A(J,M)*B(I,L,N)))                                                        
      DO 1000 P12=1,3                                                           
      P=P12                                                                     
      SUM=SUM+EL(N,P)*(S(L,J)*S(M,K)*PIJKLM(R,I,L,M,N,P)+S(L,I)*S(M,K)*         
     & PIJKLM(R,J,L,M,N,P)+S(L,I)*S(M,J)*PIJKLM(R,K,L,M,N,P)+2.*EL(L,           
     & M)*(A(I,L)*A(J,N)*B(K,M,P)+A(I,L)*A(K,N)*B(J,M,P)+A(J,L)*A(K,N           
     & )*B(I,M,P)))                                                             
 1000 CONTINUE                                                                  
      TKTLS=SUM                                                                 
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
                           SUBROUTINE VXV(V1,V2,V3)                             
C                                                                               
C                               BY L. F. FINGER                                 
C                                                                               
      DIMENSIONV1(3),V2(3),V3(3)                                                
C                                                                               
      V3(1)=V1(2)*V2(3)-V1(3)*V2(2)                                             
      V3(2)=V1(3)*V2(1)-V1(1)*V2(3)                                             
      V3(3)=V1(1)*V2(2)-V1(2)*V2(1)                                             
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION AXA(R1,I,J,R2,K,L)                                               
C                                                                               
C     BY L. W. FINGER                                                           
C                                                                               
      AXA=AIJ(R1,I,J)*AIJ(R2,K,L)                                               
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
      FUNCTION PIJKLM(R,I,J,K,L,M)                                              
C                                                                               
C     BY L. F. FINGER                                                           
C                                                                               
      DIMENSION R(3)                                                            
C                                                                               
      D(J,K,L,M)=DIJKLM(R,I,J,K,L,M)                                            
      PIJKLM=12.*D(J,K,L,M)                                                     
      RETURN                                                                    
C                                                                               
      END                                                                       
                        SUBROUTINE TRNSPS(A,B)                                  
C                                                                               
C                  TRANSPOSE MATRIX A AND STORE IN B.                           
C                                                                               
C                           BY L. W. FINGER                                     
C                                                                               
      DIMENSION A(3,3),B(3,3)                                                   
C                                                                               
      DO 100 I=1,3                                                              
      DO 100 J=1,3                                                              
 100  B(J,I)=A(I,J)                                                             
      RETURN                                                                    
C                                                                               
      END                                                                       
C                                                                               
                          SUBROUTINE MULMM(A,B,C)                               
C                                                                               
C                   MULTIPLY 2 3X3 MATRICES TOGETHER.                           
C                                                                               
C                            BY L. W. FINGER                                    
C                                                                               
      DIMENSION A(3,3),B(3,3),C(3,3),D(3),E(3)                                  
C                                                                               
      DO 10 I=1,3                                                               
      DO 10 J=1,3                                                               
      C(I,J)=0.                                                                 
      DO 10 K=1,3                                                               
 10   C(I,J)=C(I,J)+A(I,K)*B(K,J)                                               
      RETURN                                                                    
C                                                                               
      ENTRY MULMV(A,D,E)                                                        
C                                                                               
C                                                                               
      DO 20 I=1,3                                                               
      E(I)=0.                                                                   
      DO 20 J=1,3                                                               
 20   E(I)=E(I)+A(I,J)*D(J)                                                     
      RETURN                                                                    
C                                                                               
      END                                                                       
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
