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