1 double precision function tripleYlmint(l1,l2,l3,m1,m2) 2! evaluates \int d\Omega Y_{l1}^{m1} Y_{l2}^{m2} Y_{l3}^{m3} 3! where m3=-m1-m2 4! Requires functions lnfactorial and lngamma 5 implicit none 6 integer l1,l2,l3,m1,m2 7 double precision wig1,wig2,xterm,pi 8 parameter (pi=3.1415926535897932384626433832795028841972d0) 9 double precision wigner3j 10 external wigner3j 11 tripleYlmint=0.d0 12 if (l3.lt.abs(l1-l2).or.l3.gt.l1+l2) return 13 if (m1.lt.-l1.or.m1.gt.l1) return 14 if (m2.lt.-l2.or.m2.gt.l2) return 15 wig1=wigner3j(l1,l2,l3,0,0) 16 wig2=wigner3j(l1,l2,l3,m1,m2) 17 xterm=sqrt(dble((2*l1+1)*(2*l2+1)*(2*l3+1))/(4.d0*pi)) 18 tripleYlmint=xterm*wig1*wig2 19 return 20 end 21 22!*********************************************************************** 23 24 double precision function wigner3j(j1,j2,j3,m1,m2) 25! evaluates | j1 j2 j3 | 26! | m1 m2 m3 | 27! where m3=-m1-m2 28 implicit none 29 integer j1,j2,j3,m1,m2,m3 30 integer it,itmax,itmin 31 double precision xterm1,xterm2,sterm 32 double precision triangle 33 double precision lnfactorial 34 external lnfactorial 35 m3=-m1-m2 36 wigner3j=0.d0 37 if (j3.lt.abs(j1-j2).or.j3.gt.j1+j2) return 38 if (m1.lt.-j1.or.m1.gt.j1) return 39 if (m2.lt.-j2.or.m2.gt.j2) return 40 if (m3.lt.-j3.or.m3.gt.j3) return 41 triangle=(lnfactorial(j1+j2-j3)+lnfactorial(j1-j2+j3) & 42& +lnfactorial(-j1+j2+j3)-lnfactorial(j1+j2+j3+1))/2.d0 43 xterm1=(lnfactorial(j1+m1)+lnfactorial(j1-m1) & 44& +lnfactorial(j2+m2)+lnfactorial(j2-m2) & 45& +lnfactorial(j3+m3)+lnfactorial(j3-m3))/2.d0 46 xterm2=0.d0 47 itmin=max(0,j2-j3-m1,j1-j3+m2) 48 itmax=min(j1+j2-j3,j1-m1,j2+m2) 49 do it=itmin,itmax 50 sterm=exp(triangle+xterm1-lnfactorial(it) & 51& -lnfactorial(j3-j2+it+m1)-lnfactorial(j3-j1+it-m2) & 52& -lnfactorial(j1+j2-j3-it)-lnfactorial(j1-it-m1) & 53& -lnfactorial(j2-it+m2)) 54 xterm2=xterm2+((-1)**it)*sterm 55 enddo 56 wigner3j=((-1)**(j1-j2-m3))*xterm2 57 return 58 end 59 60!*********************************************************************** 61 62 function spharm(ll,mm,rr) 63! Spherical harmonic Y_ll^mm for radius vector rr 64! cosine of polar angle theta = rr(3)/|rr| 65! azimuthal angle phi = atan2(rr(2),rr(1)) 66 implicit none 67 double complex :: spharm 68 double precision :: rr(3) 69 integer :: ll,mm 70 double precision :: costheta,cosphi,sinphi,rmag 71 double complex :: expphi 72 integer :: im,mmm 73 double precision :: fact,prefact,pi,lpolyn 74 parameter (pi=3.1415926535897932384626433832795028841972d0) 75 external lpolyn 76 77 mmm=abs(mm) 78 if (ll.lt.mmm) then 79 write(6,*) 'bad input parameters in function spharm' 80 stop 81 endif 82 fact=1.d0 83 do im=ll-mmm+1,ll+mmm 84 fact=fact*im 85 enddo 86 costheta=rr(3)/sqrt(dot_product(rr,rr)) 87 rmag=sqrt(dot_product(rr(1:2),rr(1:2))) 88 if (rmag.le.0.d0) then 89 if (mmm.eq.0) then 90 expphi=(1.d0,0.d0) 91 else 92 expphi=(0.d0,0.d0) 93 endif 94 else 95 cosphi=rr(1)/rmag 96 sinphi=rr(2)/rmag 97 expphi=cosphi+(0.d0,1.d0)*sinphi 98 endif 99 prefact=sqrt((2.d0*ll+1.d0)/(fact*4.d0*pi)) 100 spharm=prefact*lpolyn(ll,mmm,costheta)*expphi**mmm 101 if (mm.lt.0) spharm=conjg(spharm)*(-1)**mmm 102 103 return 104 end 105 106!*********************************************************************** 107 108 function lpolyn(ll,mm,xx) 109! Legendre Polynomial P_ll^mm (xx) 110 implicit none 111 double precision lpolyn,xx 112 integer ll,mm 113 double precision term,fact,lpstart,lpl,lplm1,lplm2 114 integer im,il 115 116 if (ll.lt.mm.or.abs(xx).gt.1.d0) then 117 write(6,*) 'bad input parameters for function lpolyn' 118 stop 119 endif 120 term=sqrt(1.d0-xx*xx) 121 lpstart=1.d0 122 do im=1,mm 123 lpstart=-lpstart*term*(2*im-1) 124 enddo 125 if (ll.eq.mm) then 126 lpolyn=lpstart 127 else 128 lplm1=lpstart 129 lplm2=0.d0 130 do il=mm+1,ll 131 lpl=(xx*(2*il-1)*lplm1-(il+mm-1)*lplm2)/dble(il-mm) 132 lplm2=lplm1 133 lplm1=lpl 134 enddo 135 lpolyn=lpl 136 endif 137 138 return 139 end 140 141!*********************************************************************** 142 143 double precision function lngamma(xx) 144 implicit double precision (a-h,o-z) 145 double precision xx,x, cc(6), gam, sumg, work, sq2pi 146 integer i 147 data cc/76.18009172947146d0,-86.50532032941677d0, & 148& 24.01409824083091d0,-1.231739572450155d0, & 149& 0.1208650973866179d-2,-0.5395239384953d-5/ 150 x=xx-1.d0 151 sq2pi=2.5066282746310005d0 152 sumg=1.000000000190015d0 153 do i=1,6 154 sumg=sumg+cc(i)/(x+i) 155 enddo 156 work=x+5.5 157 lngamma=(x+0.5)*log(work)-work+log(sq2pi*sumg) 158 return 159 end 160 161!*********************************************************************** 162 163 double precision function lnfactorial(ii) 164 implicit none 165 integer ii 166 double precision factstore(0:100),lngamma 167 external lngamma 168 save factstore 169 data factstore/101*-1.d0/ 170 if (ii.le.100) then 171 if (factstore(ii).lt.0.d0) factstore(ii)=lngamma(dble(ii+1)) 172 lnfactorial=factstore(ii) 173 else 174 lnfactorial=lngamma(dble(ii+1)) 175 endif 176 return 177 end 178 179!*********************************************************************** 180 181 function sphbesselj(xx,nn) 182 implicit none 183 double precision sphbesselj,xx 184 integer nn 185 integer ii,np,sigfig 186 double precision tiny,dj0,dj1,djnn,djn,djnp1,djnm1,scale 187 188 if (nn.lt.0) then 189 write(6,*) 'bad argument for order n of j_n(xx) in sphbesselj' 190 write(6,*) 'nn = ',nn 191 stop 192 endif 193 sigfig=12 194 tiny=1.d-10 195 if (xx.lt.tiny) then 196 dj0=1.d0-xx**2/6.d0 197 else 198 dj0=sin(xx)/xx 199 endif 200 if (nn.eq.0) then 201 sphbesselj=dj0 202 return 203 endif 204 if (xx.lt.tiny) then 205 sphbesselj=1.d0 206 do ii=1,nn 207 sphbesselj=sphbesselj*xx/dble(2*ii+1) 208 enddo 209 elseif (xx.lt.2*nn) then 210 djn=tiny 211 djnp1=0.d0 212 np=int(sqrt(dble(nn*sigfig**2)))+1 213 do ii=nn+np,1,-1 214 djnm1=(2*ii+1)*djn/xx-djnp1 215 if (ii.eq.nn) djnn=djn 216 djnp1=djn 217 djn=djnm1 218 enddo 219 scale=djn/dj0 220 sphbesselj=djnn/scale 221 else 222 if (xx.lt.tiny) then 223 dj1=xx/3.d0-4.d0*xx**3/30.d0 224 else 225 dj1=sin(xx)/(xx**2)-cos(xx)/xx 226 endif 227 djnm1=dj0 228 djn=dj1 229 do ii=1,nn-1 230 djnp1=(2*ii+1)*djn/xx-djnm1 231 djnm1=djn 232 djn=djnp1 233 enddo 234 sphbesselj=djn 235 endif 236 237 return 238 end 239