***** ****************************************************************** * Name: SHY, Spherical Harmonics Y * Author: Norbert E. Ligterink, University of Pittsburgh * Date: 5 April 2003 (revision 19 March 2004) * Version: 1.1 (after comments from Maxim Tafipolsky) * Usage: use as you like, no warranty, acknowledgments appreciated ***** ****************************************************************** * * Using the formula 3.153 of Biederharn and Louck (Angular Momentum * in Quantum Physics) p.71, it generates solid spherical harmonics * SY (which are related to ordinary spherical harmonics Y through: * SY(l,m) = r^l Y(l,m) where r = SQRT(X**2+Y**2+Z**2) ) * which are used to calculate matrix elements of the type: * * * * normalized, so: = delta_LL' delta_MM' * * The routines stand alone, generate high-precision answers, and * are pretty fast. However, accuracy and speed can be improved * slightly by implementing the Gamma function, and maybe tabulating * its values. (Here it is not done so for portability.) * At the moment deviations are in the 15th-16th decimal place. * ***** ****************************************************************** * * To generate solid spherical harmonics in Cartesian coordinates: * SUBROUTINE Y(L,M,COEFF,LENGTH,INDEX,T,N) (see there) * * In: L, M * Out: COEFF(1:LENGTH2,1:LENGTH = * Int d Omega Y^*_LM x^X(1) y^X(2) z^X(3) Y_L'M' * * ***** ****************************************************************** * Just a little program to test. Please, use the subroutines in * any way you like. * PROGRAM fshy ***** ****************************************************************** IMPLICIT NONE INTEGER T, N, L, M, LL, MM, X(3) PARAMETER ( T = 3, N = 256) COMPLEX*16 MATEL, DM * 111 CONTINUE * PRINT *,"Give L,M (negative L will exit)" READ *,L,M IF(L .LT. 0)THEN GOTO 666 ELSE PRINT *,"Y(L=",L," M=",M,") =" CALL PRNTY(L,M,T,N) GOTO 111 ENDIF * 666 CONTINUE * PRINT *,"Give (negative L will exit)" READ *,L,M,X(1),X(2),X(3),LL,MM IF(L .LT. 0)THEN GOTO 999 ELSE DM = MATEL(L,M,LL,MM,X,T) PRINT *,"<",L,M,"|x^",X(1),"y^",X(2),"z^", c X(3),"|",LL,MM,">=",DM GOTO 666 ENDIF 999 CONTINUE * END ***** ****************************************************************** * The subroutine Y returns, given L,M, the solid spherical harmonic * r^l Y_lm in the Biederharn and Louck definition of p.71 * * LENGTH * Y = SUM COEFF(j) x^INDEX(1,j) y^INDEX(2,j) y^INDEX(3,j) * j=1 * SUBROUTINE Y(L,MM,COEFF,LENGTH,INDEX,T,N) IMPLICIT NONE INTEGER T,N INTEGER L,M,MM,LENGTH,INDEX(T,N), CM INTEGER I,J,K,IND,MAX,PHASE REAL*8 GAMMA(3*L+1),BINO(3*L+1,3*L+1),PI,PREFAC COMPLEX*16 COEFF(N) * * * because Y takes ABS of M, we use a dummy variable. * M = MM * PI = 4.0D0*ATAN(1.0D0) MAX = L+ABS(M)+1 * * Gamma function * GAMMA(1) = 1.0D0 DO I=1, MAX GAMMA(I+1) = DBLE(I)*GAMMA(I) ENDDO * * Binomial coefficients * DO I=0, MAX DO J=0, I BINO(I+1,J+1) = (GAMMA(I+1)/GAMMA(I-J+1))/GAMMA(J+1) ENDDO ENDDO * * normalization prefactor * PREFAC = SQRT(DBLE(2*L+1)*GAMMA(L+M+1)*GAMMA(L-M+1)/(4.0D0*PI)) * * negative values of M affect only the phase, controlled by CM * the phase is (i)^PHASE * * note: Y(lm) = (-1)^m Y*(l-m) (BL definition) * IF(M .LT. 0)THEN M = -M CM = 1 ELSE CM = 0 ENDIF * * IND = 1 * DO K=0, (L-M)/2 DO I=0, K DO J=0, M IF(CM .EQ. 1)THEN PHASE = MOD(4*L-2*K-J,4) ELSE PHASE = MOD(2*K+2*M+J,4) ENDIF COEFF(IND) = ((0.0D0,1.0D0)**PHASE)* c COMPLEX(PREFAC*BINO(M+1,J+1)*BINO(K+1,I+1)/ c ((2**(2*K+M))*(GAMMA(K+M+1)*GAMMA(K+1)*GAMMA(L-M-2*K+1))), c 0.0D0) INDEX(1,IND) = 2*I+M-J INDEX(2,IND) = 2*K-2*I+J INDEX(3,IND) = L-2*K-M IND = IND + 1 IF(IND .GT. N)THEN PRINT *,"Error, array overflow, increase N" STOP ENDIF ENDDO ENDDO ENDDO * LENGTH = IND - 1 * END ***** ****************************************************************** * Just a simple print routine to STDIO to generate functions to * paste in other programs (C-style power this case) * * calls Y() * SUBROUTINE PRNTY(L,M,T,N) IMPLICIT NONE INTEGER I,T,N INTEGER L,M,LENGTH,INDEX(T,N) COMPLEX*16 COEFF(N) * * CALL Y(L,M,COEFF,LENGTH,INDEX,T,N) * DO I=1, LENGTH-1 PRINT *,COEFF(I),"*x^(",INDEX(1,I),")*y^(",INDEX(2,I), c ")*z^(",INDEX(3,I),") +" ENDDO * I = LENGTH PRINT *,COEFF(I),"*x^(",INDEX(1,I),")*y^(",INDEX(2,I), c ")*z^(",INDEX(3,I),")" * END ***** ****************************************************************** * MATEL calculates the matrix element : * * = MATEL * * T = 3 (dimension of X) * * calls Y() and ANGULARINT() * COMPLEX*16 FUNCTION MATEL(LA,MA,LB,MB,X,T) IMPLICIT NONE INTEGER LA,MA,LB,MB,T,N PARAMETER ( N = 256 ) INTEGER PRFA,LENA,IXA(T,N) INTEGER PRFB,LENB,IXB(T,N) INTEGER A,B,X(T),CT,ST,CP,SP REAL*8 ANGULARINT COMPLEX*16 COEFA(N), COEFB(N), SUM * CALL Y(LA,MA,COEFA,LENA,IXA,T,N) CALL Y(LB,MB,COEFB,LENB,IXB,T,N) * SUM = COMPLEX(0.0D0,0.0D0) * * The angular integral for each term is: theta=[0:pi] phi=[0:2 pi] * cos^ct theta sin^st theta cos^cp phi sin^sp phi * DO A=1, LENA DO B=1, LENB SP = X(2) + IXB(2,B) + IXA(2,A) * * making sure sin^sp phi is even on the domain * IF(MOD(SP,2) .EQ. 0)THEN ST = X(1) + X(2) + IXA(1,A)+IXB(1,B) + IXA(2,A)+IXB(2,B) + 1 CT = X(3) + IXA(3,A) + IXB(3,B) CP = X(1) + IXB(1,B) + IXA(1,A) SUM = SUM + 2.0D0*COEFB(B)*CONJG(COEFA(A))* c COMPLEX(ANGULARINT(CT,ST)*ANGULARINT(CP,SP),0.D0) ENDIF ENDDO ENDDO * MATEL = SUM * END ***** ****************************************************************** * * returns the integral int_0^pi dx cos^N(x) sin^M(x) * REAL*8 FUNCTION ANGULARINT(N,M) IMPLICIT NONE INTEGER N,M,I REAL*8 GAMMA(N+M+4), PI * PI = 4.0D0*ATAN(1.0D0) * * GAMMA(i) is actually Gamma[i/2] * GAMMA(2) = 1.0D0 GAMMA(1) = SQRT(PI) * DO I=1, (N+M+2)/2 GAMMA(2*I+2) = DBLE(I)*GAMMA(2*I) GAMMA(2*I+1) = 0.5D0*DBLE(2*I-1)*GAMMA(2*I-1) ENDDO * IF(M .EQ. 0)THEN IF(MOD(N,2) .EQ. 0)THEN ANGULARINT = 2.0D0*GAMMA(1)*GAMMA(N+3)/(DBLE(N+1)*GAMMA(2+N)) ELSE ANGULARINT = 0.0D0 ENDIF ELSE IF(MOD(N,2) .EQ. 0)THEN ANGULARINT = GAMMA(M+1)*GAMMA(N+1)/GAMMA(N+M+2) ELSE ANGULARINT = 0.0D0 ENDIF ENDIF * END ***** ******************************************************************