       IMPLICIT DOUBLE PRECISION(A-H,O-Z)                               
C                                                                       
C                                                                       
C                                                                       
C  A CODE TO CALCULATE ANGULAR MOMENTUM  COUPLING COEFFICIENTS          
C                                                                       
C   VERSION 3.0 , BY N M CLARKE, UNIVERSITY OF BIRMINGHAM 1995          
CC                                                                      
C                                                                       
      WRITE(6,101)                                                      
  101 FORMAT(//'  ANGMOM - angular momentum calculations'//             
     > '         - by N M Clarke, Univ. of Birmingham 1995, vs 4.0'//)  
  100 WRITE(6,102)                                                      
  102 FORMAT('  Pick calculation from menu'//                           
     *       '  For                              Enter'//               
     *       '  CLEBSCH-GORDAN COEFFICIENT        1'/                   
     *       '  WIGNER 3-J SYMBOL                 2'/                   
     *       '  RACAH COEFFICIENT                 3'/                   
     *       '  WIGNER 6-J SYMBOL                 4'/                   
     *       '  WIGNER 9-J SYMBOL                 5'/                   
     *       '  ROTATION MATRIX                   6'/                   
     *       '  SPHERICAL HARMONIC                7'/                   
     *       '  REDUCED M.E. SPHERICAL TENSORS    8'/
     *       '  ASSOCIATED LEGENDRE FUNCTIONS     9'/                   
     *       '  EXIT                              0'//)                 
      READ(5,1)IA                                                       
    1 FORMAT(I1)                                                        
      IF(IA.LE.0.OR.IA.GT.9) GO TO 99                                   
      IF( IA.EQ.1 ) CALL CG                                             
      IF(IA.EQ.2 ) CALL W3J                                             
      IF(IA.EQ.3  )CALL RAC                                             
      IF(IA.EQ.4 )CALL W6J                                              
      IF(IA.EQ.5 ) CALL W9J                                             
      if (ia.eq.6)call romat                                            
      if(ia.eq.7) call fylm                                             
      if(ia.eq.8) call sphtens
      if(ia.eq.9) call alegen                                          
      GO TO 100                                                         
  99  STOP                                                              
      END                                                               
C                                                                       
C
      subroutine alegen
c
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C
C     ASSOCIATED LEGENDRE FUNCTIONS OF +VE OR -VE DEGREE
C     AND +VE OR -VE ORDER
      DATA PI/3.1415926535898D0/
      RAD=PI/1.8D2
 100   CONTINUE
      WRITE(6,*)' LEGENDRE FUNCTIONS P(L,M,Theta)'
      WRITE(6,*)' L and M may be  +VE or -VE, Theta in degrees'
  22  WRITE(6,*)' INPUT L'
      READ(5,*)L
  23  WRITE(6,*)' INPUT M'
      READ(5,*)M
  24  WRITE(6,*)' INPUT Thetastart,Thetastop,Thetastep (Degs.)'
      READ(5,*)THST,THSTOP,THSTEP
c
      write(6,98)
  98  format('  L   M     Theta(Deg.)    P(L,M,Theta)'
     >'       [P(L,M,Theta)]**2'/)
c
      THETA=THST
11    THRAD=THETA*RAD
      CALL LGNDR(L,M,THRAD,PLM)
      PLMSQ=PLM*PLM
       WRITE(6,101) L,M,THETA,PLM,PLMSQ
101   FORMAT(i3,1x,i3,2x,f10.4,6x,1pd14.6,6x,d14.6)
      THETA=THETA+THSTEP
      IF(THETA.LE.THSTOP) GO TO 11
      WRITE(6,*)' INPUT   1  TO EXIT'
      WRITE(6,*)'         2  FOR NEW DEGREE (L)'
      WRITE(6,*)'         3  FOR NEW ORDER  (M)'
      WRITE(6,*)'         4  FOR NEW ANGLES'
      READ(5,*)K
      GOTO(99,22,23,24),K
 99   return
      END
C
C     ASSOCIATED LEGENDRE FUNCTIONS PLM(THETA)        
C
C
      SUBROUTINE LGNDR(LL,MM,THETA,PLM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PLM=0.0D0
      MABS=IABS(MM)
      LABS=IABS(LL)
      IF(LL.LT.0)LABS=IABS(LL+1)
      IF(MABS.GT.LABS) RETURN
      MPLUS=MABS+1    
      LPLUS=LABS+1
      Y=DCOS(THETA)
      Z=DSIN(THETA) 
      IX=0
      DO 100 M=1,MPLUS
      LX=M-1
      L2=0
      P3=1.0
      FL1=LX
      IF(LX.EQ.0) GO TO 41
      DO 40 LT=1,LX
      FL1=FL1+1.0
      P3=P3*FL1*Z/2.0
   40 CONTINUE
   41 P2=0.0
      FL2=FL1+1.0
      FL3=1.0
      DO 90 LT=1,LPLUS
      IX1=IX+LT
      IF(L2-LX)50,70,60
   50 PLM     =0.0D0
      GO TO 75
   60 P3=(FL2*Y*P2-FL1*P1)/FL3
      FL1=FL1+1.0
      FL2=FL2+2.0
      FL3=FL3+1.0
   70 IF(ABS(P3).LT.1.0E-18) P3=0.0
      PLM     =P3
      P1=P2
      P2=P3
   75 L2=L2+1
   90 CONTINUE
      IX=IX+LPLUS
  100 CONTINUE
      IF(MM.LT.0) THEN
      PLM=PLM*(-1)**MABS* YXFCT(LABS+MABS,LABS-MABS)    
      ENDIF
      RETURN
      END
C
C
                                                                       
C                                                                       
c                                                                       
      subroutine sphtens                                                
c                                                                       
      implicit double precision(a-h,o-z)                                
c                                                                       
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
      data sqr4pi/3.544907702d0/                                        
c                                                                       
c     calculation of spherical tensor matrix elements                   
c      as Bell & Satchler Nuc. data Tables  A9 1971                     
c                                                                       
        write(6,1)                                                      
 1     format('**********************************************',         
     >'*********'//                                                     
     > ' Single particle reduced matrix'                                
     >' elements of spherical tensors'/                                 
     > '  W K Bell & G R Satchler, Nucl. Data Tables'                   
     >' 9, 2 (1971) 147'//                                              
     *'< lf,jf || T ( L, S, J ) ||  li,ji > = GSF* GLF'//)              
c                                                                       
 100  write(6,*)' INPUT lf,jf, L,S,J, li,ji'                            
      read(5,*),lm2,aj2,l,is,j,lm1,aj1                                  
c                                                                       
      ll=2*l                                                            
       ks=2*is                                                          
        jj=2*j                                                          
c                                                                       
         ll1=2*lm1                                                      
          j1=nint(2.0d0*aj1)                                            
           is1=iabs(j1-ll1)                                             
c                                                                       
           ll2=2*lm2                                                    
            j2=nint(2.0d0*aj2)                                          
             is2=iabs(j2-ll2)                                           
c                                                                       
c                                                                       
      CALL GS(LL,KS,JJ, 0,LL1,J1,IS1,LL2,J2,IS2,GSF)                    
      CALL GL(SQR4PI,LL,LL1,LL2,GLF)                                    
                                                                        
c                                                                       
      tlsj=gsf*glf                                                      
c                                                                       
       WRITE(6,90)lm2,aj2,l,is,j,lm1,aj1,gsf,glf,tlsj                   
  90  FORMAT(/                                                          
     >'<',i2,','f4.1,' || T (',i2,',',i2,',',i2,') || '                 
     > i2,',',f4.1,' >'/                                                
     > ' = ( GSF= ',1pd13.5,')*( GLF= ',d13.5,')'/                      
     > ' = <lf,jf ||T(L,S,J)|| li,ji> =',d13.5//)                       
c                                                                       
c                                                                       
      write(6,99)                                                       
  99  format('  Enter E or e  to EXIT; any other key to continue'/)     
      READ(5,4)IE                                                       
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      end                                                               
c                                                                       
c                                                                       
c                                                                       
      SUBROUTINE GS(LL,KS,JJ,KT,LL1,J1,IS1,LL2,J2,IS2,GSF)              
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      HJJ=REAL(JJ+1)                                                    
      HLL=REAL(LL+1)                                                    
      HSV=REAL(KS+1)                                                    
      HJ1=REAL(J1+1)                                                    
      GSF=DSQRT(2.0D0*HJJ*HLL*HSV*HJ1)*                                 
     * WINEJ(LL2,LL1,LL,IS2,IS1,KS,J2,J1,JJ)                            
      RETURN                                                            
      END                                                               
      SUBROUTINE GL(SQR4PI,LL,LL1,LL2,GLF)                              
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      GLF=SQRT(DBLE(REAL(LL1+1)))/SQR4PI*VCC(LL1,LL,LL2,0,0)*           
     *    PHASEF((LL+LL1-LL2)/4)                                        
      RETURN                                                            
      END                                                               
c                                                                       
c                                                                       
c                                                                       
c                                                                       
      SUBROUTINE CG                                                     
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
 100  WRITE(6,1)                                                        
    1 FORMAT('  INPUT J1,M1,J2,M2,J3 '/)                                
      READ(5,*)AJ1,AM1,AJ2,AM2,AJ3                                      
      J1=IDINT(AJ1+AJ1)                                                 
      J2=IDINT(AJ2+AJ2)                                                 
      J3=IDINT(AJ3+AJ3)                                                 
      M1=IDINT(AM1+AM1)                                                 
      M2=IDINT(AM2+AM2)                                                 
      AM3=AM2+AM1                                                       
      ANS= VCC(J1,J2,J3,M1,M2)                                          
      WRITE(6,3)                                                        
    3 FORMAT(/'  CLEBSCH-GORDAN COEFFICIENT'/)                          
      WRITE(6,2)AJ1,AM1,AJ2,AM2,AJ3,AM3,ANS                             
   2  FORMAT(' <',F4.1,',',F4.1,':',F4.1,',',F4.1,'|',                  
     *F4.1,',',F4.1,'>  = ',F15.9//                                     
     *       '  Enter E or e  to EXIT; any other key to continue'/)     
      READ(5,4)IE                                                       
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      END                                                               
C                                                                       
C                                                                       
C                                                                       
      SUBROUTINE W3J                                                    
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
 100  WRITE(6,1)                                                        
    1 FORMAT('  INPUT J1,J2,J3'/)                                       
      READ(5,*)AJ1,AJ2,AJ3                                              
      WRITE(6,11)                                                       
   11 FORMAT('  INPUT M1,M2'/)                                          
      READ(5,*) AM1,AM2                                                 
      J1=IDINT(AJ1+AJ1)                                                 
      J2=IDINT(AJ2+AJ2)                                                 
      J3=IDINT(AJ3+AJ3)                                                 
      M1=IDINT(AM1+AM1)                                                 
      M2=IDINT(AM2+AM2)                                                 
      AM3=-(AM2+AM1)                                                    
      M3=AM3+AM3                                                        
      ANS= VCC(J1,J2,J3,M1,M2)*PHASEF((J1-J2-M3)/2)/SQRT(AJ3+AJ3+1.0)   
      WRITE(6,3)                                                        
    3 FORMAT(/'  WIGNER 3-J SYMBOL'/)                                   
      WRITE(6,2)AJ1,AJ2,AJ3,AM1,AM2,AM3,ANS                             
   2  FORMAT(' (',F4.1,' ',F4.1,' ',F4.1,')'/                           
     *       ' (',F4.1,' ',F4.1,' ',F4.1,')   = ',F15.9//               
     *       '  Enter E or e  to EXIT; any other key to continue'/)     
      READ(5,4)IE                                                       
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      END                                                               
C                                                                       
C                                                                       
C                                                                       
      SUBROUTINE RAC                                                    
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
 100  WRITE(6,1)                                                        
    1 FORMAT('  INPUT J1,J2,J3'/)                                       
      READ(5,*)AJ1,AJ2,AJ3                                              
      WRITE(6,11)                                                       
   11 FORMAT('  INPUT J4,J5,J6'/)                                       
      READ(5,*) AJ4,AJ5,AJ6                                             
      J1=IDINT(AJ1+AJ1)                                                 
      J2=IDINT(AJ2+AJ2)                                                 
      J3=IDINT(AJ3+AJ3)                                                 
      J4=IDINT(AJ4+AJ4)                                                 
      J5=IDINT(AJ5+AJ5)                                                 
      J6=IDINT(AJ6+AJ6)                                                 
      ANS= RACAH(J1,J2,J3,J4,J5,J6)                                     
      WRITE(6,3)                                                        
    3 FORMAT(/'  RACAH COEFFICIENT'/)                                   
      WRITE(6,2)AJ1,AJ2,AJ3,AJ4,AJ5,AJ6,ANS                             
   2  FORMAT(' W(',F4.1,',',F4.1,',',F4.1,',',F4.1,'|',                 
     *  F4.1,',',F4.1,')  = ',F15.9//                                   
     *       '  Enter E or e  to EXIT; any other key to continue'/)     
       READ(5,4)IE                                                      
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      END                                                               
C                                                                       
C                                                                       
C                                                                       
      SUBROUTINE W6J                                                    
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
 100  WRITE(6,1)                                                        
    1 FORMAT('  INPUT J1,J2,J3'/)                                       
      READ(5,*)AJ1,AJ2,AJ3                                              
      WRITE(6,11)                                                       
   11 FORMAT('  INPUT J4,J5,J6'/)                                       
      READ(5,*) AJ4,AJ5,AJ6                                             
      J1=IDINT(AJ1+AJ1)                                                 
      J2=IDINT(AJ2+AJ2)                                                 
      J3=IDINT(AJ3+AJ3)                                                 
      J4=IDINT(AJ4+AJ4)                                                 
      J5=IDINT(AJ5+AJ5)                                                 
      J6=IDINT(AJ6+AJ6)                                                 
      ANS= RACAH(J1,J2,J3,J4,J5,J6)*PHASEF((J1+J2+J5+J4)/2)             
      WRITE(6,3)                                                        
    3 FORMAT(/'  WIGNER 6-J SYMBOL'/)                                   
      WRITE(6,2)AJ1,AJ2,AJ3,AJ4,AJ5,AJ6,ANS                             
   2  FORMAT( ' {',F4.1,' ',F4.1,' ',F4.1,'}'/                          
     *        ' {',F4.1,' ',F4.1,' ',F4.1,'}  = ',F15.9//               
     *       '  Enter E or e  to EXIT; any other key to continue'/)     
      READ(5,4)IE                                                       
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      END                                                               
C                                                                       
C                                                                       
C                                                                       
      SUBROUTINE W9J                                                    
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
 100  WRITE(6,1)                                                        
    1 FORMAT('  INPUT J1,J2,J3'/)                                       
      READ(5,*)AJ1,AJ2,AJ3                                              
      WRITE(6,11)                                                       
   11 FORMAT('  INPUT J4,J5,J6'/)                                       
      READ(5,*) AJ4,AJ5,AJ6                                             
      WRITE(6,12)                                                       
   12 FORMAT('  INPUT J7,J8,J9'/)                                       
      READ(5,*) AJ7,AJ8,AJ9                                             
      J1=IDINT(AJ1+AJ1)                                                 
      J2=IDINT(AJ2+AJ2)                                                 
      J3=IDINT(AJ3+AJ3)                                                 
      J4=IDINT(AJ4+AJ4)                                                 
      J5=IDINT(AJ5+AJ5)                                                 
      J6=IDINT(AJ6+AJ6)                                                 
      J7=IDINT(AJ7+AJ7)                                                 
      J8=IDINT(AJ8+AJ8)                                                 
      J9=IDINT(AJ9+AJ9)                                                 
      ANS= WINEJ(J1,J2,J3,J4,J5,J6,J7,J8,J9)                            
      WRITE(6,13)                                                       
   13 FORMAT(/'  WIGNER 9-J SYMBOL'/)                                   
      WRITE(6,2)AJ1,AJ2,AJ3,AJ4,AJ5,AJ6,AJ7,AJ8,AJ9,ANS                 
   2  FORMAT( ' {',F4.1,' ',F4.1,' ',F4.1,'}'/                          
     *        ' {',F4.1,' ',F4.1,' ',F4.1,'}'/                          
     *        ' {',F4.1,' ',F4.1,' ',F4.1,'}  = ',F15.9//               
     *       '  Enter E or e  to EXIT; any other key to continue'/)     
      READ(5,4)IE                                                       
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      END                                                               
C                                                                       
C                                                                       
C                                                                       
C     START OF ANGULAR MOMENTUM FUNCTIONS AND SUBROUTINES               
C                                                                       
C                                                                       
      DOUBLE PRECISION FUNCTION VCC(JX1,JX2,JX3,MX1,MX2)                
C                                                                       
C      CALCULATES THE CLEBSCH-GORDON COEFF  ( J1,M1,J2,M2:J3,M3)        
C                                                                       
C                                                                       
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION F(82),FACT(81)                                          
      COMMON /FACTRL/ F                                                 
      EQUIVALENCE (FACT(1),F(2))                                        
C                                                                       
      VCC=0.0D0                                                         
      J1=JX1                                                            
      J2=JX2                                                            
      J3=JX3                                                            
      M1=MX1                                                            
      M2=MX2                                                            
      IF(J1.LT.J2) GO TO 20                                             
      IF(J3.LT.J2) GO TO 30                                             
      ICNTR=0                                                           
      GO TO 40                                                          
   20 IF(J3.LT.J1) GO TO 30                                             
      ICNTR=-1                                                          
      IT=J1                                                             
      J1=J2                                                             
      J2=IT                                                             
      IT=M1                                                             
      M1=M2                                                             
      M2=IT                                                             
      GO TO 40                                                          
   30 ICNTR=1                                                           
      IT=J2                                                             
      J2=J3                                                             
      J3=IT                                                             
      M2=-M1-M2                                                         
   40 CONTINUE                                                          
      JZ1=(J1+J2-J3)/2                                                  
      IF(JZ1.LT.0) GO TO 150                                            
      JZ2=(J1+J3-J2)/2                                                  
      IF(JZ2.LT.0) GO TO 150                                            
      JZ3=(J2+J3-J1)/2                                                  
      IF(JZ3.LT.0) GO TO 150                                            
      IF(J1-IABS(M1).LT.0) GO TO 150                                    
      IF(J2-IABS(M2).LT.0) GO TO 150                                    
      IF(J3-IABS(M1+M2).LT.0) GO TO 150                                 
      JT1=(J1-J3+M2)/2                                                  
      JT2=(J2-J3-M1)/2                                                  
      NUMIN=MAX0(JT1,JT2,0)                                             
      JT3=(J1-M1)/2                                                     
      JT4=(J2+M2)/2                                                     
      NUMAX=MIN0(JT3,JT4,JZ1)                                           
      JT5=(J2-M2)/2                                                     
      IF(NUMAX.LT.NUMIN) GO TO 150                                      
      J4=J1/2                                                           
      J5=J3/2                                                           
      PHAS=PHASEF(NUMIN)                                                
      FCTOR=YXFCT(J4,(J1+M1)/2)+YXFCT(J4,JT3)+YXFCT((J1+J2+J3)/2+1,JZ2) 
     1  +YXFCT(J5,(J3+M1+M2)/2)+YXFCT(J5,(J3-M1-M2)/2)                  
     2  +FACT(JZ1)+FACT(JZ3)+FACT(JT4)+FACT(JT5)+LOG(REAL(J3+1))        
      FCTOR=0.5D0*FCTOR                                                 
      DO 100 NU=NUMIN,NUMAX                                             
      TERM=FCTOR    +(YXFCT(JT3-NU,J4)+YXFCT(NU-JT2,J5))                
     1  -(FACT(JT4-NU)+FACT(NU-JT1)+FACT(JZ1-NU)+FACT(NU))              
      VCC=VCC+PHAS*DEXP(TERM)                                           
      PHAS=-PHAS                                                        
  100 CONTINUE                                                          
      IF(ICNTR) 120,150,110                                             
  110 VCC=VCC*SQRT(REAL(J2+1)/REAL(J3+1))*PHASEF(JT3)                   
      GO TO 150                                                         
  120 VCC=VCC*PHASEF(JZ1)                                               
  150 RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
C     RACAH COEFFICIENTS                                                
C                                                                       
C                                                                       
      DOUBLE PRECISION FUNCTION RACAH(J1,J2,J3,J4,J5,J6)                
C                                                                       
C                                                                       
C     CALCULATES  W(J1,J2,J3,J4,J5,J6)                                  
C                                                                       
C                                                                       
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION F(82),FACT(81)                                          
      COMMON /FACTRL/ F                                                 
      EQUIVALENCE (FACT(1),F(2))                                        
C                                                                       
      RACAH=0.0D0                                                       
      Z1=DELR(J1,J2,J5)                                                 
      IF(Z1.LE.-999.0) GO TO 90                                         
      Z1=DELR(J3,J4,J5)+Z1                                              
      IF(Z1.LE.-999.0) GO TO 90                                         
      Z2=DELR(J1,J3,J6)                                                 
      IF(Z2.LE.-999.0) GO TO 90                                         
      Z2=DELR(J2,J4,J6)+Z2                                              
      IF(Z2.LE.-999.0) GO TO 90                                         
      Z1=0.5*(Z1+Z2)                                                    
      JT1=(J1+J2+J5)/2                                                  
      JT2=(J3+J4+J5)/2                                                  
      JT3=(J1+J3+J6)/2                                                  
      JT4=(J2+J4+J6)/2                                                  
      JZ1=(J1+J2+J3+J4)/2                                               
      JZ2=(J1+J4+J5+J6)/2                                               
      JZ3=(J2+J3+J5+J6)/2                                               
      NUMIN=MAX0(JT1,JT2,JT3,JT4)                                       
      NUMAX=MIN0(JZ1,JZ2,JZ3)                                           
      IF(NUMAX.LT.NUMIN) GO TO 90                                       
      PHASE=PHASEF(NUMIN+JZ1)                                           
      DO 80 NU=NUMIN,NUMAX                                              
      JY1=NU-JT1                                                        
      JY2=NU-JT2                                                        
      JY3=NU-JT3                                                        
      JY4=JZ1-NU                                                        
      JY5=JZ2-NU                                                        
      JY6=JZ3-NU                                                        
      FCTOR=FACT(JY1)+FACT(JY2)+FACT(JY3)+YXFCT(NU+1,NU-JT4)            
     1  +FACT(JY4)+FACT(JY5)+FACT(JY6)                                  
      ARG=Z1-FCTOR                                                      
      IF(ARG.GT.-40.0) RACAH=RACAH+PHASE*DEXP(ARG)                      
      PHASE=-PHASE                                                      
   80 CONTINUE                                                          
   90 RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
C                                                                       
C     WIGNER 9-J COEFFICIENTS                                           
C                                                                       
C                                                                       
      DOUBLE PRECISION FUNCTION WINEJ(J1,J2,J3,J4,J5,J6,J7,J8,J9)       
C                                                                       
C     CALCULATES {J1,J2,J3}                                             
C                {J4,J5,J6}                                             
C                {J7,J8,J9                                              
C                                                                       
C                                                                       
                                                                        
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      WINEJ=0.0                                                         
      NUMIN=MAX0(IABS(J1-J9),IABS(J2-J6),IABS(J4-J8))                   
      NUMAX=MIN0(J1+J9,J2+J6,J4+J8)                                     
      IF(NUMAX.LT.NUMIN) GO TO 40                                       
      DO 20 NUX=NUMIN,NUMAX,2                                           
      NU=NUX                                                            
      PROD=RACAH(J1,J4,J9,J8,J7,NU)*RACAH(J2,J5,NU,J4,J8,J6)            
     1    *RACAH(J9,NU,J3,J2,J1,J6)*REAL(NU+1)                          
      WINEJ=WINEJ+PROD                                                  
   20 CONTINUE                                                          
      WINEJ=WINEJ*PHASEF((J1+J3+J5+J8)/2+J2+J4+J9)                      
   40 RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
CC     PHASEF  GIVES  I**L FACTORS                                      
C                                                                       
C                                                                       
      DOUBLE PRECISION FUNCTION PHASEF(N)                               
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      PHASEF=REAL(1-2*IABS(N-2*(N/2)))                                  
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
C                                                                       
CC     YXFACT CALCULATES LOG FACTORIALS  N!/ M!                         
C                                                                       
C                                                                       
      DOUBLE PRECISION FUNCTION YXFCT(M,N)                              
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C     COMPUTES LN(NFACT/MFACT)                                          
C                                                                       
C     A TABLE OF LOG(N FACTORIAL), N = 0 TO 81.                         
C                                                                       
      DOUBLE PRECISION A,B                                              
      DIMENSION F(82),A(52),B(30)                                       
      COMMON /FACTRL/ F                                                 
      DATA A/0.0D0, 0.0D0                                               
     1    ,6.9314718055995D-01,1.7917594692281D+00,3.1780538303479D+00  
     2    ,4.7874917427820D+00,6.5792512120101D+00,8.5251613610654D+00  
     3    ,1.0604602902745D+01,1.2801827480081D+01,1.5104412573076D+01  
     4    ,1.7502307845874D+01,1.9987214495662D+01,2.2552163853123D+01  
     5    ,2.5191221182739D+01,2.7899271383841D+01,3.0671860106081D+01  
     6    ,3.3505073450137D+01,3.6395445208033D+01,3.9339884187199D+01  
     7    ,4.2335616460753D+01,4.5380138898477D+01,4.8471181351835D+01  
     8    ,5.1606675567764D+01,5.4784729398112D+01,5.8003605222981D+01  
     9    ,6.1261701761002D+01,6.4557538627006D+01,6.7889743137182D+01  
     *    ,7.1257038967168D+01,7.4658236348830D+01,7.8092223553315D+01  
     1    ,8.1557959456115D+01,8.5054467017582D+01,8.8580827542198D+01  
     2    ,9.2136175603687D+01,9.5719694542143D+01,9.9330612454787D+01  
     3    ,1.0296819861451D+02,1.0663176026064D+02,1.1032063971476D+02  
     4    ,1.1403421178146D+02,1.1777188139975D+02,1.2153308151544D+02  
     5    ,1.2531727114936D+02,1.2912393363913D+02,1.3295257503562D+02  
     6    ,1.3680272263733D+02,1.4067392364823D+02,1.4456574394634D+02  
     7    ,1.4847776695177D+02,1.5240959258450D+02/                     
      DATA B/1.5636083630308D+02                                        
     1    ,1.6033112821663D+02,1.6432011226320D+02,1.6832744544843D+02  
     2    ,1.7235279713916D+02,1.7639584840700D+02,1.8045629141754D+02  
     3    ,1.8453382886145D+02,1.8862817342367D+02,1.9273904728784D+02  
     4    ,1.9686618167289D+02,2.0100931639928D+02,2.0516819948264D+02  
     5    ,2.0934258675254D+02,2.1353224149456D+02,2.1773693411395D+02  
     6    ,2.2195644181913D+02,2.2619054832373D+02,2.3043904356578D+02  
     7    ,2.3470172344282D+02,2.3897838956183D+02,2.4326884900298D+02  
     8    ,2.4757291409619D+02,2.5189040220972D+02,2.5622113555001D+02  
     9    ,2.6056494097186D+02,2.6492164979855D+02,2.6929109765102D+02  
     *    ,2.7367312428569D+02,2.7806757344037D+02/                     
C                                                                       
      DATA K/0/                                                         
      IF (K .NE. 0) GO TO 19                                            
      DO 10 I=1,52                                                      
 10   F(I)=A(I)                                                         
      DO 11 I=1,30                                                      
 11   F(I+52)=B(I)                                                      
      K=1                                                               
 19   YXFCT=0.0d0                                                       
      NUMAX=M-N                                                         
      IF(NUMAX) 30,100,20                                               
   20 ICTRL=0                                                           
      L1=N                                                              
      L2=M                                                              
      GO TO 40                                                          
   30 ICTRL=1                                                           
      NUMAX=-NUMAX                                                      
      L1=M                                                              
      L2=N                                                              
   40 CONTINUE                                                          
      IF(L1.GT.81) GO TO 45                                             
      IF(L2.GT.81) GO TO 42                                             
      YXFCT=F(L2+1)-F(L1+1)                                             
      GO TO 90                                                          
   42 CONTINUE                                                          
      YXFCT=F(82)-F(L1+1)                                               
      L1=81                                                             
      NUMAX=L2-L1                                                       
   45 CONTINUE                                                          
      FCTOR=L1                                                          
      TERM=1.0d0                                                        
      DO 50 NU=1,NUMAX                                                  
      FCTOR=FCTOR+1.0d0                                                 
      TERM=TERM*FCTOR                                                   
      IF(TERM.LT.1.E20) GO TO 50                                        
      YXFCT=YXFCT+LOG(TERM)                                             
      TERM=1.0d0                                                        
   50 CONTINUE                                                          
      YXFCT=YXFCT+LOG(TERM)                                             
   90 CONTINUE                                                          
      IF(ICTRL.EQ.0) YXFCT=-YXFCT                                       
  100 RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
C                                                                       
CC    DELR CALCULATES THE TRIANGULAR FACTORS                            
C                                                                       
C                                                                       
      DOUBLE PRECISION FUNCTION DELR(J1,J2,J3)                          
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION F(82),FACT(81)                                          
      COMMON /FACTRL/ F                                                 
      EQUIVALENCE (FACT(1),F(2))                                        
C                                                                       
      JZ1=(J1+J2-J3)/2                                                  
      IF(JZ1.LT.0) GO TO 130                                            
      JZ2=(J1-J2+J3)/2                                                  
      IF(JZ2.LT.0) GO TO 130                                            
      JZ3=(J2+J3-J1)/2                                                  
      IF(JZ3.LT.0) GO TO 130                                            
      JZ4=(J1+J2+J3)/2+1                                                
      IF(JZ3.LT.JZ2) GO TO 80                                           
      IF(JZ3.LT.JZ1) GO TO 70                                           
      DELR=YXFCT(JZ4,JZ3)+FACT(JZ1)+FACT(JZ2)                           
      GO TO 150                                                         
   70 DELR=YXFCT(JZ4,JZ1)+FACT(JZ2)+FACT(JZ3)                           
      GO TO 150                                                         
   80 IF(JZ2.LT.JZ1) GO TO 70                                           
      DELR=YXFCT(JZ4,JZ2)+FACT(JZ1)+FACT(JZ3)                           
      GO TO 150                                                         
  130 DELR=-1000.                                                       
  150 RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
C                                                                       
C                                                                       
C                                                                       
      SUBROUTINE romat                                                  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
 100  WRITE(6,1)                                                        
    1 FORMAT('  INPUT J,M,N,ALPHA,BETA,GAMMA'
     >' - angles in degrees'/)                        
      READ(5,*)AJ,AM,AN,alpha,beta,gamma                                
       call djmn(aj,am,an,alpha,beta,gamma,d,dr,di)                     
      WRITE(6,3)                                                        
    3 FORMAT(/'  ROTATION MATRIX'/)                                     
      write(6,21)aj,am,an,d                                             
  21  format(' d( J=',f4.1,', M=',f5.1,', N=',f5.1,')'/                 
     > ' = ',1pe12.4/)                                                  
      WRITE(6,2)aj,am,an,alpha,beta,gamma,dr,di                         
   2  FORMAT(' D( J=',F4.1,', M=',F5.1,', N=',F5.1')'/                  
     >       '  (Alpha=',f7.2,', Beta=',f7.2,', Gamma=',f7.2,')'//      
     >       ' = ', 1pe15.7,' + i ',e15.7//                             
     *       '  Enter E or e  to EXIT; any other key to continue'/)     
      READ(5,4)IE                                                       
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      END                                                               
C                                                                       
C                                                                       
                                                                        
c                                                                       
      subroutine djmn(aj,am,an,alpha,beta,gamma,d,dr,di)                
      implicit double precision (a-h,o-z)                               
c                                                                       
       data rad/0.017453293d0/                                          
c                                                                       
c                                                                       
c    calculates the rotation matrix Djmn(alpha,beta,gamma)              
c                                                                       
      d=0.0d0                                                           
      den1= aj+am                                                       
      den2=aj-an                                                        
      den3= -(an-am)                                                    
c                                                                       
c      find maxixmum and minimum value of t                             
c                                                                       
       itmax=nint(min(den1,den2))                                       
       itmin=nint(max(den3,0.0))                                        
c        write(6,*) itmin,itmax                                         
cc                                                                      
c      get angle factors                                                
c                                                                       
      beta2=0.5d0*beta*rad                                              
       amcn= (alpha*am +gamma*an)*rad                                   
      cb=cos(beta2)                                                     
      sb=sin(beta2)                                                     
c      write(6,4) cb,sb                                                 
c  4   format('cb=',1pe12.5,' sb=',e12.5)                               
c                                                                       
c                                                                       
c     get numerator for djmn(beta)                                      
c                                                                       
       fnum1= yxfct(1,nint(aj+am))                                      
       fnum2=yxfct(1,nint(aj-am))                                       
       fnum3=yxfct(1,nint(aj+an))                                       
       fnum4=yxfct(1,nint(aj-an))                                       
c       write(6,5)fnum1,fnum2,fnum3,fnum4                               
c  5     format( 1pe12.5,3(1x,e12.5))                                   
      fnum=exp(0.5d0*(fnum1+fnum2+fnum3+fnum4))                         
c      write(6,*)fnum                                                   
                                                                        
                                                                        
c                                                                       
c      sum over t                                                       
c                                                                       
      sum=0.0d0                                                         
       do 100 i=itmin,itmax                                             
      fden1=yxfct(nint(aj+am-i),1)                                      
      fden2=yxfct(nint(aj-an-i),1)                                      
      fden3=yxfct(i,1)                                                  
      fden4=yxfct(nint(i+an-am),1)                                      
      term= (-1)**i * exp(fden1+fden2+fden3+fden4)                      
      term=term* (cb)**(2*aj+am-an-2*i) * (sb)**(2*i+an-am)             
c       write(6,3) i,term                                               
c  3   format('i=',i3,' term=',1pe12.5)                                 
      sum=sum+term                                                      
 100  continue                                                          
       d=sum*fnum                                                       
        dr=d*cos(amcn)                                                  
        di= -d*sin(amcn)                                                
c                                                                       
c                                                                       
      return                                                            
      end                                                               
c                                                                       
c                                                                       
c                                                                       
      subroutine fylm                                                   
      implicit double precision(a-h,o-z)                                
c                                                                       
      LOGICAL*1 IE,EE,EEE                                               
      DATA EE/'E'/,EEE/'e'/                                             
c                                                                       
 100  WRITE(6,1)                                                        
    1 FORMAT('  INPUT L (or k), M (or q), THETA, PHI '/)                
      READ(5,*)alk,amq,theta,phi                                        
c                                                                       
c                                                                       
       call ylm(alk,amq,theta,phi,ylmr,ylmi)                            
                                                                        
      WRITE(6,3)                                                        
    3 FORMAT(/'  SPHERICAL HARMONIC'/)                                  
      write(6,21)alk,amq,theta,phi,ylmr,ylmi                            
  21  format(' Y(L=',f4.1,', M=',f5.1,', Theta=',f7.2,                  
     > ', Phi=',f7.2,')'//                                              
     >       ' = ', 1pe15.7,' + i ',e15.7//                             
     *       '  Enter E or e  to EXIT; any other key to continue'/)     
      READ(5,4)IE                                                       
    4 FORMAT(A1)                                                        
      IF(IE.EQ.EE.OR.IE.EQ.EEE)RETURN                                   
      GOTO 100                                                          
      END                                                               
                                                                        
c                                                                       
c                                                                       
      subroutine ylm( alk,amq,theta,phi,ylmr,ylmi)                      
      implicit double precision(a-h,o-z)                                
c                                                                       
       data pi/3.141592654d0/                                           
                                                                        
       sign=+1.0d0                                                      
       amqq=abs(amq)                                                    
c                                                                       
        alpha=  -phi                                                    
         beta=theta                                                     
       call djmn(alk,amqq,0.0,alpha,beta,0.0,d,dr,di)                   
c                                                                       
c                                                                       
c                                                                       
       fac= sqrt(( 2.0d0*alk +1.0)/(4.0d0*pi))                          
       ylmr=fac*dr                                                      
       ylmi=fac*di                                                      
      if(amq.lt.0.0) then                                               
       ylmr= (-1)**amqq *ylmr                                           
        ylmi= (-1)**amqq * (-ylmi)                                      
      endif                                                             
c                                                                       
      return                                                            
      end                                                               
c                                                                       
c                                                                       
                                                                        
                                                                        
