C C Fortran code gmm01TrA.f for calculating scattering of electromagnetic C radiation by a randomly oriented aggregate of homogeneous spheres C (modified from gmm01f.f) C Yu-lin Xu C 7/2002: released to public C C----------------------------------------------------------------------- C C For questions/comments/suggestions/bugs/problems please contact C Yu-lin Xu at xu@astro.ufl.edu or yxu@erc.ufl.edu C C----------------------------------------------------------------------- C C October 2002: Thanks Piotr Flatau for debugging. C C----------------------------------------------------------------------- PROGRAM gmm01TrA implicit double precision (a-h,o-z) include 'gmm01f.par' C----------------------------------------------------------------------- C Two parameters in "gmm01f.par": nLp, np C nLp --- the maximum number of spheres, must be equal to or greater C than the sphere-number in an aggregate actually calculated C np --- the maximum scattering order in the incident and scattered C field expansions, must be equal to or greater than the C highest scattering order required in actual calculations; C usually, this can be estimated by Wiscombe's criterion C regarding the largest individual sphere in an aggregate C an example of "gmm01f.par": parameter (nLp=6,np=10) C C C ****** This code is more memory-demanding than the single-orientation C code gmm01f.f. The computer memory required is at the level of C (nLp*np^2)^2. C This code is written in double precision arithmetic. An individual C sphere can have a size parameter way beyond ~200. There is no limit C to the overall dimension of an aggregate. The largest individual C sphere size and the maximum number of spheres that an aggregate can C have depends on the availability of computer memory. ****** C C----------------------------------------------------------------------- C C This code requires two input files: gmmTr.in and a user-specified C datafile for aggregate configuration C C ****** An example for the input and output files ****** C C (gmmTr.in) C gold2.k C 0 C 1.d-12 200 C -1 C 0 90. C C (gold2.k) C 389.3652205 C 2 C 0. 0. 0. 6.2 0.52362125 1.63011385 C 0. 0. 14.8 6.2 0.52362125 1.63011385 C C (gmm01TrA.out) C gmm01TrA.out input file: gold2.k xv: 0.126 xs: 0.141 C C C 0.17001E+03 0.16940E+03 0.61371E+00 0.87511E+00 0.17001E+03 0.39623E-02 C C 0.88687E+00 0.88367E+00 0.32014E-02 0.45650E-02 0.12685E-04 0.39623E-02 C C 0.70391E+00 0.70137E+00 0.25410E-02 0.36233E-02 0.10068E-04 0.39623E-02 C C s.a. C 0.0 0.37795E-04 0.0000 0.18527E-04 0.37062E-06 0.37062E-06 0.18527E-04 C 90.0 0.19077E-04 0.9412 0.18330E-04 0.37284E-06 0.18654E-06 0.18769E-06 C 180.0 0.37012E-04 0.0000 0.18134E-04 0.37157E-06 0.37157E-06 0.18134E-04 C C Scattering matrix (4X4 for each scattering angle): C 0.0 0.1889741E-04 0.0000000E+00 0.3855997E-20 0.1214528E-19 C 0.0000000E+00 0.1815616E-04 0.3855997E-20 0.1214528E-19 C -0.9639992E-21 -0.2494377E-19 0.1815616E-04 -0.1623600E-21 C 0.3036320E-20 0.6408125E-22 0.4531011E-21 0.1889741E-04 C 90.0 0.9538570E-05 -0.8978043E-05 0.1410391E-20 0.3185720E-20 C -0.9164338E-05 0.8979190E-05 0.5112782E-19 0.8300842E-20 C 0.1199542E-19 -0.1247253E-19 -0.1510833E-07 0.1523647E-07 C 0.3839321E-22 -0.1357579E-20 -0.1524935E-07 -0.1335003E-07 C 180.0 0.1850583E-04 0.0000000E+00 0.1005018E-18 0.1096689E-19 C 0.0000000E+00 0.1776269E-04 0.1005018E-18 0.1096689E-19 C 0.2512546E-19 0.6226571E-22 -0.1776269E-04 0.1530543E-21 C -0.2741722E-20 -0.2450208E-20 -0.4427613E-21 -0.1850583E-04 C C----------------------------------------------------------------------- C parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2,np2=2*np) parameter (NXMAX=3000,nangmax=1801) C----------------------------------------------------------------------- C NXMAX - the maximum dimension of an auxiliary array in calculating C Mie scattering coefficients, must not be less than C 1.1*np*|m|+10 C where |m| is the largest refractive index of all spheres) C----------------------------------------------------------------------- parameter (ni0=np*(np+1)*(np2+1)/3+np*np) parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6) parameter (nrc=4*np*(np+1)*(np+2)/3+np) parameter (nij=nLp*(nLp-1)/2) integer u,v,u0,nmax(nLp),uvmax(nLp),ind(nLp) double precision k,r0(6,nLp),x(nLp),dang(nangmax),r00(3,nLp), + rsr0(NXMAX),rsi0(NXMAX),rsx0(NXMAX),px0(NXMAX),w1(np),w2(np), + w3(np),w4(np),mue(4,4,nangmax),besj(0:np2+1),besy(0:np2+1), + i11(nangmax),i21(nangmax),i22(nangmax),i12(nangmax), + inat(nangmax),pol(nangmax),cexti(nLp),cabsi(nLp),cscai(nLp), + assymi(nLp),cpri(nLp),drot(nrc,nij),c0i(nLp),c1i(nLp), + bes0(nij),confg(5,nij),taup(2),taupj(2),taupg(2),taupjg(2), + tau0p(2),tau1p(2),tau2p(2),tau0pg(2),tau1pg(2),tau2pg(2), + tau0pj(2),tau1pj(2),tau2pj(2),tau0pjg(2),tau1pjg(2), + tau2pjg(2),tau20(2,2),tau11(2,2),tau02(2,2),tau20g(2,2), + tau11g(2,2),tau02g(2,2),taum(2,2),taumg(2,2),w01s(np2+1), + wmf1(np,np,0:np,0:np2),wm1(np,np,0:np,0:np2),wcf(np2+1), + wsdt(np,np,-np:np,0:np2) complex*16 A,B,cmz,Aj,Bj,A2,B2,Aj2,Bj2,A0,B0,ephi,ci,cin,ci0, + A1,B1,Aj1,Bj1,Aj0,Bj0,cmzj,cmzg,Ag,Bg, + A0p(2,2,2,2),A1p(2,2,2,2),B0p(2,2,2,2),B1p(2,2,2,2), + A0pg(2,2,2,2),A1pg(2,2,2,2),B0pg(2,2,2,2),B1pg(2,2,2,2), + atr0(ni0,nij),btr0(ni0,nij),atr(2,np,nmp),at(nmp),bt(nmp), + atr1(ni0,nij),btr1(ni0,nij),ek(np,nij),ref(nLp),ref0(nLp), + p0(nLp,nmp),q0(nLp,nmp),an(np),bn(np),aMie(nLp,np), + bMie(nLp,np),B2i(nLp),at0(nLp,nmp),bt0(nLp,nmp),at1(nmp), + bt1(nmp),as(nLp,nmp),bs(nLp,nmp),as2(nLp,nmp),bs2(nLp,nmp), + tta(nLp,nLp,nmp),ttb(nLp,nLp,nmp),tta0(nLp,nLp,nmp), + ttb0(nLp,nLp,nmp),asr(nLp,nLp,nmp),bsr(nLp,nLp,nmp), + as0(nLp,nLp,nmp),bs0(nLp,nLp,nmp),asc(nLp,nLp,nmp), + bsc(nLp,nLp,nmp),as1(nLp,nLp,nmp),bs1(nLp,nLp,nmp), + ast(nLp,nLp,nmp),bst(nLp,nLp,nmp),asp(nLp,nLp,nmp), + bsp(nLp,nLp,nmp),asv(nLp,nLp,nmp),bsv(nLp,nLp,nmp), + atj(nmp),btj(nmp),pct(nLp,nLp,2,2,nmp,nmp), + A1m(0:np,np,np,2,2),A2m(0:np,np,np,2,2), + B11n(np,np,2,2),B12n(np,np,2,2),B13n(np,np,2,2), + B21n(np,np,2,2),B22n(np,np,2,2),B23n(np,np,2,2), + fhmf1(0:np2,np,0:np,2),fmf1(0:np2,np,0:np,2), + fhm1v(0:np2,np,0:np,2),fm1v(0:np2,np,0:np,2), + fhm1q(0:np2,np,0:np,2),fm1q(0:np2,np,0:np,2), + fhmf1vq(0:np2,np,0:np,2),fmf1vq(0:np2,np,0:np,2), + fhas(0:np2,2,2),fnhs(0:np2,2,2) CHARACTER FLNAME*20,fileout*20,fileout1*19, + tailn*3,cnr2*2,cnr3*3,cnr1*1,flout*22 COMMON/MIESUB/ twopi,pih common/rot/bcof(0:np+2),dc(-np:np,0:nmp) common/fnr/fnr(0:4*(np+1)) common/pitau/pi(nmp0),tau(nmp0) common/tran/atr common/ig0/iga0(ni0) common/g0/ga0(ng0) common/cofmnv0/cof0(ni0) common/crot/cofsr(nmp) pih = dacos(0.d0) twopi = 4.d0*pih pione = 2.d0*pih ci=dcmplx(0.d0,1.d0) cin=dcmplx(0.d0,-1.d0) ci0=dcmplx(0.d0,0.d0) gcs=0.d0 gcv=0.d0 eps=1.d-20 fint=0.01d0 C------------------------------------------------------------------------ C fint: the interaction index in the range of [0,1] (normally fint=0.01) C In the scattering calculations, a quantity "f" for a pair of component C spheres is defined by f=(r_i+r_j)/d_{ij}, where (r_i,r_j) are the radii C of the two spheres and d_{ij} is the separation distance of the two C sphere-centers. When f', nL write(6,*) 'Change nLp in gmm01f.par, recompile, then try again' stop endif if(nL.eq.1) then idMie=1 endif C------------------------------------------------------------------------ C input the configuration and particle parameters for the aggregate C each line includes 6 numbers: C x-, y-, z-coordinates of the sphere-center, the radius of the sphere C in the same unit of the incident wavelength, the real and imaginary C parts of the refractive index C------------------------------------------------------------------------ do 1 i=1,nL read(2,*,err=10) (r0(j,i),j=1,6) x0=r0(1,i) y0=r0(2,i) z0=r0(3,i) r00(1,i)=x0 r00(2,i)=y0 r00(3,i)=z0 if(r0(6,i).gt.0.d0) r0(6,i)=-r0(6,i) if(r0(5,i).eq.1.d0.and.r0(6,i).eq.0.d0) goto 1 gcs=gcs+r0(4,i)*r0(4,i) gcv=gcv+r0(4,i)*r0(4,i)*r0(4,i) 1 continue close(2) gcsr=dsqrt(gcs) gcvr=gcv**(1.d0/3.d0) goto 11 10 write(6,*) 'fatal error in the input file' stop 11 k=twopi/w xv=k*gcvr xs=k*gcsr write(6,'(a,f7.3,a,f7.3)') ' volume-equiv. xv: ',xv, + ' surface-equiv. xs: ',xs write(6,'(/)') fileout='gmm01TrA.out' fileout1=fileout if(idMie.eq.1) then fileout='M'//fileout1 write(6,*) '*** Calculating only coherent Mie-scattering ***' write(6,*) '*** No interaction included ********************' endif do i=1,nL x(i)=k*r0(4,i) ref(i)=dcmplx(r0(5,i),r0(6,i)) ref0(i)=dcmplx(r0(5,i),-r0(6,i)) enddo do j=1,np do i=1,nL aMie(i,j)=0.d0 bMie(i,j)=0.d0 enddo enddo nmax0=1 do i=1,nL if(i.eq.1) goto 12 if(x(i).eq.x(i-1).and.ref(i).eq.ref(i-1)) then nmax(i)=nmax(i-1) uvmax(i)=uvmax(i-1) goto 15 endif 12 write(6,'(a,i3,a,f10.4)') 'sphere #',i, + ' individual size parameter: ',x(i) c c calculating Mie-scattering coefficients for each spheres c the ratio method of Wang and van der Hulst is used in calculating c Riccati-Bessel functions [see Wang and van der Hulst, Appl. Opt. c 30, 106 (1991), Xu, Gustafson, Giovane, Blum, and Tehranian, c Phys. Rev. E 60, 2347 (1999)] c call abMiexud(x(i),ref(i),np,NXMAX,nmax(i),an,bn,NADD, + rsr0,rsi0,rsx0,px0,w1,w2,w3,w4,eps) if(nmax(i).gt.np) then write(6,*) ' Parameter np too small, must be >',nmax(i) write(6,*) ' Please change np in gmm01f.par, recompile,' write(6,*) ' then try again' stop endif uvmax(i)=nmax(i)*(nmax(i)+2) write(6,'(a,1x,i4)') + ' Actual single-sphere expansion truncation:',nmax(i) do j=1,nmax(i) temp1=an(j) temp2=bn(j) if(j.eq.1.or.j.eq.nmax(i)) + write(6,'(i10,4e15.7)') j,temp1, + dimag(an(j)),temp2,dimag(bn(j)) enddo 15 do j=1,nmax(i) aMie(i,j)=an(j) bMie(i,j)=bn(j) enddo if(nmax(i).gt.nmax0) nmax0=nmax(i) enddo iram=0 write(6,'(/)') write(6,*) 'input sphere-positions: ' i=1 write(6,'(i5,3f14.5)') i,r0(1,1),r0(2,1),r0(3,i) i=nL write(6,'(i5,3f14.5)') i,r0(1,i),r0(2,i),r0(3,i) C------------------------------------------------------------------------ C calculating constants and Gaunt coefficients C------------------------------------------------------------------------ n0=nmax0+2 fnr(0)=0.d0 do n=1,4*(nmax0+1) fnr(n)=dsqrt(dble(n)) enddo bcof(0)=1.d0 do n=0,n0-1 bcof(n+1)=fnr(n+n+2)*fnr(n+n+1)*bcof(n)/fnr(n+1)/fnr(n+1) enddo C C the formulation used here for the calculation of Gaunt coefficients C can be found in Bruning and Lo, IEEE Trans. Anten. Prop. Ap-19, 378 C (1971) and Xu, J. Comput. Appl. Math. 85, 53 (1997), J. Comput. Phys. C 139, 137 (1998) C call cofsrd(nmax0) call cofd0(nmax0) call cofnv0(nmax0) call gau0(nmax0) C------------------------------------------------------------------------ C calculating rotational and translation coefficients C------------------------------------------------------------------------ do i=1,nL-1 do j=i+1,nL ij=(j-1)*(j-2)/2+j-i x0=r0(1,i)-r0(1,j) y0=r0(2,i)-r0(2,j) z0=r0(3,i)-r0(3,j) call carsphd(x0,y0,z0,d,xt,sphi,cphi) temp=(r0(4,i)+r0(4,j))/d if(temp.lt.1.d0) goto 151 temp1=r0(4,i) if(r0(4,j).lt.r0(4,i)) temp1=r0(4,j) temp1=(r0(4,i)+r0(4,j)-d)/temp1 if(temp1.gt.0.001d0) then write(6,'(/)') write(6,*) 'Fatal error in aggregate configuration:' write(6,*) ' SPHERES OVERLAPPED' write(6,'(e12.4)') temp1 write(6,'(i6,4e12.4)') i,(r0(ilist,i),ilist=1,4) write(6,'(i6,4e12.4)') j,(r0(ilist,j),ilist=1,4) write(6,*) 'This code calculates aggregate-scattering' write(6,*) 'of non-intersecting spheres. The ratio of' write(6,*) '(r1+r2)/d must not exceed 1, where r1+r2' write(6,*) 'is the sum of the radii of any pair of' write(6,*) 'spheres and d is the separation distance' write(6,*) 'between the two sphere centers. This ratio' write(6,*) 'for the pair of spheres listed is greater' write(6,*) 'than 1 and the overlapped portion exceeds' write(6,*) '0.001 of the (smaller) sphere-radius as ' write(6,*) 'shown above. The sequence numbers and the' write(6,*) 'Cartesian coordinates of the two sphere' write(6,*) 'centers, the radii of the two spheres are' write(6,*) 'also shown. Please correct the input data' write(6,*) 'for the aggregate configuration and then' write(6,*) 'run again.' stop endif 151 confg(1,ij)=x0 confg(2,ij)=y0 confg(3,ij)=z0 confg(4,ij)=d confg(5,ij)=temp if(temp.le.fint) goto 16 ephi=dcmplx(cphi,sphi) nlarge=max(nmax(i),nmax(j)) do m=1,nlarge ek(m,ij)=ephi**m enddo xd=k*d nbes=2*nlarge+1 call besseljd(nbes,xd,besj) call besselyd(nbes,xd,besy) C C calculating (-1)^{m+u)d_{mu}^n where d_{mu}^n is the "reduced rotation C matrix elements" C call rotcoef(xt,nlarge) irc=0 do n=1,nlarge n1=n*(n+1) do u=-n,n do m=-n,n imn=n1+m irc=irc+1 drot(irc,ij)=dc(u,imn) enddo enddo enddo itrc=0 nsmall=min(nmax(i),nmax(j)) C C the formulation used here for the calculation of vector translation C coefficients are from Cruzan, Q. Appl. Math. 20, 33 (1962) and C Xu, J. Comput. Phys. 139, 137 (1998) C do m=-nsmall,nsmall n1=max(1,iabs(m)) do n=n1,nlarge do v=n1,nlarge itrc=itrc+1 call cofxuds0(nmax0,m,n,v,besj,besy, + atr0(itrc,ij),btr0(itrc,ij), + atr1(itrc,ij),btr1(itrc,ij)) enddo enddo enddo 16 continue enddo enddo if(idMie.eq.1) then do j=1,nL do iuv=1,uvmax(j) v=dsqrt(dble(iuv)) do i=1,nL do imn=1,uvmax(i) pct(j,i,1,1,imn,iuv)=0.d0 pct(j,i,1,2,imn,iuv)=0.d0 pct(j,i,2,1,imn,iuv)=0.d0 pct(j,i,2,2,imn,iuv)=0.d0 if(j.eq.i.and.imn.eq.iuv) then pct(j,i,1,1,imn,iuv)=aMie(j,v) pct(j,i,2,2,imn,iuv)=bMie(j,v) endif enddo enddo enddo enddo goto 1800 endif C-------------------------------------------------------------------- c begins iteration process to solve T-matrix columns of (u,v,q) c BI-CGSTAB [see van der Vorst, SIAM J. Sci. Stat. Comput. 13, 631 c (1992); Gutknecht, SIAM J. Sci. comput. 14, pp. 1020-1033 (1993)] C-------------------------------------------------------------------- do i=1,nL-1 do j=i+1,nL ij=(j-1)*(j-2)/2+j-i d=confg(4,ij) xd=k*d bes0(ij)=dsin(xd)/xd enddo enddo write(6,*) 'Starting Bi-CGSTAB to solve T-matrix' n0=nmax0*(nmax0+2) do 1001 iuv=1,n0 v=dsqrt(dble(iuv)) iuvc=v*v c iuvc=1 u=iuv-v*(v+1) do 1002 iq=1,2 iuv1=-u+v+v*v do imn=1,iuvc-1 n=dsqrt(dble(imn)) m=imn-n*n-n imn1=-m+n*n+n cz=(-1)**(m+u+n+v) do j=1,nL do i=1,nL asr(j,i,imn)=cz*pct(i,j,iq,1,iuv1,imn1) bsr(j,i,imn)=cz*pct(i,j,iq,2,iuv1,imn1) enddo enddo enddo do j=1,nL ind(j)=0 do i=1,nL do imn=1,iuvc-1 as(i,imn)=asr(j,i,imn) bs(i,imn)=bsr(j,i,imn) enddo do imn=iuvc,uvmax(i) as(i,imn)=ci0 bs(i,imn)=ci0 enddo enddo call transT(nL,r0,nmax,uvmax,fint,atr0,btr0, + ek,drot,as,bs,as2,bs2,ind,confg,iuvc,2) do i=1,nL do imn=iuvc,uvmax(i) tta0(j,i,imn)=as2(i,imn) ttb0(j,i,imn)=bs2(i,imn) enddo enddo enddo do j=1,nL do i=1,nL do imn=iuvc,n0 asr(j,i,imn)=ci0 bsr(j,i,imn)=ci0 enddo enddo enddo do j=1,nL if(iuv.le.uvmax(j)) then if(iq.eq.1) then asr(j,j,iuv)=aMie(j,v) else bsr(j,j,iuv)=bMie(j,v) endif endif enddo do j=1,nL if(iq.eq.1) then c0i(j)=aMie(j,v)*dconjg(aMie(j,v)) else c0i(j)=bMie(j,v)*dconjg(bMie(j,v)) endif enddo niter=1 do j=1,nL do i=1,nL do imn=1,iuvc-1 as(i,imn)=ci0 bs(i,imn)=ci0 enddo do imn=iuvc,uvmax(i) as(i,imn)=asr(j,i,imn) bs(i,imn)=bsr(j,i,imn) enddo enddo call transT(nL,r0,nmax,uvmax,fint,atr0,btr0, + ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1) do i=1,nL do imn=iuvc,uvmax(i) tta(j,i,imn)=as2(i,imn)+tta0(j,i,imn) ttb(j,i,imn)=bs2(i,imn)+ttb0(j,i,imn) enddo enddo enddo do 611 i=1,nL c1i(i)=0.d0 do 6111 imn=iuvc,uvmax(i) n=dsqrt(dble(imn)) do 6112 j=1,nL if(iuv.gt.uvmax(j)) goto 6112 as1(j,i,imn)=0.d0 bs1(j,i,imn)=0.d0 if(imn.eq.iuv) then ijmax=max(i,j) ijmin=min(i,j) if(ijmax.eq.ijmin) then cz=1.d0 else ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin cz=bes0(ij) endif if(iq.eq.1) then as1(j,i,imn)=as1(j,i,imn)+cz*aMie(i,n) else bs1(j,i,imn)=bs1(j,i,imn)+cz*bMie(i,n) endif endif do 6113 jj=1,nL ijmax=max(jj,j) ijmin=min(jj,j) if(ijmax.eq.ijmin) then cz=1.d0 else ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin cz=bes0(ij) endif A0=asr(jj,i,imn) B0=bsr(jj,i,imn) A0=A0+aMie(i,n)*tta(jj,i,imn) B0=B0+bMie(i,n)*ttb(jj,i,imn) as1(j,i,imn)=as1(j,i,imn)-cz*A0 bs1(j,i,imn)=bs1(j,i,imn)-cz*B0 6113 continue A=as1(j,i,imn) B=bs1(j,i,imn) c1i(i)=c1i(i)+A*dconjg(A) c1i(i)=c1i(i)+B*dconjg(B) 6112 continue 6111 continue 611 continue temp=0.d0 B0=0.d0 do 612 i=1,nL cext0=c1i(i)/c0i(i) if(cext0.lt.small) ind(i)=1 if(ind(i).gt.0) goto 612 if(cext0.gt.temp) temp=cext0 B0=B0+c1i(i) 612 continue if(temp.lt.small) then do i=1,nL do imn=1,uvmax(i) do j=1,nL pct(j,i,1,iq,imn,iuv)=asr(j,i,imn) pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn) enddo enddo enddo goto 1002 endif do 613 i=1,nL if(ind(i).gt.0) goto 613 do imn=1,iuvc-1 do j=1,nL asp(j,i,imn)=ci0 bsp(j,i,imn)=ci0 as0(j,i,imn)=ci0 bs0(j,i,imn)=ci0 enddo enddo do imn=iuvc,uvmax(i) do j=1,nL asp(j,i,imn)=as1(j,i,imn) bsp(j,i,imn)=bs1(j,i,imn) as0(j,i,imn)=as1(j,i,imn) bs0(j,i,imn)=bs1(j,i,imn) enddo enddo 613 continue do j=1,nL do i=1,nL do imn=1,iuvc-1 as(i,imn)=ci0 bs(i,imn)=ci0 enddo do imn=iuvc,uvmax(i) as(i,imn)=asp(j,i,imn) bs(i,imn)=bsp(j,i,imn) enddo enddo call transT(nL,r0,nmax,uvmax,fint,atr0,btr0, + ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1) do i=1,nL do imn=iuvc,uvmax(i) tta(j,i,imn)=as2(i,imn) ttb(j,i,imn)=bs2(i,imn) enddo enddo enddo A0=0.d0 do 614 i=1,nL if(ind(i).gt.0) goto 614 do 6141 imn=iuvc,uvmax(i) n=dsqrt(dble(imn)) do 6142 j=1,nL if(iuv.gt.uvmax(j)) goto 6142 ast(j,i,imn)=0.d0 bst(j,i,imn)=0.d0 do 6143 jj=1,nL ijmax=max(jj,j) ijmin=min(jj,j) if(ijmax.eq.ijmin) then cz=1.d0 else ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin cz=bes0(ij) endif Aj2=asp(jj,i,imn) Bj2=bsp(jj,i,imn) Aj2=Aj2+aMie(i,n)*tta(jj,i,imn) Bj2=Bj2+bMie(i,n)*ttb(jj,i,imn) ast(j,i,imn)=ast(j,i,imn)+cz*Aj2 bst(j,i,imn)=bst(j,i,imn)+cz*Bj2 6143 continue A0=A0+dconjg(as0(j,i,imn))*ast(j,i,imn) A0=A0+dconjg(bs0(j,i,imn))*bst(j,i,imn) 6142 continue 6141 continue 614 continue if(cdabs(A0).lt.1.d-200) then do i=1,nL do imn=1,uvmax(i) do j=1,nL pct(j,i,1,iq,imn,iuv)=asr(j,i,imn) pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn) enddo enddo enddo goto 1002 endif Aj=B0/A0 62 do 621 i=1,nL if(ind(i).gt.0) goto 621 do imn=iuvc,uvmax(i) do 6211 j=1,nL if(iuv.gt.uvmax(j)) then asv(j,i,imn)=0.d0 bsv(j,i,imn)=0.d0 go to 6211 endif asv(j,i,imn)=asp(j,i,imn)-Aj*ast(j,i,imn) bsv(j,i,imn)=bsp(j,i,imn)-Aj*bst(j,i,imn) 6211 continue enddo 621 continue do j=1,nL do i=1,nL do imn=1,iuvc-1 as(i,imn)=ci0 bs(i,imn)=ci0 enddo do imn=iuvc,uvmax(i) as(i,imn)=asv(j,i,imn) bs(i,imn)=bsv(j,i,imn) enddo enddo call transT(nL,r0,nmax,uvmax,fint,atr0,btr0, + ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1) do i=1,nL do imn=iuvc,uvmax(i) tta(j,i,imn)=as2(i,imn) ttb(j,i,imn)=bs2(i,imn) enddo enddo enddo A2=0.d0 B2=0.d0 do 622 i=1,nL if(ind(i).gt.0) goto 622 do 6221 imn=iuvc,uvmax(i) n=dsqrt(dble(imn)) do 6222 j=1,nL if(iuv.gt.uvmax(j)) goto 6222 asc(j,i,imn)=0.d0 bsc(j,i,imn)=0.d0 do 6223 jj=1,nL ijmax=max(jj,j) ijmin=min(jj,j) if(ijmax.eq.ijmin) then cz=1.d0 else ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin cz=bes0(ij) endif Aj2=asv(jj,i,imn) Bj2=bsv(jj,i,imn) Aj2=Aj2+aMie(i,n)*tta(jj,i,imn) Bj2=Bj2+bMie(i,n)*ttb(jj,i,imn) asc(j,i,imn)=asc(j,i,imn)+cz*Aj2 bsc(j,i,imn)=bsc(j,i,imn)+cz*Bj2 6223 continue A2=A2+dconjg(asc(j,i,imn))*asv(j,i,imn) A2=A2+dconjg(bsc(j,i,imn))*bsv(j,i,imn) B2=B2+dconjg(asc(j,i,imn))*asc(j,i,imn) B2=B2+dconjg(bsc(j,i,imn))*bsc(j,i,imn) 6222 continue 6221 continue 622 continue if(cdabs(B2).lt.1.d-200) then do i=1,nL do imn=1,uvmax(i) do j=1,nL pct(j,i,1,iq,imn,iuv)=asr(j,i,imn) pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn) enddo enddo enddo goto 1002 endif Bj=A2/B2 do 623 i=1,nL if(ind(i).gt.0) goto 623 do imn=iuvc,uvmax(i) do 6231 j=1,nL if(iuv.gt.uvmax(j)) then asp(j,i,imn)=ci0 bsp(j,i,imn)=ci0 go to 6231 endif asp(j,i,imn)=asv(j,i,imn)-Bj*asc(j,i,imn) bsp(j,i,imn)=bsv(j,i,imn)-Bj*bsc(j,i,imn) 6231 continue enddo 623 continue do 624 i=1,nL c1i(i)=0.d0 if(ind(i).gt.0) goto 624 do 6241 imn=iuvc,uvmax(i) do 6242 j=1,nL if(iuv.gt.uvmax(j)) goto 6242 Aj2=Aj*as1(j,i,imn)+Bj*asv(j,i,imn) Bj2=Aj*bs1(j,i,imn)+Bj*bsv(j,i,imn) asr(j,i,imn)=asr(j,i,imn)+Aj2 bsr(j,i,imn)=bsr(j,i,imn)+Bj2 c1i(i)=c1i(i)+Aj2*dconjg(Aj2) c1i(i)=c1i(i)+Bj2*dconjg(Bj2) 6242 continue 6241 continue 624 continue cext0=0.d0 cext1=0.d0 do 625 i=1,nL if(ind(i).gt.0) goto 625 cext0=cext0+c0i(i) cext1=cext1+c1i(i) 625 continue temp=cext1/cext0 if(iuv.eq.1.and.iq.eq.1) then if(niter/20*20.eq.niter.or.niter.eq.1) then write(6,'(a11,i4,2x,e15.7)') + 'iteration #',niter,temp endif endif if(temp.lt.small) then if(u.eq.v) then write(6,'(3i4,a12,i4,2x,e15.5)') u,v,iq, + ' iteration ',niter,temp endif do i=1,nL do imn=1,uvmax(i) do j=1,nL pct(j,i,1,iq,imn,iuv)=asr(j,i,imn) pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn) enddo enddo enddo goto 1002 endif if(niter.gt.MXINT) then write(6,'(a9,i5,i4,i3)') 'Caution: ',u,v,iq write(6,*) '*** Maximum iterations exceeded ***' write(6,*) '*** Solution may be inaccurate ***' goto 1002 endif B2=0.d0 do 626 i=1,nL if(ind(i).gt.0) goto 626 B2i(i)=0.d0 do imn=iuvc,uvmax(i) do 6261 j=1,nL if(iuv.gt.uvmax(j)) goto 6261 Aj2=dconjg(as0(j,i,imn))*asp(j,i,imn) Bj2=dconjg(bs0(j,i,imn))*bsp(j,i,imn) B2i(i)=B2i(i)+Aj2 B2i(i)=B2i(i)+Bj2 6261 continue enddo B2=B2+B2i(i) 626 continue A0=B0*Bj if(cdabs(A0).lt.1.d-200) then do i=1,nL do imn=1,uvmax(i) do j=1,nL pct(j,i,1,iq,imn,iuv)=asr(j,i,imn) pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn) enddo enddo enddo goto 1002 endif A0=-Aj*B2/A0 do 627 i=1,nL if(ind(i).gt.0) goto 627 do imn=iuvc,uvmax(i) do 6271 j=1,nL if(iuv.gt.uvmax(j)) goto 6271 as1(j,i,imn)=as1(j,i,imn)-Bj*ast(j,i,imn) bs1(j,i,imn)=bs1(j,i,imn)-Bj*bst(j,i,imn) as1(j,i,imn)=asp(j,i,imn)-A0*as1(j,i,imn) bs1(j,i,imn)=bsp(j,i,imn)-A0*bs1(j,i,imn) 6271 continue enddo 627 continue B0=B2 do j=1,nL do i=1,nL do imn=1,iuvc-1 as(i,imn)=ci0 bs(i,imn)=ci0 enddo do imn=iuvc,uvmax(i) as(i,imn)=as1(j,i,imn) bs(i,imn)=bs1(j,i,imn) enddo enddo call transT(nL,r0,nmax,uvmax,fint,atr0,btr0, + ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1) do i=1,nL do imn=iuvc,uvmax(i) tta(j,i,imn)=as2(i,imn) ttb(j,i,imn)=bs2(i,imn) enddo enddo enddo A0=0.d0 do 629 i=1,nL if(ind(i).gt.0) goto 629 do 6291 imn=iuvc,uvmax(i) n=dsqrt(dble(imn)) do 6292 j=1,nL if(iuv.gt.uvmax(j)) goto 6292 ast(j,i,imn)=0.d0 bst(j,i,imn)=0.d0 do 6293 jj=1,nL ijmax=max(jj,j) ijmin=min(jj,j) if(ijmax.eq.ijmin) then cz=1.d0 else ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin cz=bes0(ij) endif Aj2=as1(jj,i,imn) Bj2=bs1(jj,i,imn) Aj2=Aj2+aMie(i,n)*tta(jj,i,imn) Bj2=Bj2+bMie(i,n)*ttb(jj,i,imn) ast(j,i,imn)=ast(j,i,imn)+cz*Aj2 bst(j,i,imn)=bst(j,i,imn)+cz*Bj2 6293 continue A0=A0+dconjg(as0(j,i,imn))*ast(j,i,imn) A0=A0+dconjg(bs0(j,i,imn))*bst(j,i,imn) 6292 continue 6291 continue 629 continue if(cdabs(A0).lt.1.d-200) then do i=1,nL do imn=1,uvmax(i) do j=1,nL pct(j,i,1,iq,imn,iuv)=asr(j,i,imn) pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn) enddo enddo enddo goto 1002 endif Aj=B0/A0 niter=niter+1 goto 62 1002 continue 1001 continue C------------------------------------------------------------------------ C calculating random-orientation averaged total and individual-particle C extinction cross-sections C------------------------------------------------------------------------ 1800 do i=1,nL ind(i)=0 cexti(i)=0.d0 cscai(i)=0.d0 cpri(i)=0.d0 enddo cext=0.d0 csca=0.d0 cpr=0.d0 n0=nmax0*(nmax0+2) do iuv=1,n0 do iq=1,2 do 1801 j=1,nL if(iuv.gt.uvmax(j)) goto 1801 do i=1,nL do imn=1,uvmax(i) as(i,imn)=pct(j,i,1,iq,imn,iuv) bs(i,imn)=pct(j,i,2,iq,imn,iuv) enddo enddo call transT(nL,r0,nmax,uvmax,fint,atr1,btr1, + ek,drot,as,bs,as2,bs2,ind,confg,1,1) do i=1,nL do imn=1,uvmax(i) A=as2(i,imn)+as(i,imn) B=bs2(i,imn)+bs(i,imn) pct(j,i,1,iq,imn,iuv)=A pct(j,i,2,iq,imn,iuv)=B enddo enddo cz=pct(j,j,iq,iq,iuv,iuv) cext=cext+cz cexti(j)=cexti(j)+cz 1801 continue enddo enddo C------------------------------------------------------------------------ C calculating random-orientation averaged asymmetry parameter and C total and individual-particle scattering cross-sections C------------------------------------------------------------------------ do j=1,nL do iuv=1,uvmax(j) v=dsqrt(dble(iuv)) u=iuv-v*v-v iuv1=iuv-2*u cb=(-1)**(v+u) jv1=v+1 iuv2=-u+jv1*(jv1+1) iuv3=-u+(v-1)*v fuv1=v*jv1 fuv1=dble(u)/fuv1 fuv2=fnr(v)*fnr(v+2)/fnr(2*v+1) fuv2=fuv2/fnr(2*v+3)/dble(jv1) fuv2=fuv2*fnr(jv1-u)*fnr(jv1+u) fuv3=fnr(v-1)*fnr(jv1)/fnr(2*v-1) fuv3=fuv3/fnr(2*v+1)/dble(v) fuv3=fuv3*fnr(v-u)*fnr(v+u) juv1=iuv1-1 juv2=iuv1+1 juv3=iuv2-1 juv4=iuv3-1 juv5=iuv2+1 juv6=iuv3+1 guv=v*(v+1) guv1=-fnr(v-u)*fnr(v+u+1)/guv guv2=-fnr(v+u)*fnr(v-u+1)/guv guv=dble(v+1)*fnr(v)*fnr(v+2)*fnr(2*v+1)*fnr(2*v+3) guv3=u*(u+1)+((v-u)*(v+u+3)+(v+u)*(v-u+1))/2 guv3=guv3*fnr(v+u+1)*fnr(v+u+2)/guv guv5=u*(u-1)+((v-u)*(v+u+1)+(v+u)*(v-u+3))/2 guv5=-guv5*fnr(v-u+1)*fnr(v-u+2)/guv if(v.gt.1) then guv=dble(v)*fnr(v-1)*fnr(v+1)*fnr(2*v-1)*fnr(2*v+1) guv4=u*(u+1)+((v-u-2)*(v+u+1)+(v+u)*(v-u+1))/2 ntemp=v-u-1 if(ntemp.lt.0) then temp=0.d0 else temp=fnr(ntemp) endif guv4=-guv4*fnr(v-u)*temp/guv guv6=u*(u-1)+((v-u)*(v+u+1)+(v+u-2)*(v-u+1))/2 ntemp=v+u-1 if(ntemp.lt.0) then temp=0.d0 else temp=fnr(ntemp) endif guv6=guv6*fnr(v+u)*temp/guv endif do i=1,nL do imn=1,uvmax(i) n=dsqrt(dble(imn)) m=imn-n*n-n imn1=-m+n*n+n sb=(-1)**(m+n) sb=cb*sb n1=n+1 n2=2*n rn=1.0d0/dble(n*n1) p=fnr(n)*fnr(n+2)/fnr(n2+1) p=p/fnr(n2+3)/dble(n1) t=fnr(n-1)*fnr(n+1)/fnr(n2-1) t=t/fnr(n2+1)/dble(n) rm=dble(m)*rn imn2=(n+1)*(n+2)+m fnp=fnr(n+m+1)*fnr(n-m+1)*p imn3=(n-1)*n+m fn=fnr(n+m)*fnr(n-m)*t jmn1=imn+1 jmn2=imn-1 jmn3=imn2+1 jmn4=imn3+1 jmn5=imn2-1 jmn6=imn3-1 gmn=n*(n+1) gmn1=-fnr(n-m)*fnr(n+m+1)/gmn gmn2=-fnr(n+m)*fnr(n-m+1)/gmn gmn=dble(n+1)*fnr(n)*fnr(n+2) gmn=gmn*fnr(2*n+1)*fnr(2*n+3) gmn3=m*(m+1)+((n-m)*(n+m+3)+(n+m)*(n-m+1))/2 gmn3=gmn3*fnr(n+m+1)*fnr(n+m+2)/gmn gmn5=m*(m-1)+((n-m)*(n+m+1)+(n+m)*(n-m+3))/2 gmn5=-gmn5*fnr(n-m+1)*fnr(n-m+2)/gmn if(n.gt.1) then gmn=dble(n)*fnr(n-1)*fnr(n+1) gmn=gmn*fnr(2*n-1)*fnr(2*n+1) gmn4=m*(m+1)+((n-m-2)*(n+m+1)+(n+m)*(n-m+1))/2 ntemp=n-m-1 if(ntemp.lt.0) then temp=0.d0 else temp=fnr(ntemp) endif gmn4=-gmn4*fnr(n-m)*temp/gmn gmn6=m*(m-1)+((n-m)*(n+m+1)+(n+m-2)*(n-m+1))/2 ntemp=n+m-1 if(ntemp.lt.0) then temp=0.d0 else temp=fnr(ntemp) endif gmn6=gmn6*fnr(n+m)*temp/gmn endif do ip=1,2 do iq=1,2 A=pct(j,i,ip,iq,imn,iuv) B=sb*pct(i,j,iq,ip,iuv1,imn1) B=dconjg(B) cz=B*A csca=csca+cz cscai(j)=cscai(j)+cz A=rm*pct(j,i,3-ip,iq,imn,iuv) if(n.eq.nmax(i)) goto 1951 A=A+fnp*pct(j,i,ip,iq,imn2,iuv) 1951 if(n.eq.1.or.iabs(m).gt.n-1) goto 1952 A=A+fn*pct(j,i,ip,iq,imn3,iuv) 1952 B=fuv1*pct(i,j,3-iq,ip,iuv1,imn1) if(v.eq.nmax(j)) goto 1953 B=B-fuv2*pct(i,j,iq,ip,iuv2,imn1) 1953 if(v.eq.1.or.iabs(u).gt.v-1) goto 1954 B=B-fuv3*pct(i,j,iq,ip,iuv3,imn1) 1954 B=sb*dconjg(B) cz=B*A cpr=cpr+cz cpri(j)=cpri(j)+cz A=ci0 if(iabs(m+1).gt.n) goto 1961 A=gmn1*pct(j,i,3-ip,iq,jmn1,iuv) 1961 if(n.eq.nmax(i)) goto 1962 A=A+gmn3*pct(j,i,ip,iq,jmn3,iuv) 1962 if(n.eq.1) goto 1963 if(iabs(m+1).gt.n-1) goto 1963 A=A+gmn4*pct(j,i,ip,iq,jmn4,iuv) 1963 B=ci0 if(iabs(u-1).gt.v) goto 1964 B=-guv2*pct(i,j,3-iq,ip,juv2,imn1) 1964 if(v.eq.nmax(j)) goto 1965 B=B+guv5*pct(i,j,iq,ip,juv5,imn1) 1965 if(v.eq.1) goto 1966 if(iabs(u-1).gt.v-1) goto 1966 B=B+guv6*pct(i,j,iq,ip,juv6,imn1) 1966 cz=dconjg(B)*A temp=0.5d0*sb*cz cpr=cpr+temp cpri(j)=cpri(j)+temp A=ci0 if(iabs(m-1).gt.n) goto 1971 A=gmn2*pct(j,i,3-ip,iq,jmn2,iuv) 1971 if(n.eq.nmax(i)) goto 1972 A=A+gmn5*pct(j,i,ip,iq,jmn5,iuv) 1972 if(n.eq.1) goto 1973 if(iabs(m-1).gt.n-1) goto 1973 A=A+gmn6*pct(j,i,ip,iq,jmn6,iuv) 1973 B=ci0 if(iabs(u+1).gt.v) goto 1974 B=-guv1*pct(i,j,3-iq,ip,juv1,imn1) 1974 if(v.eq.nmax(j)) goto 1975 B=B+guv3*pct(i,j,iq,ip,juv3,imn1) 1975 if(v.eq.1) goto 1976 if(iabs(u+1).gt.v-1) goto 1976 B=B+guv4*pct(i,j,iq,ip,juv4,imn1) 1976 cz=dconjg(B)*A temp=0.5d0*sb*cz cpr=cpr+temp cpri(j)=cpri(j)+temp enddo enddo enddo enddo enddo enddo if(idscmt.lt.0) goto 2000 C------------------------------------------------------------------------ C calculating random-orientation averaged Mueller matrix elements C------------------------------------------------------------------------ write(6,'(/)') write(6,*) 'Calculating (Mueller) scattering matrix' write(6,'(/)') n0=nmax0*(nmax0+2) nmax2=2*nmax0 do n=1,nmax0 do v=1,nmax0 do m=0,n it=n+v+1 call xuwigd(n,v,m,-1,w01s,wcf,it,nmf1) do ii=1,nmf1 it=it-1 wmf1(n,v,m,it)=w01s(ii) enddo it=n+v+1 call xuwigd(n,v,m,1,w01s,wcf,it,nm1) do ii=1,nm1 it=it-1 wm1(n,v,m,it)=w01s(ii) enddo enddo enddo enddo do n=1,nmax0 do jn=1,nmax0 do ip=1,2 do jp=1,2 B11n(n,jn,ip,jp)=ci0 B12n(n,jn,ip,jp)=ci0 B13n(n,jn,ip,jp)=ci0 B21n(n,jn,ip,jp)=ci0 B22n(n,jn,ip,jp)=ci0 B23n(n,jn,ip,jp)=ci0 do m=0,nmax0 A1m(m,n,jn,ip,jp)=ci0 A2m(m,n,jn,ip,jp)=ci0 enddo enddo enddo enddo enddo do 1900 ids=-nmax2,nmax2 do n=1,nmax0 do 19001 v=1,nmax0 if(iabs(ids).gt.n+v) goto 19001 fv0=0.5d0*fnr(2*v+1) do 19002 ms=-n,n mw=ms+ids if(iabs(mw).gt.v) goto 19002 it=n+v+1 call xuwigd(n,v,ms,-mw,w01s,wcf,it,nvs) fv=((-1)**ms)*fv0 do ii=1,nvs it=it-1 wsdt(n,v,ms,it)=fv*w01s(ii) enddo 19002 continue 19001 continue enddo do i=1,nL do j=1,nL do it=0,nmax2 do n=1,nmax0 do m=0,n do ip=1,2 fhmf1(it,n,m,ip)=ci0 fmf1(it,n,m,ip)=ci0 fhm1v(it,n,m,ip)=ci0 fm1v(it,n,m,ip)=ci0 fhm1q(it,n,m,ip)=ci0 fm1q(it,n,m,ip)=ci0 fhmf1vq(it,n,m,ip)=ci0 fmf1vq(it,n,m,ip)=ci0 enddo enddo enddo enddo do 19003 n=1,nmax0 if(n.gt.nmax(i)) goto 19003 in0=n*(n+1) cv=1.d0 do 19004 v=1,nmax0 if(v.gt.nmax(j)) goto 19004 if(iabs(ids).gt.n+v) goto 19004 cv=-cv iv0=v*(v+1) do it=0,n+v do ip=1,2 do iq=1,2 fhas(it,ip,iq)=ci0 fnhs(it,ip,iq)=ci0 enddo enddo enddo do 1901 ms=-n,n mw=ms+ids if(iabs(mw).gt.v) goto 1901 it=n+v+1 fn=(-1)**(n+v+ids) isn=in0+ms isf=in0-ms iwv=iv0+mw iwf=iv0-mw nvs=n+v-max(iabs(n-v),iabs(ids))+1 do ii=1,nvs it=it-1 do ip=1,2 do iq=1,2 A=pct(i,j,iq,ip,iwf,isf) A=fn*wsdt(n,v,ms,it)*A fhas(it,ip,iq)= + fhas(it,ip,iq)+A B=pct(j,i,ip,iq,isn,iwv) B=wsdt(n,v,ms,it)*B fnhs(it,ip,iq)= + fnhs(it,ip,iq)+B enddo enddo enddo 1901 continue itmin=max(iabs(ids),iabs(n-v)) do it=itmin,n+v do m=0,n if(iabs(m-1).gt.it) then cwmf1=0.d0 else cwmf1=wmf1(n,v,m,it) endif if(m+1.gt.it) then cwm1=0.d0 else cwm1=wm1(n,v,m,it) endif do ip=1,2 do iq=1,2 cq=(-1)**iq A=cwmf1*fhas(it,ip,iq) A=dconjg(A) B=cv*A fhmf1(it,n,m,ip)= + fhmf1(it,n,m,ip)+A fhmf1vq(it,n,m,ip)= + fhmf1vq(it,n,m,ip)+cq*B A=cwm1*fhas(it,ip,iq) A=dconjg(A) B=cv*A fhm1v(it,n,m,ip)= + fhm1v(it,n,m,ip)+B fhm1q(it,n,m,ip)= + fhm1q(it,n,m,ip)+cq*A A=cwmf1*fnhs(it,ip,iq) B=cv*A fmf1(it,n,m,ip)= + fmf1(it,n,m,ip)+A fmf1vq(it,n,m,ip)= + fmf1vq(it,n,m,ip)+cq*B A=cwm1*fnhs(it,ip,iq) B=cv*A fm1v(it,n,m,ip)= + fm1v(it,n,m,ip)+B fm1q(it,n,m,ip)= + fm1q(it,n,m,ip)+cq*A enddo enddo enddo enddo 19004 continue 19003 continue do 1904 n=1,nmax0 if(n.gt.nmax(i)) goto 1904 itmin=iabs(ids) if(itmin.gt.n+nmax0) goto 1904 do 1905 jn=1,nmax0 if(jn.gt.nmax(j)) goto 1905 if(itmin.gt.jn+nmax0) goto 1905 itmax=min(n,jn)+nmax0 do 1906 it=itmin,itmax gt=2*it+1 do m=0,min(n,jn) do ip=1,2 do jp=1,2 A=fhmf1(it,n,m,ip) B=fmf1(it,jn,m,jp) A0=A*B A=fhm1q(it,n,m,ip) B=fm1q(it,jn,m,jp) A=gt*(A0+A*B) A1m(m,n,jn,ip,jp)= + A1m(m,n,jn,ip,jp)+A fn=(-1)**(ip+jp+n+jn) A=fhmf1vq(it,n,m,ip) B=fmf1vq(it,jn,m,jp) A0=A*B A=fhm1v(it,n,m,ip) B=fm1v(it,jn,m,jp) A=fn*gt*(A0+A*B) A2m(m,n,jn,ip,jp)= + A2m(m,n,jn,ip,jp)+A if(m.gt.0) goto 19061 temp=((-1)**ip)*gt A=fhm1q(it,n,0,ip) B=fmf1(it,jn,2,jp) A=temp*A*B B11n(n,jn,ip,jp)= + B11n(n,jn,ip,jp)+A nj=ip+n+it temp=((-1)**nj)*gt A=fhmf1vq(it,n,1,ip) B=fmf1(it,jn,1,jp) A=temp*A*B B12n(n,jn,ip,jp)= + B12n(n,jn,ip,jp)+A nj=ip+n+jn temp=((-1)**nj)*gt A=fhmf1vq(it,n,2,ip) B=fm1v(it,jn,0,jp) A=temp*A*B B13n(n,jn,ip,jp)= + B13n(n,jn,ip,jp)+A nj=jp+n+jn temp=((-1)**nj)*gt A=fhm1v(it,n,0,ip) B=fmf1vq(it,jn,2,jp) A=temp*A*B B21n(n,jn,ip,jp)= + B21n(n,jn,ip,jp)+A nj=jp+jn+it temp=((-1)**nj)*gt A=fhmf1(it,n,1,ip) B=fmf1vq(it,jn,1,jp) A=temp*A*B B22n(n,jn,ip,jp)= + B22n(n,jn,ip,jp)+A temp=((-1)**jp)*gt A=fhmf1(it,n,2,ip) B=fm1q(it,jn,0,jp) A=temp*A*B B23n(n,jn,ip,jp)= + B23n(n,jn,ip,jp)+A 19061 continue enddo enddo enddo 1906 continue 1905 continue 1904 continue enddo enddo 1900 continue do 1910 ia=1,nang iang=2*nang-ia dang(ia)=sang*dble(ia-1) dang(iang)=180.0d0-dang(ia) theta=dang(ia)*pione/180.0d0 xt=dcos(theta) call tipitaud(nmax0,xt) do i=1,4 do j=1,4 mue(i,j,ia)=0.d0 mue(i,j,iang)=0.d0 enddo enddo do i=1,2 do j=1,2 do ik=1,2 do jk=1,2 A0p(i,j,ik,jk)=ci0 A1p(i,j,ik,jk)=ci0 A0pg(i,j,ik,jk)=ci0 A1pg(i,j,ik,jk)=ci0 B0p(i,j,ik,jk)=ci0 B1p(i,j,ik,jk)=ci0 B0pg(i,j,ik,jk)=ci0 B1pg(i,j,ik,jk)=ci0 enddo enddo enddo enddo do n=1,nmax0 gn=(-1)**n itau0=(n-1)*(n+2)/2+1 tau0p(1)=0.5d0*tau(itau0) tau0p(2)=0.d0 tau0pg(1)=-gn*tau0p(1) tau0pg(2)=0.d0 itau=itau0+1 gn=-gn tau1p(1)=tau(itau) tau1p(2)=pi(itau) tau1pg(1)=-gn*tau1p(1) tau1pg(2)=gn*tau1p(2) itau=itau0+2 gn=-gn tau2p(1)=tau(itau) tau2p(2)=pi(itau) tau2pg(1)=-gn*tau2p(1) tau2pg(2)=gn*tau2p(2) do jn=1,nmax0 jtau0=(jn-1)*(jn+2)/2+1 gn=(-1)**jn tau0pj(1)=0.5d0*tau(jtau0) tau0pj(2)=0.d0 tau0pjg(1)=-gn*tau0pj(1) tau0pjg(2)=0.d0 jtau=jtau0+1 gn=-gn tau1pj(1)=tau(jtau) tau1pj(2)=pi(jtau) tau1pjg(1)=-gn*tau1pj(1) tau1pjg(2)=gn*tau1pj(2) jtau=jtau0+2 gn=-gn tau2pj(1)=tau(jtau) tau2pj(2)=pi(jtau) tau2pjg(1)=-gn*tau2pj(1) tau2pjg(2)=gn*tau2pj(2) do ik=1,2 do jk=1,2 tau20(ik,jk)=tau2p(ik)*tau0pj(jk) tau11(ik,jk)=tau1p(ik)*tau1pj(jk) tau02(ik,jk)=tau0p(ik)*tau2pj(jk) tau20g(ik,jk)=tau2pg(ik)*tau0pjg(jk) tau11g(ik,jk)=tau1pg(ik)*tau1pjg(jk) tau02g(ik,jk)=tau0pg(ik)*tau2pjg(jk) enddo enddo do ip=1,2 do jp=1,2 A=B11n(n,jn,ip,jp) B=B12n(n,jn,ip,jp) cmz=B13n(n,jn,ip,jp) Aj=B21n(n,jn,ip,jp) Bj=B22n(n,jn,ip,jp) cmzj=B23n(n,jn,ip,jp) do ik=1,2 do jk=1,2 A1=A*tau02(ik,jk) A2=B*tau11(ik,jk) Aj2=cmz*tau20(ik,jk) B1=A1-A2+Aj2 A1=Aj*tau02(ik,jk) A2=Bj*tau11(ik,jk) Aj2=cmzj*tau20(ik,jk) B2=A1-A2+Aj2 A0=B1+B2 B0p(ip,jp,ik,jk)=B0p(ip,jp,ik,jk)+A0 A0=B1-B2 B1p(ip,jp,ik,jk)=B1p(ip,jp,ik,jk)+A0 A1=A*tau02g(ik,jk) A2=B*tau11g(ik,jk) Aj2=cmz*tau20g(ik,jk) B1=A1-A2+Aj2 A1=Aj*tau02g(ik,jk) A2=Bj*tau11g(ik,jk) Aj2=cmzj*tau20g(ik,jk) B2=A1-A2+Aj2 A0=B1+B2 B0pg(ip,jp,ik,jk)=B0pg(ip,jp,ik,jk)+A0 A0=B1-B2 B1pg(ip,jp,ik,jk)=B1pg(ip,jp,ik,jk)+A0 enddo enddo enddo enddo do m=0,min(n,jn) itau=itau0+m gmn=(-1)**(m+n) taup(1)=tau(itau) taup(2)=pi(itau) if(m.eq.0) then taup(1)=0.5d0*taup(1) taup(2)=0.d0 endif taupg(1)=-gmn*taup(1) taupg(2)=gmn*taup(2) gmnj=(-1)**(jn+m) jtau=jtau0+m taupj(1)=tau(jtau) taupj(2)=pi(jtau) if(m.eq.0) then taupj(1)=0.5d0*taupj(1) taupj(2)=0.d0 endif taupjg(1)=-gmnj*taupj(1) taupjg(2)=gmnj*taupj(2) do ik=1,2 do jk=1,2 taum(ik,jk)=taup(ik)*taupj(jk) taumg(ik,jk)=taupg(ik)*taupjg(jk) enddo enddo do ip=1,2 do jp=1,2 A=A1m(m,n,jn,ip,jp) Aj=A B=A2m(m,n,jn,ip,jp) A=A+B Aj=Aj-B do ik=1,2 do jk=1,2 B0=A*taum(ik,jk) B1=Aj*taum(ik,jk) A0p(ip,jp,ik,jk)= + A0p(ip,jp,ik,jk)+B0 A1p(ip,jp,ik,jk)= + A1p(ip,jp,ik,jk)+B1 B0=A*taumg(ik,jk) B1=Aj*taumg(ik,jk) A0pg(ip,jp,ik,jk)= + A0pg(ip,jp,ik,jk)+B0 A1pg(ip,jp,ik,jk)= + A1pg(ip,jp,ik,jk)+B1 enddo enddo enddo enddo enddo enddo enddo do ip=1,2 do jp=1,2 temp=A0p(ip,jp,ip,jp)+A0p(ip,jp,3-ip,3-jp) mue(1,1,ia)=mue(1,1,ia)+temp temp=B0p(ip,jp,3-ip,3-jp)-B0p(ip,jp,ip,jp) mue(1,2,ia)=mue(1,2,ia)+temp A0=A1p(ip,jp,ip,jp)+A1p(ip,jp,3-ip,3-jp) A0=A0-B1p(ip,jp,ip,jp)+B1p(ip,jp,3-ip,3-jp) temp=-dimag(A0) mue(1,3,ia)=mue(1,3,ia)+temp temp=-A0 mue(1,4,ia)=mue(1,4,ia)+temp temp=A0p(ip,jp,ip,jp)-A0p(ip,jp,3-ip,3-jp) mue(2,1,ia)=mue(2,1,ia)+temp temp=B0p(ip,jp,ip,jp)+B0p(ip,jp,3-ip,3-jp) mue(2,2,ia)=mue(2,2,ia)-temp A0=A1p(ip,jp,ip,jp)+A1p(ip,jp,3-ip,3-jp) A0=A0+B1p(ip,jp,ip,jp)-B1p(ip,jp,3-ip,3-jp) temp=-dimag(A0) mue(2,3,ia)=mue(2,3,ia)+temp temp=-A0 mue(2,4,ia)=mue(2,4,ia)+temp temp=0.5d0*dimag(A1p(ip,jp,ip,3-jp)) mue(3,1,ia)=mue(3,1,ia)+temp temp=-0.5d0*dimag(B1p(ip,jp,ip,3-jp)) mue(3,2,ia)=mue(3,2,ia)+temp A0=A0p(ip,jp,ip,3-jp)-A0p(ip,jp,3-ip,jp) A0=A0-B0p(ip,jp,ip,3-jp)-B0p(ip,jp,3-ip,jp) temp=A0 mue(3,3,ia)=mue(3,3,ia)+temp temp=-dimag(A0) mue(3,4,ia)=mue(3,4,ia)+temp temp=-0.5d0*A1p(ip,jp,ip,3-jp) mue(4,1,ia)=mue(4,1,ia)+temp temp=0.5d0*B1p(ip,jp,ip,3-jp) mue(4,2,ia)=mue(4,2,ia)+temp A0=A0p(ip,jp,ip,3-jp)+A0p(ip,jp,3-ip,jp) A0=A0-B0p(ip,jp,ip,3-jp)+B0p(ip,jp,3-ip,jp) temp=dimag(A0) mue(4,3,ia)=mue(4,3,ia)+temp temp=A0 mue(4,4,ia)=mue(4,4,ia)+temp if(ia.eq.iang) goto 19101 temp=A0pg(ip,jp,ip,jp)+A0pg(ip,jp,3-ip,3-jp) mue(1,1,iang)=mue(1,1,iang)+temp temp=B0pg(ip,jp,3-ip,3-jp)-B0pg(ip,jp,ip,jp) mue(1,2,iang)=mue(1,2,iang)+temp A0=A1pg(ip,jp,ip,jp)+A1pg(ip,jp,3-ip,3-jp) A0=A0-B1pg(ip,jp,ip,jp)+B1pg(ip,jp,3-ip,3-jp) temp=-dimag(A0) mue(1,3,iang)=mue(1,3,iang)+temp temp=-A0 mue(1,4,iang)=mue(1,4,iang)+temp temp=A0pg(ip,jp,ip,jp)-A0pg(ip,jp,3-ip,3-jp) mue(2,1,iang)=mue(2,1,iang)+temp temp=B0pg(ip,jp,ip,jp)+B0pg(ip,jp,3-ip,3-jp) mue(2,2,iang)=mue(2,2,iang)-temp A0=A1pg(ip,jp,ip,jp)+A1pg(ip,jp,3-ip,3-jp) A0=A0+B1pg(ip,jp,ip,jp)-B1pg(ip,jp,3-ip,3-jp) temp=-dimag(A0) mue(2,3,iang)=mue(2,3,iang)+temp temp=-A0 mue(2,4,iang)=mue(2,4,iang)+temp temp=0.5d0*dimag(A1pg(ip,jp,ip,3-jp)) mue(3,1,iang)=mue(3,1,iang)+temp temp=-0.5d0*dimag(B1pg(ip,jp,ip,3-jp)) mue(3,2,iang)=mue(3,2,iang)+temp A0=A0pg(ip,jp,ip,3-jp)-A0pg(ip,jp,3-ip,jp) A0=A0-B0pg(ip,jp,ip,3-jp)-B0pg(ip,jp,3-ip,jp) temp=A0 mue(3,3,iang)=mue(3,3,iang)+temp temp=-dimag(A0) mue(3,4,iang)=mue(3,4,iang)+temp temp=-0.5d0*A1pg(ip,jp,ip,3-jp) mue(4,1,iang)=mue(4,1,iang)+temp temp=0.5d0*B1pg(ip,jp,ip,3-jp) mue(4,2,iang)=mue(4,2,iang)+temp A0=A0pg(ip,jp,ip,3-jp)+A0pg(ip,jp,3-ip,jp) A0=A0-B0pg(ip,jp,ip,3-jp)+B0pg(ip,jp,3-ip,jp) temp=dimag(A0) mue(4,3,iang)=mue(4,3,iang)+temp temp=A0 mue(4,4,iang)=mue(4,4,iang)+temp 19101 continue enddo enddo temp=mue(1,1,ia)+mue(1,2,ia)+mue(2,1,ia)+mue(2,2,ia) i22(ia)=0.5d0*temp temp=mue(1,1,ia)-mue(1,2,ia)-mue(2,1,ia)+mue(2,2,ia) i11(ia)=0.5d0*temp temp=mue(1,1,ia)+mue(1,2,ia)-mue(2,1,ia)-mue(2,2,ia) i21(ia)=0.5d0*temp temp=mue(1,1,ia)-mue(1,2,ia)+mue(2,1,ia)-mue(2,2,ia) i12(ia)=0.5d0*temp if(ia.eq.iang) goto 1910 temp=mue(1,1,iang)+mue(1,2,iang)+mue(2,1,iang)+mue(2,2,iang) i22(iang)=0.5d0*temp temp=mue(1,1,iang)-mue(1,2,iang)-mue(2,1,iang)+mue(2,2,iang) i11(iang)=0.5d0*temp temp=mue(1,1,iang)+mue(1,2,iang)-mue(2,1,iang)-mue(2,2,iang) i21(iang)=0.5d0*temp temp=mue(1,1,iang)-mue(1,2,iang)+mue(2,1,iang)-mue(2,2,iang) i12(iang)=0.5d0*temp 1910 continue cbak=i11(2*nang-1) do i=1,nang2 inat(i)=i11(i)+i22(i)+i12(i)+i21(i) pol(i)=(i11(i)+i12(i)-i22(i)-i21(i))/inat(i) enddo 2000 cz=twopi/k**2 csca=csca*cz cext=cext*cz cabs=cext-csca cpr=cpr*cz assym=cpr/csca cbak=2.d0*cbak*cz write(6,'(5x,a6,7x,a6,7x,a6,7x,a6,7x,a5,6x,a12)') +'','','','','','' write(6,'(2x,6e13.5)') cext,cabs,csca,cbak,cext-cpr,assym cscax=0.d0 cextx=0.d0 cprx=0.d0 do i=1,nL cscax=cscax+cscai(i) cextx=cextx+cexti(i) cprx=cprx+cpri(i) enddo cscax=cscax*cz cextx=cextx*cz cprx=cprx*cz cabsx=cextx-cscax assymx=cprx/cscax cbakx=cbak write(6,'(2x,6e13.5)') cextx,cabsx,cscax,cbakx,cextx-cprx,assymx do i=1,nL cexti(i)=cexti(i)*cz cscai(i)=cscai(i)*cz cabsi(i)=cexti(i)-cscai(i) cpri(i)=cpri(i)*cz assymi(i)=cpri(i)/csca write(6,'(i5,5e15.6)') i,cexti(i),cabsi(i), + cscai(i),cpri(i),assymi(i) enddo c flout='cr'//fileout c open(33,file=flout,status='unknown') c write(33,'(a20,a47)') c + flout,'(Total and individual-particle cross sections)' c write(33,'(a32,2x,a22)') c + 'input sphere-aggregate filename:',FLNAME c write(33,'(12x,a4,11x,a4,11x,a4,11x,a3,8x,a12)') c + 'Cext','Cabs','Csca','Cpr','' c write(33,'(a5,5e15.6)') 'total',cext,cabs,csca,cext-cpr, c + assym c do i=1,nL c write(33,'(i5,5e15.6)') i,cexti(i),cabsi(i),cscai(i), c + cpri(i),assymi(i) c enddo c close(33) cz=pione*gcvr*gcvr cscav=csca/cz cextv=cext/cz cprv=cpr/cz cbakv=cbak/cz cabsv=cextv-cscav temp=gcvr*gcvr/gcs cscas=cscav*temp cexts=cextv*temp cabss=cabsv*temp cbaks=cbakv*temp cprs=cprv*temp 222 format(6e13.5) 221 format(3x,a7,6x,a7,6x,a7,6x,a7,6x,a6,5x,a12) write(6,221) +'','','','','','' write(6,222) cextv,cabsv,cscav,cbakv,cprv,assym write(6,221) +'','','','','','' write(6,222) cexts,cabss,cscas,cbaks,cprs,assym open(12,file=fileout,status='unknown') write(12,'(/)') write(12,'(1x,a15,a13,a18,a4,f8.3,a5,f8.3)') fileout, + ' input file: ',FLNAME,' xv:',xv,' xs:',xs if(idscmt.lt.0) then write(12,'(/)') write(12,*) '*** backscattering and scatttering matrix', + ' are not calculated ***' endif write(12,'(/)') write(12,221) +'','','','','','' write(12,222) + cext,cabs,csca,cbak,cext-cpr,assym write(12,221) +'','','','','','' write(12,222) cextv,cabsv,cscav,cbakv,cprv,assym write(12,221) +'','','','','','' write(12,222) cexts,cabss,cscas,cbaks,cprs,assym if(idscmt.lt.0) goto 2001 write(12,'(/)') write(12,'(2x,a4,4x,a7,4x,a6,4x,a7,6x,a7,6x,a7,6x,a7)') +'s.a.','','','','','','' do i=1,nang2 write(12,'(f6.1,e13.5,f8.4,4e13.5)') + dang(i),inat(i),pol(i),i11(i),i21(i),i12(i),i22(i) enddo write(12,'(/)') write(12,'(1x,a50)') + 'Scattering matrix (4X4 for each scattering angle):' do i=1,nang2 write(12,'(f7.1,4e16.7)') + dang(i),mue(1,1,i),mue(1,2,i),mue(1,3,i),mue(1,4,i) write(12,'(7x,4e16.7)') + mue(2,1,i),mue(2,2,i),mue(2,3,i),mue(2,4,i) write(12,'(7x,4e16.7)') + mue(3,1,i),mue(3,2,i),mue(3,3,i),mue(3,4,i) write(12,'(7x,4e16.7)') + mue(4,1,i),mue(4,2,i),mue(4,3,i),mue(4,4,i) enddo 2001 close(12) STOP END SUBROUTINE abMiexud(X,REFREL,NP,NMAX,NM,AN,BN,NADD, + RSR,RSI,RSX,PX,AR,AI,BR,BI,EPS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer NMAX,NM,NADD,NSTOP,NX,K,N COMPLEX*16 REFREL,AN(NP),BN(NP) DOUBLE PRECISION AR(NP),AI(NP),BR(NP),BI(NP), + RSR(NMAX),RSI(NMAX),RSX(NMAX),PX(NMAX) IF(EPS.GT.1.D0.OR.EPS.LT.0.D0) EPS=1.0D-20 IF(NADD.NE.0) EPS=0.D0 CTC=EPS XM=DBLE(REFREL) YM=DIMAG(REFREL) XMX=X*XM YMX=X*YM RP2=XMX*XMX+YMX*YMX NSTOP=X+4.D0*X**.3333D0 NSTOP=NSTOP+2 NM=NSTOP+NADD XN=DSQRT(XM**2+YM**2)*X NX=1.1D0*XN+10.D0 if(NX-NM.lt.10) NX=NM+10 write(6,*) 'Wiscombe criterion:',NSTOP write(6,*) 'NADD:',NADD write(6,*) 'NX:',NX IF(NX.GT.NMAX) THEN WRITE(6,*) 'Parameter NXMAX too small' WRITE(6,*) ' NXMAX must be greater than', NX WRITE(6,*) 'Please correct NXMAX in main code,' WRITE(6,*) ' recompile, then try again' STOP ENDIF IF(NM.GT.NP) THEN WRITE(6,*) 'Parameter np too small' WRITE(6,*) ' np must be greater than', NM WRITE(6,*) 'Please correct np in gmm01f.par,' WRITE(6,*) ' recompile the code, then try again' STOP ENDIF C DOWN RECURSION FOR RATIOS RSR,RSI,RSX,PNR,PNI,PX PNX=X/DBLE(2*NX+3) PNR=XMX/DBLE(2*NX+3) PNI=YMX/DBLE(2*NX+3) DO 5 K=1,NX N=NX-K+1 CN=DBLE(N) ALN=(2.D0*CN+1.D0)*XMX/RP2-PNR BEN=(2.D0*CN+1.D0)*YMX/RP2+PNI RSR(N)=-CN*XMX/RP2+ALN RSI(N)=CN*YMX/RP2-BEN PZD=ALN*ALN+BEN*BEN PNR=ALN/PZD PNI=BEN/PZD RSX(N)=(CN+1.D0)/X-PNX IF(N.EQ.1) GO TO 20 PNX=X/(2.D0*CN+1.D0-PNX*X) PX(N)=PNX 5 CONTINUE 20 SNM1X=DSIN(X) CNM1X=DCOS(X) IF(X-0.1D0)21,22,22 21 SNX=X**2./3.D0-X**4./30.D0+X**6./840.D0-X**8./45360.D0 GO TO 23 22 SNX=SNM1X/X-CNM1X 23 CNX=CNM1X/X+SNM1X DO 10 N=1,NX PX(N)=SNX !preparing for the calculation of Cabs C=DBLE(N) DCNX=CNM1X-C*CNX/X DSNX=RSX(N)*SNX C CALCULATION OF EXTERIOR COEFFICIENTS AN AND BN ANNR=RSR(N)*SNX-XM*DSNX ANNI=RSI(N)*SNX-YM*DSNX TA1=RSR(N)*SNX-RSI(N)*CNX TA2=RSI(N)*SNX+RSR(N)*CNX ANDR=TA1-XM*DSNX+YM*DCNX ANDI=TA2-XM*DCNX-YM*DSNX AND=ANDR*ANDR+ANDI*ANDI BNNR=(XM*RSR(N)-YM*RSI(N))*SNX-DSNX BNNI=(XM*RSI(N)+YM*RSR(N))*SNX TB1=RSR(N)*SNX-RSI(N)*CNX TB2=RSR(N)*CNX+RSI(N)*SNX BNDR=XM*TB1-YM*TB2-DSNX BNDI=XM*TB2+YM*TB1-DCNX BND=BNDR*BNDR+BNDI*BNDI AR(N)=(ANNR*ANDR+ANNI*ANDI)/AND AI(N)=(ANNI*ANDR-ANNR*ANDI)/AND BR(N)=(BNNR*BNDR+BNNI*BNDI)/BND BI(N)=(BNNI*BNDR-BNNR*BNDI)/BND C MIE SERIES CONVERGENCE TEST IS MADE BY TESTING AN'S AND BN'S TI=AR(N)*AR(N)+AI(N)*AI(N)+BR(N)*BR(N)+BI(N)*BI(N) TI=TI/(AR(1)*AR(1)+AI(1)*AI(1)+BR(1)*BR(1)+BI(1)*BI(1)) IF(TI-CTC) 16,18,18 18 IF(NM-N) 15,15,6 6 IF(N-NX)7,15,15 7 M=N+1 SNX=PX(M)*SNX CNM2X=CNM1X CNM1X=CNX CNX=(2.D0*C+1.D0)*CNM1X/X-CNM2X 10 CONTINUE GO TO 15 16 WRITE(6,*) '*** NOTE THAT THE FIELD-EXPANSION TRANCATION' WRITE(6,*) '*** IS DETERMINED BY eps GIVEN IN THE INPUT' WRITE(6,*) '*** FILE gmm01f.in' WRITE(6,*) '*** IN CASE YOU NEED A HIGHER ORDER, eps MUST' WRITE(6,'(a,e9.1)') + ' *** BE SMALLER THAN THE CURRENT VALUE',EPS 15 NM=N DO I=1,NM AN(I)=DCMPLX(AR(I),-AI(I)) BN(I)=DCMPLX(BR(I),-BI(I)) ENDDO RETURN END subroutine cofsrd(nmax) include 'gmm01f.par' parameter (nmp=np*(np+2)) double precision cofsr(nmp),lnfacd,c common/crot/cofsr i=0 do n=1,nmax do m=-n,n i=i+1 c=lnfacd(dble(n-m))-lnfacd(dble(n+m)) cofsr(i)=0.5d0*c c c=0.5d0*c c cofsr(i)=dexp(c) enddo enddo return end subroutine cofd0(nmax) implicit double precision (a-h,o-z) include 'gmm01f.par' parameter (nmp=np*(np+2)) parameter (ni0=np*(np+1)*(2*np+1)/3+np*np) integer v double precision lnfacd common/cofmnv0/cof0(ni0) common/crot/cofsr(nmp) common/fnr/fnr(0:4*(np+1)) i=0 sm=-0.5d0*dble((-1)**nmax) do m=-nmax,nmax ns=max(1,iabs(m)) sm=-sm do n=ns,nmax inm=n*(n+1)-m do v=ns,nmax i=i+1 ivm=v*(v+1)+m c=cofsr(inm)+cofsr(ivm) c=sm*dexp(c) c0=fnr(2*n+1)*fnr(2*v+1) c1=fnr(n)*fnr(v)*fnr(n+1)*fnr(v+1) c0=c0/c1 cof0(i)=c*c0 enddo enddo enddo return end subroutine cofnv0(nmax) include 'gmm01f.par' integer n,v double precision c1,lnfacd,cnv(np,np) common/cfnv/cnv do n=1,nmax do v=n,nmax c1=lnfacd(dble(2*n))+lnfacd(dble(2*v)) c1=c1-lnfacd(dble(2*n+2*v)) c1=c1+2.d0*lnfacd(dble(n+v)) c1=c1-lnfacd(dble(n))-lnfacd(dble(v)) cnv(n,v)=c1 enddo enddo return end C subroutine gau0.f generates tabulated values for C Gaunt coefficients up to n=v=n_max subroutine gau0(nmax) include 'gmm01f.par' parameter (ni0=np*(np+1)*(2*np+1)/3+np*np) parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6) integer v,qmax,uvmax,iga0(ni0) double precision ga0(ng0) common/g0/ga0 common/ig0/iga0 na=0 uvmax=nmax*(nmax+2) i=0 do m=-nmax,nmax ns=max(1,iabs(m)) do n=ns,nmax do v=ns,nmax call gxurcd0(-m,n,v,qmax,na) i=i+1 iga0(i)=na na=na+qmax+1 enddo enddo enddo return end c transforms the rectangular coordinates (x,y,z) c to spherical coordinates (r,theta,phi) subroutine carsphd(x,y,z,r,xt,sphi,cphi) double precision x,y,z,r,xt,sphi,cphi r=dsqrt(x*x+y*y+z*z) if(r.eq.0.d0) then xt=1.d0 sphi=0.d0 cphi=1.d0 return endif xt=z/r if(y.eq.0.d0.and.x.eq.0.d0) then sphi=0.d0 cphi=1.d0 return endif sphi=dsqrt(x*x+y*y) cphi=x/sphi sphi=y/sphi return end c subroutine besseljd.f (in double precision arithmetic) c returns an array of the spherical Bessel function of the c first kind with a real argument: j_0,j_1,j_2,...,j_{NC} c uses Ru Wang's ratio method for the downward recursive c calculation of the Riccati-Bessel function Psi_n(z)=z j_n(z) c [see Xu et al., Physical Review E, v.60, 2347-2365 (1999)] SUBROUTINE besseljd(NC,X,BESJ) INTEGER NC,NX,K,N DOUBLE PRECISION X,BESJ(0:NC),PN,CN,X2 DO K=1,NC BESJ(K)=0.D0 ENDDO IF(DABS(X).LT.1.D-12) THEN BESJ(0)=1.D0 BESJ(1)=1.D0/3.D0*X RETURN ENDIF c down-recursively calculating an array of the ratio functions c P_{NC},P_{NC-1},...,P(2) stored in the same array for j_n, c starting with an asymptotic value P_{NX+1}=X/(2 NX+3) at the c highest order NX+1, where NX=NC+1.1X+10 NX=1.1D0*X+10.D0 NX=NC+NX PN=X/DBLE(2*NX+3) DO 5 K=1,NX-1 N=NX-K+1 CN=DBLE(N) PN=X/(DBLE(2*N+1)-PN*X) IF(N.GT.NC) GOTO 5 BESJ(N)=PN 5 CONTINUE C calculating j_0(x) and j_1(x) IF(DABS(X)-0.1D0) 10,11,11 10 X2=X*X BESJ(0)=1.D0-X2/72.D0 BESJ(0)=1.D0-X2/42.D0*BESJ(0) BESJ(0)=1.D0-X2/20.D0*BESJ(0) BESJ(0)=1.D0-X2/6.D0*BESJ(0) BESJ(1)=1.D0/45360.D0-X2/3991680.D0 BESJ(1)=1.D0/840.D0-X2*BESJ(1) BESJ(1)=1.D0/30.0d0-X2*BESJ(1) BESJ(1)=X*(1.D0/3.0d0-X2*BESJ(1)) GOTO 12 11 BESJ(0)=DSIN(X)/X BESJ(1)=(DSIN(X)/X-DCOS(X))/X c calculating j_2,...,j_{NC} 12 DO 20 N=2,NC BESJ(N)=BESJ(N)*BESJ(N-1) 20 CONTINUE RETURN END c sub. besselyd.f (in double precision arithmetic) c returns an array of the spherical Bessel function of c the second kind with a real argument: y_0,y_1,...,y_n subroutine besselyd(n,x,besy) integer i,n double precision x,besy(0:n),besyn,x2 if(x.eq.0.d0) then write(6,*) 'bad argument in sub. besselyd' stop endif if(dabs(x)-0.1d0)10,11,11 10 x2=x*x besyn=1.d0-x2/72.d0 besyn=1.d0-x2/42.d0*besyn besyn=1.d0-x2/20.d0*besyn besyn=1.d0-x2/6.d0*besyn besy(0)=1.d0-x2/56.d0 besy(0)=1.d0-x2/30.d0*besy(0) besy(0)=1.d0-x2/12.d0*besy(0) besy(0)=1.d0-0.5d0*x2*besy(0) besy(0)=-besy(0)/x goto 12 11 besyn=dsin(x)/x besy(0)=-dcos(x)/x 12 besy(1)=besy(0)/x-besyn do i=2,n besy(i)=dble(2*i-1)/x*besy(i-1)-besy(i-2) enddo return end c subroutine rotcoef is originally written by Mackowski, Fuller, and c Mishchenko (taken from the code scsmfo1b.for released to public by c the authors at 8/1999) c the rotational coefficients are required for the implementation of c Mackowski's three-step numerical technique in subroutine rtr.f for c decomposition of translation matrix into rotational and axial c translational parts [see Mackowski, Proc. R. Soc. Lond. A 433, 599 c (1991)] c Yu-lin Xu 12/2000 subroutine rotcoef(cbe,nmax) include 'gmm01f.par' parameter (nmp=np*(np+2)) implicit double precision (a-h,o-z) double precision dk0(-2*np:2*np),dk01(-2*np:2*np) common/rot/bcof(0:np+2),dc(-np:np,0:nmp) common/fnr/fnr(0:4*(np+1)) sbe=dsqrt((1.d0+cbe)*(1.d0-cbe)) cbe2=.5d0*(1.d0+cbe) sbe2=.5d0*(1.d0-cbe) in=1 dk0(0)=1.d0 sben=1.d0 dc(0,0)=1.d0 dk01(0)=0.d0 do n=1,nmax nn1=n*(n+1) in=-in sben=sben*sbe/2.d0 dk0(n)=dble(in)*sben*bcof(n) dk0(-n)=dble(in)*dk0(n) dk01(n)=0.d0 dk01(-n)=0.d0 dc(0,nn1+n)=dk0(n) dc(0,nn1-n)=dk0(-n) do k=-n+1,n-1 kn=nn1+k dkt=dk01(k) dk01(k)=dk0(k) dk0(k)=(cbe*dble(n+n-1)*dk01(k)-fnr(n-k-1)*fnr(n+k-1)*dkt) 1 /(fnr(n+k)*fnr(n-k)) dc(0,kn)=dk0(k) enddo im=1 do m=1,n im=-im fmn=1.d0/fnr(n-m+1)/fnr(n+m) m1=m-1 dkm0=0.d0 do k=-n,n kn=nn1+k dkm1=dkm0 dkm0=dc(m1,kn) if(k.eq.n) then dkn1=0.d0 else dkn1=dc(m1,kn+1) endif dc(m,kn)=(fnr(n+k)*fnr(n-k+1)*cbe2*dkm1 1 -fnr(n-k)*fnr(n+k+1)*sbe2*dkn1 1 -dble(k)*sbe*dc(m1,kn))*fmn dc(-m,nn1-k)=dble((-1)**(k)*im)*dc(m,kn) enddo enddo enddo return end c subroutine cofxuds0.f returns the two classes of vector c (axial) translation coefficients for a given combination of c (m,n,m,v) and a given dimensionless translation distance kd cu uses subroutine gid0.f subroutine cofxuds0(nmax,m,n,v,sja,sya,A,B,Aj,Bj) implicit double precision (a-h,o-z) include 'gmm01f.par' parameter (ni0=np*(np+1)*(2*np+1)/3+np*np) parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6) integer v,p,qmax double precision sja(0:n+v+1),sya(0:n+v+1) complex*16 A,B,Aj,Bj,signz common/ig0/iga0(ni0) common/g0/ga0(ng0) common/cofmnv0/cof0(ni0) fa(m,p)=dble(-2*m*p*(p-1)) fb(n,v,p)=dble(p*p-(n+v+1)*(n+v+1)) + *dble(p*p-(n-v)*(n-v))/dble(4*p*p-1) if(iabs(m).gt.n.or.iabs(m).gt.v) then write(6,*) '|m|>n or v in subroutine cofxuds0.f' stop endif A=0.d0 B=0.d0 Aj=0.d0 Bj=0.d0 call gid0(nmax,m,n,v,id) c=cof0(id) ig=iga0(id) nv2=v*(v+1)+n*(n+1) signz=dcmplx(0.d0,1.d0)**(n+v) p=n+v+2 qmax=min(n,v) do i=1,qmax+1 p=p-2 cp=dble(nv2-p*(p+1))*ga0(ig+i) sj=sja(p) sy=sya(p) A=A+dcmplx(sj,sy)*signz*cp Aj=Aj+sj*signz*cp signz=-signz enddo A=A*c Aj=Aj*c if(m.eq.0) return signz=dcmplx(0.d0,1.d0)**(n+v+1) p=n+v do 20 i=1,qmax p=p-2 signz=-signz if(i.eq.1) then cp=dble(2*p+3)*fa(m,p+3) cp=cp*ga0(ig+1)/dble((p+3)*(p+2)) goto 21 endif if(i.eq.qmax) then if(p.eq.0) goto 22 nv2=p*(p+1) cp=dble(2*p+3)*fa(m,p+1) cp=-cp*ga0(ig+i+1)/dble(nv2) goto 21 endif 22 c4=fa(m,p+2) cp=-dble((p+1)*(p+2))*fb(n,v,p+2)*ga0(ig+i) cp=cp+dble((p+2)*(p+1))*fb(n,v,p+1)*ga0(ig+i+1) cp=cp*dble(2*p+3)/c4 21 sj=sja(p+1) sy=sya(p+1) B=B+dcmplx(sj,sy)*signz*cp Bj=Bj+sj*signz*cp 20 continue B=B*c Bj=Bj*c return end c subroutine tipitaud.f c calculates pi(m,n) & tau(m,n) up to a specified nmax for all c m=0,1,...n at a given x c pi(m,n) & tau(m,n) calculated are normalized by c C_mn=[(2n+1)(n-m)!/n/(n+1)/(n+m)!]^(1/2) c Yu-lin Xu 12/2000 subroutine tipitaud(nmax,x) include 'gmm01f.par' parameter (nmp0=(np+1)*(np+4)/2) implicit double precision (a-h,o-z) common/fnr/fnr(0:4*(np+1)) common/pitau/pi(nmp0),tau(nmp0) nt=(nmax+1)*(nmax+4)/2 ! calculates pi up to nmax+1 if(nt.gt.nmp0.or.dabs(x).gt.1.d0) then write(6,*) 'dimension or argument wrong in sub. tipitaud' write(6,*) 'argument: ',x stop endif sx=dsqrt(1.d0-x*x) pi(1)=0.d0 ! pi(0,1) pi(0,n)=0 when m=0 pi(2)=dsqrt(.75d0) ! pi(1,1) pi(3)=0.d0 ! pi(0,2) t125=dsqrt(1.25d0) pi(4)=t125*x ! pi(1,2) pi(5)=t125*sx ! pi(2,2) imn=5 do i=3,nmax+1 n=i imn=imn+1 pi(imn)=0.d0 ! pi(0,n)=0 do 11 j=2,n m=j-1 imn=imn+1 i1=(n-2)*(n+1)/2+m+1 if(m.eq.n-1) then pi(imn)=fnr(n-1)*fnr(2*n+1)/fnr(n+1)*x*pi(i1) goto 11 endif i2=(n-3)*n/2+m+1 t=fnr(n)*fnr(2*n-3) t=fnr(n-2)*fnr(n-m-1)*fnr(n+m-1)/t pi(imn)=fnr(2*n-1)*x*pi(i1)-t*pi(i2) t=fnr(n+1)*fnr(n-m)*fnr(n+m) t=fnr(n-1)*fnr(2*n+1)/t pi(imn)=t*pi(imn) 11 continue imn=imn+1 i1=(n-2)*(n+1)/2+n t=fnr(n-1)*fnr(n+1) t=dsqrt(.5d0)*fnr(n)*fnr(2*n+1)/t pi(imn)=t*sx*pi(i1) enddo tx=x*sx tau(1)=-dsqrt(1.5d0)*sx ! tau(0,1) tau(2)=pi(2)*x ! tau(1,1) tau(3)=-dsqrt(7.5d0)*tx ! tau(0,2) tau(4)=t125*(2.d0*x*x-1.d0) ! tau(1,2) tau(5)=t125*tx ! tau(2,2) imn=5 do i=3,nmax n=i do 21 j=1,n+1 m=j-1 imn=imn+1 if(m.eq.0) then i1=(n-2)*(n+1)/2+1 i2=(n-3)*n/2+1 t=fnr(2*n-3) t=fnr(n-2)*fnr(n)/t tau(imn)=fnr(2*n-1)*x*tau(i1)-t*tau(i2) t=fnr(n-1)*fnr(n+1) t=fnr(2*n+1)/t tau(imn)=t*tau(imn) goto 21 endif i1=n*(n+3)/2+m+1 t=fnr(n)*fnr(2*n+3) t=fnr(n+2)*fnr(2*n+1)*fnr(n-m+1)*fnr(n+m+1)/t tau(imn)=t*pi(i1)-dble(n+1)*x*pi(imn) tau(imn)=tau(imn)/dble(m) 21 continue enddo return end c lnfacd.f (double precision function) c returns ln(z!) z>-1.0 c based on Lanczos' method [see Xu, Journal of Computational c Physics, v.139, 137-165 (1998)] double precision function lnfacd(z) integer i double precision z,a,b,cp,c0(11) data c0/0.16427423239836267d5, -0.48589401600331902d5, + 0.55557391003815523d5, -0.30964901015912058d5, + 0.87287202992571788d4, -0.11714474574532352d4, + 0.63103078123601037d2, -0.93060589791758878d0, + 0.13919002438227877d-2,-0.45006835613027859d-8, + 0.13069587914063262d-9/ a=1.d0 cp=2.5066282746310005d0 b=z+10.5d0 b=(z+0.5d0)*dlog(b)-b do i=1,11 z=z+1.d0 a=a+c0(i)/z enddo lnfacd=b+dlog(cp*a) return end c gxurcd0.f to compute Gaunt coefficients a(-m,n,m,v,p) cu uses lnfacd.f to compute ln(z!) subroutine gxurcd0(m,n,v,qmax,na) implicit double precision (a-h,o-z) include 'gmm01f.par' parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6) integer v,qmax,p double precision lnfacd,cnv(np,np),ga0(ng0) common/cfnv/cnv common/g0/ga0 fb(n,v,p)=dble(p-(n+v+1))*dble(p+(n+v+1)) + *dble(p-(n-v))*dble(p+(n-v)) + /(dble(2*p+1)*dble(2*p-1)) if(iabs(m).gt.n.or.iabs(m).gt.v) then write(6,*) 'warning: |m|>n or v in gxurcd0' qmax=-1 return endif qmax=min(n,v) nq=qmax+1 if(n.le.v) then c1=cnv(n,v) else c1=cnv(v,n) endif c1=c1-lnfacd(dble(n-m))-lnfacd(dble(v+m)) ga0(na+1)=dexp(c1) if(qmax.lt.1) return p=n+v do 8 i=2,nq p=p-2 if(m.eq.0) then c1=fb(n,v,p+1) c2=fb(n,v,p+2) goto 2 endif c1=fb(n,v,p+1) c2=dble(4*m*m)+fb(n,v,p+2)+fb(n,v,p+3) if(i.eq.2) goto 2 c3=-fb(n,v,p+4) goto 4 2 ga0(na+i)=c2*ga0(na+i-1)/c1 goto 8 4 ga0(na+i)=(c2*ga0(na+i-1)+c3*ga0(na+i-2))/c1 8 continue return end subroutine gid0(nmax,m,n,iv,id) nt=nmax*(nmax+1)*(2*nmax+1)/3+nmax*nmax ns=max(1,iabs(m)) nc0=nmax-iabs(m) id=nc0*(nc0+1)*(2*nc0+1)/6 if(m) 10,11,12 10 id=id+(n-ns)*(nc0+1)+iv-ns+1 return 11 id=id+(n-ns)*nmax+iv return 12 id=id+(nmax-n)*(nc0+1)+nmax-iv id=nt-id return end function plgndrd(l,mr,x) integer l,mr,m,index,i,ll double precision x,fact,pll,pmm,pmmp1,somx2 double precision plgndrd,lnfacd m=mr index=0 if(m.lt.0) then index=1 m=-m endif if(dabs(x).gt.1.d0) then write(6,*) 'n,m,x: ', l,(-2*index+1)*m,x pause 'bad arguments in plgndrd' end if if(m.gt.l) then plgndrd=0.d0 return end if pmm=1.d0 if(m.gt.0) then somx2=dsqrt((1.d0-x)*(1.d0+x)) fact=1.d0 do 11 i=1,m pmm=-pmm*fact*somx2 fact=fact+2.d0 11 continue endif if(l.eq.m) then plgndrd=pmm else pmmp1=x*(2*m+1)*pmm if(l.eq.m+1) then plgndrd=pmmp1 else do 12 ll=m+2,l pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m) pmm=pmmp1 pmmp1=pll 12 continue plgndrd=pll endif endif plgndrd=-plgndrd if(m/2*2.eq.m) plgndrd=-plgndrd if(index.gt.0) then fact=lnfacd(dble(l-m))-lnfacd(dble(l+m)) fact=dexp(fact) plgndrd=-plgndrd*fact if(m/2*2.eq.m) plgndrd=-plgndrd endif return end C subroutine rtrT(anpt,nodrj,nodri,ekt,drot,ij1,ij2,ii1,ii2) implicit double precision (a-h,o-z) include 'gmm01f.par' parameter(nmp=np*(np+2)) parameter (nrc=4*np*(np+1)*(np+2)/3+np,nij=nLp*(nLp-1)/2) double precision drot(nrc) complex*16 anpt(2,nmp),ant(2,2*np),amt(2,-np:np),a,b, + ekt(np),ek(-np:np),atr(2,np,nmp) common/tran/atr nj1=dsqrt(dble(ij1)) nj2=dsqrt(dble(ij2)) ni1=dsqrt(dble(ii1)) ni2=dsqrt(dble(ii2)) ek(0)=1.d0 nmax=max(nj2,ni2) do m=1,nmax ek(m)=ekt(m) ek(-m)=dconjg(ek(m)) enddo irc=0 do n=1,nj1-1 do k=-n,n do m=-n,n irc=irc+1 enddo enddo enddo do n=nj1,nj2 n1=n*(n+1) do m=-n,n amt(1,m)=0.d0 amt(2,m)=0.d0 enddo do 10 k=-n,n kn=n1+k if(kn.lt.ij1) then irc=irc+2*n+1 goto 10 endif if(kn.gt.ij2) goto 11 a=ek(k)*anpt(1,kn) b=ek(k)*anpt(2,kn) do m=-n,n irc=irc+1 amt(1,m)=amt(1,m)+a*drot(irc) amt(2,m)=amt(2,m)+b*drot(irc) enddo 10 continue 11 do m=-n,n imn=n1+m anpt(1,imn)=amt(1,m) anpt(2,imn)=amt(2,m) enddo enddo mmax=min(nj2,ni2) do m=-mmax,mmax n1=max(1,iabs(m)) n1j=max(n1,nj1) do n=n1j,nj2 imn=n*(n+1)+m do ip=1,2 ant(ip,n)=anpt(ip,imn) enddo enddo n1i=max(n1,ni1) do n=n1i,ni2 imn=n*(n+1)+m a=0.d0 b=0.d0 do 20 l=n1j,nj2 ml=l*(l+1)+m a=a+atr(1,n,ml)*ant(1,l) 1 +atr(2,n,ml)*ant(2,l) b=b+atr(1,n,ml)*ant(2,l) 1 +atr(2,n,ml)*ant(1,l) 20 continue anpt(1,imn) = a anpt(2,imn) = b enddo enddo in=1 irc=0 do n=1,ni1-1 do k=-n,n do m=-n,n irc=irc+1 enddo enddo enddo do n=ni1,ni2 in=-in n1=n*(n+1) do m=-n,n amt(1,m)=0.d0 amt(2,m)=0.d0 enddo sik=-in do 30 k=-n,n sik=-sik kn=n1+k a=sik*anpt(1,kn) b=sik*anpt(2,kn) do m=-n,n irc=irc+1 amt(1,m)=amt(1,m)+a*drot(irc) amt(2,m)=amt(2,m)+b*drot(irc) enddo 30 continue sik=-in do 31 m=-n,n sik=-sik imn=n1+m if(imn.lt.ii1.or.imn.gt.ii2) goto 31 anpt(1,imn)=amt(1,m)*ek(-m)*sik anpt(2,imn)=amt(2,m)*ek(-m)*sik 31 continue enddo return end subroutine transT(nL,r0,nmax,uvmax,fint,atr0,btr0,ek, + drot,as,bs,as1,bs1,ind,confg,iuvc,isw) implicit double precision (a-h,o-z) include 'gmm01f.par' parameter (nmp=np*(np+2)) parameter (ni0=np*(np+1)*(2*np+1)/3+np*np) parameter (nrc=4*np*(np+1)*(np+2)/3+np) parameter (nij=nLp*(nLp-1)/2) integer v,nmax(nLp),uvmax(nLp),ind(nLp) double precision r0(6,nLp),drot(nrc,nij),confg(5,nij) complex*16 atr(2,np,nmp),atr0(ni0,nij),btr0(ni0,nij), + ek(np,nij),as(nLp,nmp),bs(nLp,nmp),as1(nLp,nmp), + bs1(nLp,nmp),at1(2,nmp) common/tran/atr do i=1,nL do imn=1,uvmax(i) as1(i,imn)=dcmplx(0.d0,0.d0) bs1(i,imn)=dcmplx(0.d0,0.d0) enddo enddo do 11 i=1,nL-1 if(ind(i).gt.0) goto 11 do 10 j=i+1,nL if(ind(j).gt.0) goto 10 ij=(j-1)*(j-2)/2+j-i x0=confg(1,ij) y0=confg(2,ij) z0=confg(3,ij) temp=confg(5,ij) if(temp.le.fint) goto 10 nlarge=max(nmax(i),nmax(j)) itrc=0 nsmall=min(nmax(i),nmax(j)) do m=-nsmall,nsmall n1=max(1,iabs(m)) do n=n1,nlarge do v=n1,nlarge itrc=itrc+1 iuv=v*(v+1)+m atr(1,n,iuv)=atr0(itrc,ij) atr(2,n,iuv)=btr0(itrc,ij) if(x0.eq.0.d0.and.y0.eq.0.d0) then if(z0.lt.0.d0) goto 20 endif goto 21 20 sic=dble((-1)**(n+v)) atr(1,n,iuv)=sic*atr(1,n,iuv) atr(2,n,iuv)=-sic*atr(2,n,iuv) 21 continue enddo enddo enddo do iuv=1,uvmax(j) at1(1,iuv)=as(j,iuv) at1(2,iuv)=bs(j,iuv) enddo If(nmax(j).lt.nlarge) then do n=nmax(j)+1,nlarge do m=-n,n iuv=n*(n+1)+m at1(1,iuv)=dcmplx(0.d0,0.d0) at1(2,iuv)=dcmplx(0.d0,0.d0) enddo enddo endif goto(22,23),isw 22 ij1=iuvc ij2=uvmax(j) ii1=iuvc ii2=uvmax(i) goto 24 23 ij1=1 ij2=iuvc-1 ii1=iuvc ii2=uvmax(i) 24 if(x0.eq.0.d0.and.y0.eq.0.d0) then call trvT(at1,nmax(j),nmax(i),ij1,ij2,ii1,ii2) else call rtrT(at1,nmax(j),nmax(i),ek(1,ij), + drot(1,ij),ij1,ij2,ii1,ii2) endif do imn=ii1,ii2 as1(i,imn)=as1(i,imn)+at1(1,imn) bs1(i,imn)=bs1(i,imn)+at1(2,imn) enddo do m=-nsmall,nsmall n1=max(1,iabs(m)) do n=n1,nlarge do v=n1,nlarge iuv=v*(v+1)+m sic=dble((-1)**(n+v)) atr(1,n,iuv)=sic*atr(1,n,iuv) atr(2,n,iuv)=-sic*atr(2,n,iuv) enddo enddo enddo do iuv=1,uvmax(i) at1(1,iuv)=as(i,iuv) at1(2,iuv)=bs(i,iuv) enddo If(nmax(i).lt.nlarge) then do n=nmax(i)+1,nlarge do m=-n,n iuv=n*(n+1)+m at1(1,iuv)=dcmplx(0.d0,0.d0) at1(2,iuv)=dcmplx(0.d0,0.d0) enddo enddo endif goto(25,26),isw 25 ii1=iuvc ii2=uvmax(i) ij1=iuvc ij2=uvmax(j) goto 27 26 ii1=1 ii2=iuvc-1 ij1=iuvc ij2=uvmax(j) 27 if(x0.eq.0.d0.and.y0.eq.0.d0) then call trvT(at1,nmax(i),nmax(j),ii1,ii2,ij1,ij2) else call rtrT(at1,nmax(i),nmax(j),ek(1,ij), + drot(1,ij),ii1,ii2,ij1,ij2) endif do imn=ij1,ij2 as1(j,imn)=as1(j,imn)+at1(1,imn) bs1(j,imn)=bs1(j,imn)+at1(2,imn) enddo 10 continue 11 continue return end subroutine trvT(anpt,nodrj,nodri,ij1,ij2,ii1,ii2) implicit double precision (a-h,o-z) include 'gmm01f.par' parameter(nmp=np*(np+2)) complex*16 anpt(2,nmp),ant(2,2*np),a,b,atr(2,np,nmp) common/tran/atr nj1=dsqrt(dble(ij1)) nj2=dsqrt(dble(ij2)) ni1=dsqrt(dble(ii1)) ni2=dsqrt(dble(ii2)) mmax=min(nj2,ni2) do m=-mmax,mmax n1=max(1,iabs(m)) n1j=max(n1,nj1) do 10 n=n1j,nj2 imn=n*(n+1)+m if(imn.lt.ij1.or.imn.gt.ij2) goto 10 do ip=1,2 ant(ip,n)=anpt(ip,imn) enddo 10 continue n1i=max(n1,ni1) do 20 n=n1i,ni2 imn=n*(n+1)+m if(imn.lt.ii1.or.imn.gt.ii2) goto 20 a=0.d0 b=0.d0 do 21 l=n1j,nodrj ml=l*(l+1)+m if(ml.lt.ij1.or.ml.gt.ij2) goto 21 a=a+atr(1,n,ml)*ant(1,l) + +atr(2,n,ml)*ant(2,l) b=b+atr(1,n,ml)*ant(2,l) + +atr(2,n,ml)*ant(1,l) 21 continue anpt(1,imn) = a anpt(2,imn) = b 20 continue enddo return end c c subroutine xuwigd.f c returns the entire set of Wigner 3jm symbols in double c precision for a given integer group (j1,j2,m1,m2) c needs lnfacd.f to compute ln(z!) c using the algorithm described in [Xu, Journal of Computational c Physics 139, 137-165 (1998)] c subroutine xuwigd(j1,j2,m1,m2,c,cf,n,kmax) include 'gmm01f.par' implicit double precision (a-h,o-z) double precision c(n),cf(n),lnfacd common/fnr/fnr(0:4*(np+1)) data small/1.d-12/ a(j1,j2,j3,m1,m2)=fnr(j3+j1-j2)*fnr(j1+j2+1+j3) + *fnr(j3+m1+m2)*fnr(j3-j1+j2) + *fnr(j1+j2+1-j3)*fnr(j3-m1-m2) b(j1,j2,j3,m1,m2)=dble(2*j3+1)*(dble((m1+m2) + *(j1*j1+j1-j2*j2-j2)) + -dble((m1-m2)*(j3*j3+j3))) do i=1,n c(i)=0.d0 enddo if(iabs(m1).gt.j1.or.iabs(m2).gt.j2) return kmax=j1+j2-max(iabs(j1-j2),iabs(m1+m2))+1 if(kmax.gt.n) stop j3=j1+j2 m3=-m1-m2 s=lnfacd(dble(2*j1))+lnfacd(dble(2*j2)) s=s-lnfacd(dble(2*j1+2*j2+1))+lnfacd(dble(j1+j2-m1-m2)) s=s-lnfacd(dble(j1-m1))-lnfacd(dble(j1+m1)) s=s+lnfacd(dble(j1+j2+m1+m2))-lnfacd(dble(j2-m2)) s=s-lnfacd(dble(j2+m2)) j=j1-j2+m1+m2 c(1)=((-1)**j)*dexp(.5d0*s) if(kmax.eq.1) return c(2)=-b(j1,j2,j3,m1,m2)*c(1) + /(dble(j3+1)*a(j1,j2,j3,m1,m2)) if(kmax.eq.2) return j3f=max(iabs(j1-j2),iabs(m1+m2)) j3=j3f if(j3.eq.0) then c(kmax)=dble((-1)**(j1-m1))/fnr(2*j1+1) else s=lnfacd(dble(j3+j1-j2))+lnfacd(dble(j3-j1+j2)) s=s+lnfacd(dble(j1+j2-j3))+lnfacd(dble(j1-m1)) s=s+lnfacd(dble(j1+m1))+lnfacd(dble(j2-m2)) s=s+lnfacd(dble(j2+m2))-lnfacd(dble(j1+j2+j3+1)) s=s+lnfacd(dble(j3-m3))+lnfacd(dble(j3+m3)) t=0.5d0*s if(j3.eq.j1-j2.or.j3.eq.m3) k=j2+m2 if(j3.eq.j2-j1.or.j3.eq.-m3) k=j1-m1 s=lnfacd(dble(j1+j2-j3-k))+lnfacd(dble(j1-m1-k)) s=s+lnfacd(dble(k))+lnfacd(dble(j2+m2-k)) s=s+lnfacd(dble(j3-j2+m1+k))+lnfacd(dble(j3-j1-m2+k)) j=j1-j2-m3 c(kmax)=((-1)**(k+j))*dexp(t-s) endif if(kmax.eq.3) return if(j3.eq.0) then c(kmax-1)=c(kmax)*dble(m1)/(fnr(j1)*fnr(j1+1)) else c(kmax-1)=-b(j1,j2,j3,m1,m2)*c(kmax)/(dble(j3)* + a(j1,j2,j3+1,m1,m2)) endif if(kmax.eq.4) return cf(kmax)=c(kmax) cf(kmax-1)=c(kmax-1) j3=j3f do i=kmax-2,3,-1 j3=j3+1 cf(i)=dble(j3+1)*a(j1,j2,j3,m1,m2)*cf(i+2) + +b(j1,j2,j3,m1,m2)*cf(i+1) cf(i)=-cf(i)/(dble(j3)*a(j1,j2,j3+1,m1,m2)) enddo j3=j1+j2 do i=3,kmax-2 j3=j3-1 c(i)=dble(j3)*a(j1,j2,j3+1,m1,m2)*c(i-2)+ + b(j1,j2,j3,m1,m2)*c(i-1) c(i)=-c(i)/(dble(j3+1)*a(j1,j2,j3,m1,m2)) cr=dabs(c(i)-cf(i)) if(dabs(c(i)).gt.1.d0) cr=cr/dabs(c(i)) if(cr.lt.small) then do j=i,kmax-2 c(j)=cf(j) enddo return endif enddo return end